A priori analysis for high-fidelity large-eddy simulation of wall-bounded transcritical turbulent flows

The Journal of Supercritical Fluids 207 (2024) 106191 Available online 19 January 2024 0896-8446/© 2024 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect The Journal of Supercritical Fluids journal homepage: www.elsevier.com/locate/supflu A priori analysis for high-fidelity large-eddy simulation of wall-bounded transcritical turbulent flows Marc Bernades∗, Lluís Jofre, Francesco Capuano Department of Fluid Mechanics, Universitat Politècnica de Catalunya ⋅BarcelonaTech (UPC), Barcelona 08034, Spain HIGHLIGHTS •Large-eddy simulation framework for wall-bounded transcritical turbulence. •Kinetic-energyand pressure-equilibriumpreserving scheme. •Scale-similarity-based subfilter-scale closure expressions. •Stable and non-dissipative high-fidelity scale-resolving simulations. GRAPHICAL ABSTRACT ARTICLE INFO Dataset link: https://github.com/marc-bernade s/LES-TT,https://gitlab.com/ProjectRHEA/flo wsolverrhea Keywords: Large-eddy simulation High-pressure Turbulence Supercritical fluids Transcritical wall-bounded flows ABSTRACT Transcritical turbulent flows are governed by the compressible Navier–Stokes equations along with a realgas equation of state. Their computation is strongly susceptible to numerical instabilities and requires kinetic-energyand pressure-equilibrium-preserving schemes to yield stable and non-dissipative scale-resolving simulations. Building upon a recently developed kinetic-energyand pressure-equilibrium-preserving discretization framework based on transporting a pressure equation, the objectives of this paper are to (i) derive a filtered set of equations suitable for large-eddy simulation, and (ii) characterize the properties of the resulting subfilterscale terms by performing a priori analyses of transcritical wall-bounded turbulence direct numerical simulation data. The filtering operation leads to three unconventional subfilter-scale terms that emerge from the pressure equation and require dedicated modeling. The subfilter-scale stress tensor is dissected in terms of magnitude, shape and orientation based on an eigendecomposition analysis, and compared with existing subfilter-scale models. A priori analyses confirm that models of eddy-viscosity type are favorable for this framework, although the tensor shape is not fully captured. Closure expressions are finally proposed and tested for the novel subfilter terms, showing acceptable performances. ∗Corresponding author. E-mail address: [email protected] (M. Bernades). https://doi.org/10.1016/j.supflu.2024.106191 Received 5 November 2023; Received in revised form 16 January 2024; Accepted 17 January 2024 The Journal of Supercritical Fluids 207 (2024) 106191 2 M. Bernades et al. 1. Introduction The study of complex turbulent flows by means of large-eddy simulation (LES) has become increasingly popular in many scientific inquiries as well as in engineering applications. The underlying filtering operation in LES enables to significantly reduce the spatio-temporal resolution requirements compared to direct numerical simulation (DNS). On the other hand, the small-scale motions and their effects on the resolved flow field are not negligible, and therefore require supplementary modeling [1,2]. The development and assessment of novel strategies and models for closing the resulting subfilter-scale (SFS) terms is a very active field of research in a wide range of multiphysics turbulent flow regimes and applications [3–5]. Research efforts are especially challenging in the case of highpressure transor supercritical flows, i.e., when fluids operate across or above their critical point, due to the additional phenomena arising from the large localized thermophysical variations in the vicinity of the pseudo-boiling region [6–8]. Transcritical flows are relevant in many engineering applications, including internal-combustion and rocket engines, among others. Furthermore, the peculiar coexistence of both gas-like and liquid-like states can be fine-tuned for several purposes. For example, they can be leveraged to achieve turbulent regimes in microfluidic devices [9–11], a concept of remarkable interest for energy applications given the enhanced mixing and transfer rates of turbulent flows [9,12]. The difficulties associated with in-vitro characterizations of this regime make high-fidelity numerical simulations an essential tool to elucidate the underlying physics of transand supercritical fluids turbulence; the LES approach is envisioned as to be particularly promising to mitigate the computational complexity of the problem. Only a limited number of studies have focused on LES modeling for trans/supercritical thermodynamic conditions. In one of the earliest works, Bellan [13] analyzed the impact of density gradients on the turbulence mixing processes near critical conditions. The study highlighted the need for appropriate SFS models under supercritical conditions. Selle et al. [14] performed a priori and a posteriori analysis of closure models for homogeneous isotropic turbulence and jets under various thermodynamic regimes. While a priori analyses were encouraging, a posteriori studies provided poor results. They proved that generally neglected subfilter terms in low-pressure cases may become important at high-pressure, hence SFS models need to incorporate the strong thermodynamic non-linearity. Later, Taşkinoğlu and Bellan [15] revisited these results and proposed the approximate deconvolution model as an alternative approach for the heat-flux subfilter term. However, the unavailability of specific SFS models for transcritical flow forced Schmitt et al. [16],Ren et al. [17],Wang et al. [18], among others, to use classical closure models such as Smagorinsky [19] or wall-adapting local eddy-viscosity (WALE) [20,21] models. The impact of classical SFS models on mixing under supercritical pressure for a combustion jet axis was studied by Petit et al. [22]. Three SFS models were tested and evaluated against experimental data: the constant and dynamic Smagorinsky, WALE and Vreman [23] models; very good agreement was found with the dynamic Smagorinsky model. Of note, the pressure-related SFS was modeled by means of the first-order Taylor expansion applied on the filtered field [24]. Analogously, Müller et al. [25] also assessed these SFS models alongside the adaptive local deconvolution method [26] for cyrogenic injection at supercritical pressures. They found that the underlying SFS modeling plays a less important role if one is only interested in the mean flow. In addition, they observed that the relative influence of physical and numerical model uncertainties in the LES of injection at high pressures strongly depends on the thermodynamic regime, i.e., the thermodynamics model is crucial for the prediction of first-order moments under transcritical conditions, whereas its effect is softer at supercritical conditions. Similarly, Borghesi and Bellan [24] analyzed multi-species high-pressure turbulent mixing, and highlighted the need for dedicated LES models. More recently, Unnikrishnan et al. [27] quantified the impact of SFS modeling on the filtered equation of state (EOS) for supercritical turbulent mixing. The direct evaluation of the filtered density (or pressure) based on the Favre-filtered thermodynamic state variables results in computed filtered quantity errors, whose magnitude is biased towards the denser states at the subfilter level. To this extent, two models were proposed, a gradient model and a PDF-based approach. Although the latter was found to be the best candidate, it entails additional information regarding the subgrid variances of temperature and mass fraction, which in turn requires supplementary models. It is clear that more efforts are needed towards a complete and reliable LES formulation for trans/supercritical turbulent flows, especially in the context of wall-bounded configurations. The overwhelming majority of the above-mentioned works is based on a rather classical discretization framework, where the conservation equations (mass, momentum, total energy) are directly evolved, and the corresponding convective terms are discretized in their conservative form. This approach requires numerical dissipation or filtering to stabilize the simulations, which would otherwise suffer from nonlinear instabilities and spurious pressure oscillations [8,28–31]. However, it is widely accepted that LES should be performed with no (or minimal) numerical dissipation, to preserve the subtle energy-cascade processes within the inertial range, where the LES filter cut-off is usually located [32]. In the case of incompressible and ideal-gas compressible flows, stable and non-dissipative simulations can be achieved using kinetic-energy-preserving (KEP) numerical schemes [33], so that dissipation, if needed, is exclusively provided by the SFS models. On the other hand, as studied by Bernades et al. [12], in the case of highpressure transcritical flows, numerical schemes, in addition to being KEP, should also discretely preserve pressure equilibrium, i.e., pressureequilibrium-preserving (PEP), to avoid spurious flow oscillations. Accordingly, a novel scheme that simultaneously satisfies both properties was recently proposed [12]. This framework differs from the numerical methods utilized for previous transcritical flow LES-based analyses. In particular, the PEP property is achieved by solving a pressure evolution equation, while the convective terms in the continuity and momentum equations are expanded according to the Kennedy-Gruber-Pirozzoli (KGP) splitting [33]. As a result, the method (i) preserves kinetic energy by convection, (ii) is locally conservative for mass and momentum, (iii) preserves pressure equilibrium, and (iv) yields stable and robust direct numerical simulations without adding any numerical diffusion to the solution or stabilization procedures. Building upon this framework, this work aims to (i) develop a filtered set of equations suitable for LES based on this novel kinetic-energyand pressure-equilibrium-preserving numerical scheme, and (ii) characterize the properties of the resulting SFS terms by means of performing term-by-term quantitative comparisons and a priori analyses of transcritical wall-bounded turbulence DNS data [8]. Given the novelty of the above-mentioned numerical formulation, in addition to the inherent challenges associated with transcritical LES, an a priori analysis of the pressure-based formulation is highly warranted to understand the relative importance of each of the unclosed terms, with the ultimate objective of deriving physics-based subfilter-scale models. The paper is organized as follows. First, in Section 2, the flow physics modeling and discretization approach for supercritical fluids are presented. Next, the LES framework for high-pressure transcritical turbulence is described in Section 3by introducing the filtering approach, defining the filtered equations of motion with the resulting subfilter terms and presenting the models employed for the study. The DNS of the transcritical channel flow source case and the filtering method applied are presented in Section 4. Section 5analyses the activity of the resolved and subfilter terms, and the subfilter stress tensor is characterized in terms of magnitude, shape and orientation based on an eigendecomposition approach based on the filtered DNS. An a priori analysis for high-fidelity scale-resolving simulations of turbulent flows is presented in Section 6. The presented SFS models are assessed based on correlation analysis and subfilter stress tensor eigendecomposition; the section also proposes closure expressions for the novel unclosed terms. Finally, Section 7reports concluding remarks and proposes future directions. The Journal of Supercritical Fluids 207 (2024) 106191 3 M. Bernades et al. 2. Physical and numerical modeling The turbulent flow motion of supercritical fluids is described by the following set of transport equations of mass, momentum and pressure 𝜕𝜌 𝜕𝑡 + ∇ ⋅(𝜌𝐮)=0,(1) 𝜕(𝜌𝐮) 𝜕𝑡 + ∇ ⋅(𝜌𝐮𝐮)= −∇𝑃+ ∇ ⋅𝝈,(2) 𝜕𝑃 𝜕𝑡 +𝐮⋅∇𝑃+𝜌𝑐2∇⋅𝐮=1 𝜌 𝛽𝑣 𝑐𝑣𝛽𝑇 (𝝈∶ ∇ ⊗𝐮− ∇ ⋅𝒒),(3) where 𝜌is the density, 𝐮is the velocity vector, 𝑃is the pressure, 𝝈=𝜇(∇𝐮+ ∇𝐮𝑇)− (2𝜇∕3)(∇ ⋅𝐮)𝑰is the viscous stress tensor with 𝜇the dynamic viscosity and 𝑰the identity matrix, 𝑐=1∕√𝜌𝛽𝑠is the speed of sound with 𝛽𝑠= −(1∕𝑣)(𝜕𝑣∕𝜕𝑃)𝑠the isentropic compressibility and 𝑣=1∕𝜌the specific volume, 𝛽𝑣= (1∕𝑣)(𝜕𝑣∕𝜕𝑇 )𝑃is the volume expansivity with 𝑇the temperature, 𝑐𝑣is the isochoric specific heat capacity, 𝛽𝑇= −(1∕𝑣)(𝜕𝑣∕𝜕𝑃)𝑇is the isothermal compressibility, and 𝒒= −𝜅∇𝑇is the Fourier heat conduction flux with 𝜅the thermal conductivity. Body forces are not taken into account as the resulting Froude number of the problem is 𝐹𝑟 =𝑢𝑏∕√𝑔𝐻 ≈55, where 𝑢𝑏is the bulk streamwise velocity, 𝑔corresponds to gravitational force and 𝐻is the height of the channel, and consequently inertial forces are roughly 3000×more important; viz. the importance of gravity scales as 1∕𝐹𝑟2[8]. 2.1. Real-gas thermodynamics The thermodynamic space of solutions for the state variables pressure 𝑃, temperature 𝑇, and density 𝜌of a monocomponent substance is described by an equation of state. One popular choice for systems at high pressures, which is used in this study, is the Peng-Robinson [34] equation of state written as 𝑃=𝑅𝑢𝑇 (𝑀∕𝜌) − 𝑏−𝑎 (𝑀∕𝜌)2+2𝑏(𝑀∕𝜌) − 𝑏2,(4) where 𝑅𝑢is the universal gas constant and 𝑀is the molar mass. The coefficients 𝑎and 𝑏take into account real-gas effects related to attractive forces and finite packing volume, respectively, and depend on the critical temperature 𝑇𝑐, critical pressure 𝑃𝑐, and acentric factor 𝜔. They are defined as 𝑎=0.457(Ru𝑇𝑐)2 𝑃𝑐[1+𝑐(1−√𝑇∕𝑇𝑐)]2,(5) 𝑏=0.078Ru𝑇𝑐 𝑃𝑐 ,(6) where coefficient 𝑐is provided by 𝑐={0.380 +1.485𝜔−0.164𝜔2+0.017𝜔3if𝜔 > 0.49, 0.375 +1.542𝜔−0.270𝜔2otherwise.(7) The Peng-Robinson real-gas equation of state needs to be supplemented with the corresponding high-pressure thermodynamic variables based on departure functions [35] calculated as a difference between two states. In particular, their usefulness is to transform thermodynamic variables from ideal-gas conditions (low pressure only temperature dependent) to supercritical conditions (high pressure). The ideal-gas parts are calculated by means of the NASA 7-coefficient polynomial [36], while the analytical departure expressions to high pressures are derived from the Peng-Robinson equation of state as detailed in Jofre and Urzay [7]. 2.2. High-pressure transport coefficients The high pressures involved in the analyses conducted in this work prevent the use of simple relations for the calculation of the dynamic viscosity 𝜇and thermal conductivity 𝜅. In this regard, standard methods for computing these coefficients for Newtonian fluids are based on the correlation expressions proposed by Chung et al. [37,38]. These correlation expressions are mainly function of critical temperature 𝑇𝑐and density 𝜌𝑐, molecular weight 𝑊, acentric factor 𝜔, association factor 𝜅𝑎and dipole moment , and the NASA 7-coefficient polynomial [36]; further details can be found in dedicated works, like for example Poling et al. [39] and Jofre and Urzay [7]. In particular, correlation results of the Chung et al. model [37,38] with respect to the NIST reference database [40] has been carefully validated by Bernades and Jofre [9]. 2.3. Numerical method Simulations of high-pressure transcritical turbulence are strongly susceptible to numerical instabilities due to the presence of non-linear thermodynamic phenomena and large density gradients, which can trigger spurious pressure oscillations that may contaminate the solution and even lead to its divergence. Consequently, it is highly beneficial that the numerical schemes utilized, in addition to being KEP, also attain the so-called PEP property [28,31], which consists in being able to maintain a constant pressure and velocity field when these are initially constant. The numerical scheme utilized in this work has been developed specifically to be simultaneously KEP and PEP. For compressible flow, a family of KEP formulations for the convective term has been recently derived [41]. Instead, the latter property is achieved by solving a pressure evolution equation. A thorough description and validation of this method can be found in Bernades et al. [12,30,42]. In brief, considering a one-dimensional, inviscid, semi-discrete version of Eq. (3) expressed as 𝑃𝑡= −𝑐𝑝− (𝜌𝑐2−𝑃)𝛿𝑥𝑢, (8) where 𝛿𝑥is a discrete centered derivative operator and 𝑐𝑝=𝛿𝑥(𝑃𝑢). The latter double product can be generally discretized as a combination of aconservative and an advective formulation in the form 𝛿𝑥(𝑃𝑢) = 𝜂𝛿𝑥(𝑃 𝑢)+(1−𝜂)(𝑝𝛿𝑥𝑢+𝑢𝛿𝑥𝑃).(9) Obviously, for constant pressure and velocity, any value of 𝜂will lead to 𝑃𝑡=0. The case 𝜂=0 was selected and the corresponding PEP scheme developed by Bernades et al. [12]. The transport equations are numerically solved by adopting a standard semi-discretization procedure; viz. they are first discretized in space and then integrated in time. In particular, spatial operators are treated using second-order central-differencing schemes, and time-advancement is performed by means of a third-order strong-stability preserving (SSP) Runge–Kutta explicit approach [43]. The convective terms are expanded according to the Kennedy-Gruber-Pirozzoli (KGP) splitting [33,41], which has been recently assessed for high-pressure supercritical fluids turbulence [42], yielding non-dissipative stable simulations of compressible discontinuity-free turbulent flows. As a result, the method utilized (i) preserves kinetic energy by convection, (ii) is locally conservative for mass and momentum, (iii) preserves pressure equilibrium, and (iv) yields stable and robust numerical simulations without adding any numerical diffusion to the solution or stabilization procedures. The lack of total energy conservation (TEC) was found not to be detrimental to the results, whereas enforcing the PEP property was of utmost importance in guaranteeing solution fidelity and stability. In fact, the development of methods that are TEC, and simultaneously KEP and PEP, is limited by the lack of a discrete chain rule for linear finitedifferencing schemes [12]. It is worth to mention that a potential alternative framework for strictly subsonic flows would be to employ the low-Mach number approximation, which has been recently utilized in conjunction with real-gas thermodynamics [44]. The discrete conservation and robustness properties of this formulation for realfluid flows deserve a dedicated investigation. Similarly, the overall efficiency of this approach, that involves solving an implicit equation for pressure, should be compared with the supposedly simpler explicit integration of the fully compressible equations. Of note, even if the application studied in this paper is subsonic, the KEP and PEP framework proposed here can be directly employed for mildly compressible or even trans-/supersonic cases, as long as the flow field is free from discontinuities. The Journal of Supercritical Fluids 207 (2024) 106191 4 M. Bernades et al. 3. Large-eddy simulation framework The large-eddy simulation framework is detailed below. First, the filtering approach is presented, then, the filtered equations of fluid motion are described and the related unclosed terms identified for this high-pressure transcritical turbulence numerical method. Next, some classical and state-of-the-art subgrid-scale models available in the literature, and that can be applied in this context, are introduced. 3.1. The filtering approach Following the classical filtering formalism [45] any flow variable 𝑓 is decomposed into a filtered (resolved) contribution 𝑓and a subfilterscale1component 𝑓′, i.e., 𝑓=𝑓+𝑓′. The filtered part 𝑓is defined as 𝑓(𝐱) = ∫𝛺 𝐺(𝐱,𝝃, 𝛥)𝑓(𝐱)𝑑𝝃,(10) where 𝐱and 𝝃are vectors in the flow domain 𝛺. The filter function 𝐺 usually depends on a parameter 𝛥, called the filter width, and satisfies the normalization condition ∫𝛺 𝐺(𝐱,𝝃, 𝛥)𝑑𝝃=1,(11) for every 𝐱and 𝛺. For compressible flows, Favre [46] introduced a related filter operation  𝑓=𝜌𝑓 𝜌,(12) which leads to the decomposition 𝑓= 𝑓+𝑓′′, where 𝑓′′ corresponds to the residual field. Typical filters commonly used in large-eddy simulation include the top-hat (also referred as spatial or box filter), Gaussian and spectral cutoff filter [47]. In three dimensions, the filter width is usually defined as 𝛥= (𝛥1𝛥2𝛥3)1∕3,(13) where the symbol 𝛥𝑖denotes the filter width along the 𝑖th direction. For homogeneous filters, i.e., for kernels that do not depend on the spatial position 𝐱, the filtering operator commutes with spatial derivatives [48, 49]. On the other hand, the Favre filter in general does not commute with partial derivatives [47]. 3.2. Filtered equations of fluid motion LES formulations based on a real-gas EOS are inherently non-trivial due to the non-linearity of the thermodynamic relations at play. In this work, inspired by recent efforts [27] and based on the proposed novel scheme [12], the LES equations based on the classical explicit lowpass filtering formalism are derived and the corresponding unclosed terms identified. Interestingly, previously unknown terms arise from the analysis, which will require specific modeling efforts. In this regard, it is also worth to emphasize that transcritical wall-bounded flow physics profoundly differs from that observed in subcritical cases, most notably due to the presence of a strong baroclinic instability and, as a consequence, to the failure of standard scaling transformations for velocity within the pseudo-boiling region. The LES equations are written for the transported variable vector Ψ=[𝜌, 𝜌𝐮, 𝑃 ], which in terms of Favre-filtered variables writes as Ψ= [𝜌, 𝜌 𝐮, 𝑃 ]. Assuming that differentiation and filtering commute [50,51], i.e, the filtered Navier–Stokes equations are valid if the filtered operator is any linear operator that commutes with partial differential operators 1In this work, the term subfilter is deliberately used in place of the (perhaps) more popular subgrid, to emphasize that a clear distinction between the filtering operation, the SFS modeling and the discretization is assumed. 𝜕𝑡 and 𝜕𝑗 [47], the LES equations describing the motion of supercritical fluid turbulence correspond to the following set of low-pass filtered equations 𝜕𝜌 𝜕𝑡 + ∇ ⋅(𝜌 𝐮)=0,(14) 𝜕(𝜌 𝐮) 𝜕𝑡 + ∇ ⋅(𝜌 𝐮 𝐮)+ ∇𝑃− ∇ ⋅ 𝝈= −𝛼1+𝛼2,(15) 𝜕𝑃 𝜕𝑡 + 𝐮⋅∇𝑃+𝜌 𝑐2∇⋅ 𝐮−1 𝜌  𝛽𝑣 𝑐𝑣 𝛽𝑇 ( 𝝈∶ ∇ ⊗ 𝐮− ∇ ⋅ 𝒒) = 𝛼3+𝛼4+𝛼5, (16) where 𝜌and 𝑃are the filtered density and pressure variables, respectively, and  𝐮corresponds to the Favre-filtered velocity vector. As usual in LES, the Eqs. (14)–(16) are presented such that the left-hand sides resemble the Navier–Stokes Eqs. (1)–(3) expressed in the filtered variables Ψ, whereas the right-hand sides of Eqs. (14)–(16) contains subfilter-scale stresses 𝛼𝑖, which read 𝛼1= ∇ ⋅𝜌𝝉,(17a) 𝛼2= ∇ ⋅(𝝈−𝜎),(17b) 𝛼3= ( 𝐮∇⋅𝑃−𝐮∇⋅𝑃),(17c) 𝛼4= (𝜌 𝑐2∇⋅ 𝐮−𝜌𝑐2∇⋅𝐮),(17d) 𝛼5=[1 𝜌 𝛽𝑣 𝑐𝑣𝛽𝑇 (𝝈∶ ∇ ⊗𝐮− ∇ ⋅𝒒) − 1 𝜌  𝛽𝑣 𝑐𝑣 𝛽𝑇 ( 𝝈∶ ∇ ⊗ 𝐮− ∇ ⋅ 𝒒)],(17e) where the variables denoted by ⋅correspond to the field function of filtered transported variables, i.e., 𝑓(Ψ). Mass equation. By definition the Favre-filtered mass Eq. (14) does not generate any SFS term, as the filtered product 𝜌𝐮is equivalent to 𝜌 𝐮. Momentum equation. The filtered momentum equation, Eq. (15), produces two unclosed terms. The first one, 𝛼1, generates due to the non-linearity of the convective term. This term is written below in the form of the so-called SFS stress tensor 𝝉 𝜌𝝉=𝜌( 𝐮𝐮 − 𝐮 𝐮),(18) which corresponds to the interaction between subfilter and resolved scales, and consequently its closure requires modeling. Second, the resulting subfilter-term from the viscous stress tensor 𝛼2, which results from the non-linearity of the viscous term and the fact that Favre-filtering operator does not commute with partial derivatives. It is typically neglected in high-Reynolds-number flows at low-pressure conditions. In fact, a priori tests confirm that it is an order of magnitude smaller than 𝛼1in subcritical conditions, where  𝝈=𝑓( 𝐮, 𝑇)[47]. Nevertheless, due to the strong non-linearities that arise in transcritical regimes in the vicinity of the pseudo-boiling line [9], this term requires closure to model the difference with respect to 𝝈=𝑓(𝐮, 𝑇 )[24]. Pressure equation. The filtered pressure equation was initially analyzed by Zang et al. [52] in the context of ideal-gas thermodynamics. In the case of real fluids, the resulting subfilter terms are (i) 𝛼3, which represents the effect of subfilter-scale turbulence on the power of pressure forces, (ii) 𝛼4is associated with flow dilatation due to compressibility effects, with the speed of sound depending on the filtered variables as 𝑐 =𝑓(𝜌, 𝑃 ), and (iii) 𝛼5contains the subfilter contribution of the viscous, Fourier  𝒒=𝑓(𝜌,  𝑇)and thermophysical quantities where  𝛽𝑣= 𝑓(𝜌, 𝑃 ), 𝛽𝑇=𝑓(𝜌, 𝑃 )and 𝑐𝑣=𝑓(𝜌, 𝑃,  𝑇). It is worth noticing that the term directly associated with the Fourier flux that appears when total energy is evolved is embedded into 𝛼5. Equation of state. The Peng-Robinson equation relates the thermodynamic variables and can be generally expressed as 𝑃=𝑃(𝜌, 𝑇 ). Hence, the filtered equation of state reads 𝑃=𝑃(𝜌, 𝑇 ) = 𝑅𝑢𝑇 (𝑀∕𝜌) − 𝑏−𝑎 (𝑀∕𝜌)2+2𝑏(𝑀∕𝜌) − 𝑏2.(19) The Journal of Supercritical Fluids 207 (2024) 106191 5 M. Bernades et al. However, this equation is not expressed as a function of the transported (filtered) variables. Hence, similarly to the filtered Navier–Stokes equations, the EOS can be rearranged based on the filtered variables, i.e., 𝑃(𝜌,  𝑇), as follows 𝑃=𝑃(𝜌,  𝑇) + 𝛼6,(20) where in this case the temperature  𝑇=𝑓(𝑃 , 𝜌)is expressed as a function of the filtered density and pressure from the transported variables vector Ψ, with the filtered coefficient related to attractive forces defined as 𝑎 =𝑓( 𝑇). As a consequence, Eq. (20) generates an additional unclosed term, 𝛼6, which can be expressed as 𝛼6=𝑃(𝜌, 𝑇 ) − 𝑃(𝜌,  𝑇).(21) 3.3. Subfilter-scale models The SFS terms appearing in the right-hand side of Eqs. (14)−(16) contain information from the unfiltered field and, thus, they require models in order to express the filtered Navier–Stokes equations in filtered variables. The novelty of the proposed pressure-based framework induces SFS terms which have not yet been characterized. In this regard, this section briefly introduces several classical subfilter-scale models typically used for super/transcritical simulations based on a careful literature survey. In particular, SFS models have been proposed for (i) the SFS stress tensor 𝝉(𝛼1), (ii) the trace of the SFS stress tensor models, and (iii) the real-gas equation of state (𝛼6). 3.3.1. Subfilter stress tensor The general objective of any SFS model is to express unclosed terms as a function of known (transported) variables, i.e., filtered variables. In the case of the residual stress tensor, this is achieved by substituting the unknown value of 𝜏𝑖𝑗 by a model 𝜏𝑆𝐹 𝑆 𝑖𝑗 . A popular and successful modeling approach is based on the eddy-viscosity concept, where the anisotropic part of the tensor 𝜏𝑎 𝑖𝑗 𝑆𝐹 𝑆 is modeled as 𝜏𝑎 𝑖𝑗 𝑆𝐹 𝑆 =𝜏𝑆𝐹 𝑆 𝑖𝑗 −1 3𝛿𝑖𝑗 𝜏𝑆𝐹 𝑆 𝑘𝑘 = −𝜈𝑆𝐹 𝑆 𝑆𝑖𝑗(𝑢),(22) where 𝜈𝑆𝐹 𝑆 is the SFS eddy-viscosity, 𝛿𝑖𝑗 is the Kronecker symbol, and 𝑆𝑖𝑗(𝑢)=(𝜕𝑢𝑖∕𝜕𝑥𝑗+𝜕𝑢𝑗∕𝜕𝑥𝑖)−2∕3𝛿𝑖𝑗 𝜕𝑢𝑘∕𝜕𝑥𝑘is the rate-of-strain tensor. As it can be observed from Eq. (22),𝜈𝑆𝐹 𝑆 governs the magnitude of the tensor, while instead, its degree of anisotropy and orientation are determined by 𝑆𝑖𝑗(𝑢), hereinafter expressed also as  𝑆𝑖𝑗. Typically, in incompressible flows, the isotropic part of the SFS tensor is not modeled and incorporated into the filtered pressure. On the other hand, closure models for the trace tensor are required to obtain the isotropic term in compressible flows. The most popular choices are the Yoshizawa [53] and Vreman et al. [54] models also presented in this section. Classical smagorinsky. The constant-coefficient Smagorinsky was one of the first developments based on an eddy-viscosity model [19,55] and is expressed as 𝜈𝑆𝐹 𝑆 = (𝐶𝑆𝛥)2|𝑆(𝑢)|,(23) where |𝑆(𝑢)|= (1∕2 𝑆𝑖𝑗  𝑆𝑖𝑗)1∕2, hereafter written also as | 𝑆|, and 𝐶𝑆 is the Smagorinsky constant whose values are typically 0.2 in isotropic turbulence and 0.1 for turbulent channel flow [47]; the latter value is selected for this work. Dynamic smagorinsky (𝜏𝑖𝑗𝑎𝑆𝐹𝑆1). In the context of eddy-viscosity models, the dynamic approach determines the suitable local value of the model coefficients by comparing the amount of eddy dissipation needed at different scales. Dynamic models are fundamentally based on the Germano identity [56], which relates the SFS stress tensor at two filter levels: (i) the one associated with 𝛥, in this context denoted as Flevel, and (ii) an additional one denoted as G-level and expressed by  (.), which usually corresponds to a filter width of 2𝛥. The consecutive applications of these two filters results in  (.), a so-called ‘‘FG-level’’ filter. The subgrid-term on FG-level is expressed as 𝑇𝑓= 𝑓(𝑥) − 𝑓( 𝑥),(24) where 𝑥is a vector function of time and space. This expression can be related to the F-level as 𝑇𝑓−𝜏𝑓=𝐿𝑓,(25) resulting in the generalized Germano identity, where 𝜏𝑓=𝑓(𝑥) − 𝑓(𝑥) and 𝐿𝑓= 𝑓(𝑥) − 𝑓( 𝑥). The left-hand side term of Eq. (25) cannot be calculated from variables on the F-level. In particular, for the subfilter stress tensor, where 𝑥=𝑥(𝜌, 𝐮)and then 𝑓(𝑥) = 𝜌𝑢𝑖𝑢𝑗, the identity reads  𝜌𝑇𝑖𝑗 − 𝜌𝜏𝑎 𝑖𝑗 =𝐿𝑖𝑗,(26) with  𝜌𝑇𝑖𝑗 and 𝐿𝑖𝑗 yielding  𝜌𝑇𝑖𝑗 = 𝜌𝑢𝑖𝑢𝑗− 𝜌𝑢𝑖 𝜌𝑢𝑗∕ 𝜌, 𝐿𝑖𝑗 = 𝜌𝑢𝑖𝜌𝑢𝑗∕𝜌− 𝜌𝑢𝑖 𝜌𝑢𝑗∕ 𝜌. (27) To this extent, adopting Smagorinsky formulae and replacing 𝐶𝑆 constant by coefficient 𝐶𝐷[56] the dynamic eddy-viscosity model results in 𝜈𝑆𝐹 𝑆 = (𝐶𝐷𝛥)2|𝑆(𝑢)|,(28) This SFS model 𝜏𝑎 𝑖𝑗 is substituted into the Germano identity. According to Germano’s notation, the quantities 𝑇𝑖𝑗 and 𝜏𝑎 𝑖𝑗 are obtained by formulating the subfilter model in FGand F-filtered quantities yielding ⎧ ⎪ ⎨ ⎪ ⎩ 𝐶𝐷𝑀𝑖𝑗 =𝐿𝑖𝑗, 𝑀𝑖𝑗 = − 𝜌(2𝛥)2|𝑆(𝐯)|𝑆𝑖𝑗(𝐯) +  𝜌𝛥2|𝑆( 𝐮)|𝑆𝑖𝑗( 𝐮), (29) where 𝐯= 𝜌𝐮∕ 𝜌represents the velocity at FG-level. Therefore, 𝐶𝐷 can be obtained with a system of six-equations solved as a least-square approach. However, artificial modification is required if it returns negative values to prevent numerical instability 𝐶𝐷=⟨𝑀𝑖𝑗𝐿𝑖𝑗⟩ ⟨𝑀𝑖𝑗𝑀𝑖𝑗 ⟩.(30) Alternative dynamic mixed model variations also exist such as the hybrid-based model proposed by Zang et al. [57] or the Clark model [58]. Anisotropic minimum-dissipation model (AMD, 𝜏𝑖𝑗𝑎𝑆𝐹 𝑆2). This category of subfilter models provides the minimum eddy dissipation required to dissipate the energy of the subfilter scales. The first model of this family was the QR model for isotropic grids [59,60]. However, a lack of consistency with the exact subfilter tensor was observed. In fact, this method relies solely on the filter width and cannot be extended to anisotropic grids. In particular, the variations along the directions with smallest grid spacings dominate, and as a result the model leads to insufficient damping in the largest grid spacing direction. Consequently, Rozema et al. [61] developed the anisotropic minimum-dissipation model (AMD) for incompressible flow, which was later extended by Abkar and Moin [62] to model the subfilter scalar flux. In this case, the energy of the subfilter scales is regularized using a modified Poincaré inequality and results in a modified constant which also incorporates the dependence on the filter width by scaling the velocity gradient. Instead, the production and dissipation are expressed in terms of the velocity gradient. The eddy-viscosity model can be expressed using the tensor invariants [63] as 𝜈𝑆𝐹 𝑆 = (𝐶 𝛥)2𝑚𝑎𝑥 [0,−(𝐼3−𝐼4)] 𝐼1−𝐼2 ,(31) The Journal of Supercritical Fluids 207 (2024) 106191 6 M. Bernades et al. 𝐶𝑤2=𝐶𝑠2⟨√2(1∕2 𝑆𝑖𝑗1∕2 𝑆𝑖𝑗)3∕2⟩ ⟨1∕2 𝑆𝑖𝑗1∕2 𝑆𝑖𝑗(1∕2 𝑆𝑑𝑖𝑗1∕2 𝑆𝑑𝑖𝑗)3∕2∕[(1∕2 𝑆𝑖𝑗1∕2 𝑆𝑖𝑗)5∕2+ (1∕2 𝑆𝑑𝑖𝑗1∕2 𝑆𝑑𝑖𝑗)5∕4]⟩ .(34) Box I. where 𝐶=0.3 for second-order finite difference numerical schemes. The invariants read 𝐼1=𝑡𝑟 [(1∕2 𝑆𝑖𝑗)2], 𝐼2=𝑡𝑟 [(1∕2 𝛺𝑖𝑗)2], 𝐼3=𝑡𝑟 [(1∕2 𝑆𝑖𝑗)3], 𝐼4=𝑡𝑟 [(1∕2 𝑆𝑖𝑗)(1∕2 𝛺𝑖𝑗)2], (32) where  𝛺𝑖𝑗 =𝜕𝑢𝑖∕𝜕𝑥𝑗−𝜕𝑢𝑗∕𝜕𝑥𝑖is the rate-of-rotation tensor. Of note, the rate-of-strain tensor employed follows the compressible-flow convention, defined in Eq. (22), where the term −2∕3∇(𝐮)is subtracted on the tensor diagonal. Nevertheless, it will be shown in the results that its weight represents (10−3) with respect to off-diagonal terms. Wall-adapting local eddy-viscosity model (WALE, 𝜏𝑖𝑗𝑎𝑆𝐹𝑆3). The WALE model [20,21] has been previously used under supercritical conditions for a posteriori LES [22,64]. It is an eddy-viscosity-based model that takes into account the local strain and rotation rates. It has a number of favorable properties: it correctly recovers the asymptotic scaling at the wall, and provides low dissipation in pure shear zones, avoiding the need for damping functions in laminar or transitional areas of the flow. The WALE model uses the velocity gradient  𝐴𝑖𝑗 =𝜕𝑢𝑖∕𝜕𝑥𝑗to represent the fluctuations at the length scale 𝛥to overcome the drawbacks of the Smagorinsky model, which is solely based on the second invariant of the symmetric part of 𝑆𝑖𝑗. The WALE model reads 𝜈𝑆𝐹 𝑆 = (𝐶𝑤𝛥)2(1∕2 𝑆𝑑𝑖𝑗 1∕2 𝑆𝑑𝑖𝑗)3∕2 (1∕2 𝑆𝑖𝑗 1∕2 𝑆𝑖𝑗)5∕2+ (1∕2 𝑆𝑑𝑖𝑗 1∕2 𝑆𝑑𝑖𝑗)5∕4,(33) where  𝑆𝑑𝑖𝑗 = (  𝐴2 𝑖𝑗 + 𝐴2 𝑗𝑖) −2∕3𝛿𝑖𝑗  𝐴2 𝑘𝑘 is the trace-less symmetric part of the square of the Favre-filtered velocity gradient tensor. The constant 𝐶𝑤is derived by assuming that the new model gives the same ensembleaverage subfilter kinetic energy dissipation as the classical Smagorinsky as given in Box I. Similarly, Rieth et al. [65] reported that 𝜎−model also accurately captures the near wall behavior with respect to dynamic Smagorinsky. Scale-similarity model (𝜏𝑖𝑗𝑎𝑆𝐹 𝑆4). Unlike the models presented above, this is not of eddy-viscosity type; rather, it assumes similarity among the energy transfer that occurs between flow structures at different scales [66,67]. More specifically, the tensor model is defined as a function of only filtered variables 𝜏𝑎 𝑖𝑗 𝑆𝐹 𝑆 =1 𝜌(𝜌𝑢𝑖𝜌𝑢𝑗∕𝜌−𝜌𝑢𝑖𝜌𝑢𝑗∕𝜌)=1 𝜌(𝜌𝑢𝑖𝑢𝑗−𝜌𝑢𝑖𝜌 𝑢𝑗∕𝜌).(35) It is known that the correlation to the exact SFS stress is relatively high, hence the model predicts structures of the stresses at the right locations. Nevertheless, the magnitude of the turbulent stresses is less accurately predicted. Due to the filtered-variables-based model, it typically underestimates the turbulent stresses in the turbulent regime, whereas in laminar regimes the model is not too dissipative. 3.3.2. Models for the trace of the SFS stress tensor In the case of compressible flows, explicit subfilter models are required for 𝜏𝑆𝐹 𝑆 𝑘𝑘 , unlike in incompressible flows where the isotropic part 𝜏𝑆𝐹 𝑆 𝑘𝑘 ∕3 is added to the filtered pressure [68]. Two models are presented in this regard. Yoshizawa model (𝜏𝑘𝑘𝑆𝐹 𝑆1). One choice is the parametrization by Yoshizawa [53] 𝜏𝑆𝐹 𝑆 𝑘𝑘 =𝐶𝐼𝛥2| 𝑆|2,(36) where 𝐶𝐼is a model coefficient that can be approximated, for wallbounded flows, as ensemble-average scalars for each direction, as proposed by Moin et al. [69] 𝐶𝐼=⟨ 𝜌𝑢𝑖𝑢𝑖− (1∕ 𝜌)(  𝜌𝑢𝑖)(  𝜌𝑢𝑖)⟩ ⟨ 𝜌 𝛥2|  𝑆|2 −𝛥2 𝜌| 𝑆|2⟩ .(37) Vreman model (𝜏𝑘𝑘𝑆𝐹 𝑆2). A different approach by Vreman et al. [54] proposed to model 𝜏𝑆𝐹 𝑆 𝑘𝑘 as 𝜏𝑆𝐹 𝑆 𝑘𝑘 =𝐶𝑤𝛥2∑ 𝑖,𝑗  𝑆2 𝑖𝑗,(38) where 𝐶𝑤=0.03, based on the filtered DNS of the problem of interest, although for homogeneous isotropic turbulence 𝐶𝑤=0.325 is proposed [70]. 3.3.3. Equation of state unclosed term model The subfilter term 𝛼6that comes from the real EOS cannot be neglected at transcritical conditions due to the strong non-linearity effects in the vicinity of the pseudo-boiling line [24,71]. This is in contrast with so-called low-pressure standard LES assumptions (SLA) [24], where the temperature is directly expressed as a function of the Favrefiltered variables, i.e., 𝑃(𝜌, 𝑇 ) = 𝑃(𝜌,  𝑇). Taylor expansion model. Selle et al. [14] performed an a priori analysis of supercritical multi-species mixing layers. It was first checked that the SLA is not suitable for such application, i.e., 𝑓=𝑓(𝛹). In addition, they evaluated the relative importance of the various unclosed SFS terms, and the term ∇[𝑃−𝑃(𝛹)]yielded an important contribution in magnitude. To this extent, they proposed a model based on the expansion of the filtered EOS in Taylor series as 𝑃(𝛹) = 𝑃(𝛹) + 𝜕𝑃 𝜕𝛹𝑚||||𝛹=𝛹 (𝛹𝑚−𝛹𝑚) +1 2 𝜕2𝑃 𝜕𝛹𝑚𝜕𝛹𝑛||||𝛹=𝛹 (𝛹𝑚−𝛹𝑚)(𝛹𝑛−𝛹𝑛) + (𝛹3),(39) which assuming that (i) 𝑃(𝛹) = 𝑃(𝛹), (ii) the partial derivatives terms of the expansion can be removed from filtering operation, and (iii) filter is a projection implying that 𝛹𝑚−𝛹𝑚=0, it then reads ⎧ ⎪ ⎨ ⎪ ⎩ 𝑃(𝛹) = 𝑃(𝛹) + 𝛿𝑃, 𝛿𝑃=1 2 𝜕2𝑃 𝜕𝛹𝑚𝜕𝛹𝑛||||𝛹=𝛹 (𝛹𝑚𝛹𝑛−𝛹𝑚𝛹𝑛),(40) which 𝛿𝑃is a second-order approximation. This addition of 𝛿𝑃was found to potentially deteriorate the correlation between the filtered pressure and the LES results, and consequently is only valid for small perturbations, i.e., 𝛹−𝛹←←→ 0. Improved LES assumptions (ILA). Borghesi and Bellan [24] performed an a priori analysis for a mixing layer with real-gas and multiple species, and proposed a so-called improved LES assumptions (ILA) model. This model relies on a combination of the scale-similarity assumption and The Journal of Supercritical Fluids 207 (2024) 106191 7 M. Bernades et al. the Germano identity for 𝑃(𝜌, 𝑇 ), the Fourier term and species-related unclosed terms. In this regard, for the context of this work, only focus is placed on the former unclosed term, denoted as 𝛼6in Section 3.2. Therefore, without loss of generality this G-level SFS flux can be defined as function of the transported fields 𝛹as 𝜗(𝛹) = 𝑃(𝛹) − 𝑃(𝛹), 𝛶 (𝛹) = 𝜗(𝛹) =  𝑃(𝛹) − 𝑃( 𝛹).(41) Following Germano identity, they can be recast as 𝐿𝑝=𝛶− 𝜗=𝜗(𝛹) −  𝜗(𝛹) =  𝑃(𝛹) − 𝑃( 𝛹),(42) which is used to evaluate the scale-similarity coefficient and within this context it can be modeled as 𝜗(𝛹) = 𝐶′[𝑃(𝛹) − 𝑃(𝛹)], 𝛶 (𝛹) = 𝐶′′ ⎡⎢⎢⎣  𝑃( 𝛹) − 𝑃(   𝛹)⎤⎥⎥⎦ .(43) Assuming 𝑐𝑝=𝐶′=𝐶′′ the Germano identity can be rewritten as ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝐿𝑝=𝑐𝑝𝑀𝑝, 𝑀𝑝=⎡⎢⎢⎣  𝑃( 𝛹) − 𝑃(   𝛹)⎤⎥⎥⎦ −[ 𝑃(𝛹) −  𝑃(𝛹)].(44) As previously introduced, the 𝑐𝑝can be found by computing the least-square method resulting in 𝑐𝑝=⟨𝑀𝑝𝐿𝑝⟩ ⟨𝑀𝑝𝑀𝑝⟩.(45) As a result, the model counterpart of the exact filtered term 𝑃(𝛹)is calculated as 𝑃(𝛹) = 𝑃(𝛹) + 𝑐𝑝(𝑃(𝛹) − 𝑃(𝛹)).(46) 4. Filtered DNS dataset Before tackling the problem from a modeling standpoint, a dataset based on direct numerical simulations of a transcritical channel flow is constructed for the subsequent analyses. In this regard, this section describes the flow setup, the filter selection and the consequent filtered datasets that will be employed for the resolved and SFS terms assessment and the a priori analysis. 4.1. DNS setup The DNS of transcritical channel flow computed by Bernades et al. [8] utilizing the in-house flow solver RHEA [72,73] is briefly described in this section. The system operates with N2at a supercritical bulk pressure of 𝑃𝑏∕𝑃𝑐=2 and confined between bottom/cold (𝑏𝑤) and top/hot (𝑡𝑤) isothermal walls, separated in this case at a distance 𝐻= 2𝛿with 𝛿=100 µm the channel half-height, at 𝑇𝑐𝑤∕𝑇𝑐=0.75 and 𝑇ℎ𝑤∕𝑇𝑐=1.5, respectively, where sub-indexes 𝑏and 𝑐correspond to bulk and critical point values. In this regard, a schematic of the problem of interest is displayed in Fig. 1(a). This configuration forces the fluid to undergo a transcritical trajectory by operating within a thermodynamic region across the pseudo-boiling line as illustrated on the 𝑃−𝑣diagram of Fig. 1(b). The friction Reynolds number (denoted by sub-index 𝜏) selected at the bottom wall is Re𝜏,𝑐𝑤 =𝜌𝑐𝑤𝑢𝜏,𝑐𝑤𝛿∕𝜇𝑐𝑤 =100 to ensure fully-developed turbulent flow conditions. The corresponding dimensional parameters are: dynamic viscosity 𝜇𝑐𝑤 =1.6⋅10−4Pa s, density 𝜌𝑐𝑤 =839.4 kg/m3, and friction velocity 𝑢𝜏,𝑐𝑤 =1.9⋅10−1m/s. The boundary conditions imposed correspond to an impermeable no-slip condition for the wall normal-direction, and periodic for the streamwise and spanwise directions. The computational domain is 4𝜋𝛿 ×2𝛿×4∕3𝜋𝛿 in the streamwise (𝑥), wall-normal (𝑦), and spanwise (𝑧) directions, respectively. The grid is uniform in the streamwise and spanwise directions with resolutions in wall units (based on 𝑐𝑤 values) equal to 𝛥𝑥+=9.8 and 𝛥𝑧+=3.3, and stretched towards the walls in the vertical direction with the first grid point at 𝑦+=𝑦𝑢𝜏,𝑐𝑤∕𝜈𝑐𝑤 =0.1 and with sizes in the range 0.4≤𝛥𝑦+≤2.3. Thus, this arrangement corresponds to a grid size of 128 ×128 ×128 points. Based on the estimates provided by Jofre and Urzay [7], the characteristic length scale for density gradients in this case is approximately 10×larger than the Kolmogorov scale, therefore the latter is the driving factor to select mesh resolution. The selected grid size is thus assumed to resolve all the relevant flow scales. The simulation strategy starts from a linear velocity profile with random fluctuations, which is advanced in time to reach turbulent steady-state conditions after approximately five flowthrough-time (FTT) units; based on the bulk velocity 𝑢𝑏and the length of the channel 𝐿𝑥=4𝜋𝛿, a FTT is defined as 𝑡𝑏=𝐿𝑥∕𝑢𝑏∼𝛿∕𝑢𝜏. In this regard, flow statistics are collected for roughly 10 FTTs once steady-state conditions are achieved. The numerical method described in Section 2.3 is used to solve the governing equations. 4.2. Filtered DNS high-pressure supercritical turbulence Filter selection. The top-hat filter has often been used in the context of a priori analysis and specifically for high-pressure turbulence of real fluids [14,24,71]. The filter can be expressed based on Taylor series expansion and directly applied to a generic discrete scalar field 𝜃as follows [74] 𝜃=𝜃+(𝜖𝛥)2 24 ∇2𝜃+(𝛥4),(47) where 𝛥is the characteristic cell size, 𝜖is the dimensionless parameter controlling the filter width (𝛥=𝜖𝛥). However, as highlighted in recent works [12,28,30], filters of this type may induce spurious oscillations in the vicinity of strongly non-linear regions, viz. across the pseudo-boiling line. Although the DNS dataset was completely free of oscillations thanks to the KEP and PEP numerical method utilized, the application of the top-hat filter resulted indeed in an increase of turbulent kinetic energy levels, which could not be mitigated by increasing the spatial order of accuracy of the filter. A family of local extrema diminishing (LED) filters is therefore finally selected to prevent overshoots in the filtered quantities. Particularly, the adaptive box filter is employed, where LED is achieved for any 𝜖≥1. This discrete filter kernel consists of an averaging of the values of a neighborhood stencil as [75] 𝜃𝑜=1 𝜖(𝜃𝑜+𝜖−1 ∑𝑝∈𝑁𝑜𝜔𝑝∑ 𝑝∈𝑁𝑜 𝜃𝑝𝜔𝑝),(48) where subscript (⋅)𝑜and (⋅)𝑝denote to current filtered position and the corresponding neighbors, respectively, and 𝜔𝑝is the volume of each of the neighbor grid cells. Of note, the width of the filter stencil changes accordingly with 𝜖. Filtered transcritical DNS dataset. The effect of the filter on the small scales is depicted in Fig. 2. The ensemble-averaged TKE (averaged along the spanwise and streamwise directions) at different values of 𝑦+(denoted as ⟨⋅⟩), normalized by the DNS value (labeled as ‖⋅‖), decreases monotonically with the filter width for both bottom/cold (cw) and top/hot walls (hw). For the hw, 𝜖=2 ensures that at least 90% of the TKE is captured. Fig. 3 qualitatively compares the DNS contours of 𝑢+with filtered DNS at  𝛥∕𝛥=2 and  𝛥∕𝛥=4. The loss of resolution on the small scales of the flow can be clearly seen, which becomes more pronounced for larger filter widths. In this case,  𝛥∕𝛥=2 holds 90% of TKE, whereas  𝛥∕𝛥=4 produces much stronger effects and will be hereinafter labeled as coarse LES, although this definition should be intended in relative rather than absolute terms. The effect of larger filter widths and larger Reynolds numbers will be explored as part of future work. For completeness, the filtered behavior on high-pressure is compared against equivalent energized flow case at low-pressure with the same volumetric input power as the high-pressure setup [8]. Fig. 4 presents the firstand second-order statistics for both DNS The Journal of Supercritical Fluids 207 (2024) 106191 8 M. Bernades et al. Fig. 1. (a) Schematic of the channel flow problem studied. (b) Transcritical trajectory (labeled with horizontal double-end arrow at system bulk pressure 𝑃𝑏=2𝑃𝑐) of the system represented on a 𝑃−𝑣diagram for N2. The auxiliary curves correspond to: vapor–liquid equilibrium (VLE) line, pseudo-boling line (Pb), critical point (CP), and isothermal trajectories at various temperatures to properly identify the fluid states, i.e., supercritical liquid-like, supercritical gas-like, liquid, vapor and 2-phase region inside the saturation curve. Fig. 2. Normalized turbulent kinetic energy vs. filter width  𝛥∕𝛥at different 𝑦+levels for (a) cold and (b) hot wall. cases at different filter levels. Two different observations can be made from these plots. First, the filtering behavior is similar on the firstorder statistics in both lowand high-pressure frameworks. Second, the second-order fluctuations are more sensitive to filter width for the transcritical case along the wall-normal direction, and in particular in the near-wall regions. This evidence reflects (i) the higher filter sensitivity due to the strong thermodynamic non-linearities, and (ii) the enhanced levels of fluctuations of the high-pressure case. The Journal of Supercritical Fluids 207 (2024) 106191 9 M. Bernades et al. Fig. 3. Snapshot of instantaneous streamwise velocity in wall units 𝑢+(cw) on a 𝑥-𝑦slice for (a) DNS and (b) filtered DNS with adaptive box filter with  𝛥∕𝛥=2 and (c)  𝛥∕𝛥=4. 5. Resolved and SFS term-by-term analysis Before conducting an a priori analysis of several existing SFS models and proposing relevant closure expressions, the various unclosed terms identified in Section 4.2 are preliminarily analyzed in terms of their relative activity, i.e., the relative magnitude of resolved and SFS contributions, and the behavior of the SFS stress tensor is dissected in terms of magnitude, shape and orientation based on an eigendecomposition analysis. Of note, the study focuses only on the outer region of the channel flow, mainly to avoid near-wall filtering and misleading conclusions from unresolved scales in the viscous layer, which are smaller near the edge of the buffer layer and grow linearly from the wall [76]. Hence, it is assumed that near-wall scales solely rely either on the wall-resolved or the wall-model LES. The grid resolution is typically imposed with a coarser mesh on the viscous layer and finer grid-spacing in the outer layer at wall height of a fraction of the boundary-layer thickness (𝛿𝐵𝐿) as 𝑦∕𝛿𝐵𝐿 ≈0.2 or 𝑦+≈50 [77], and as a result approximately 80% of the boundary layer is treated by LES. Consequently, in this DNS, the offwall region at 𝑦∕𝛿𝐵𝐿 ≥0.2 corresponds to 𝑦∕𝛿=0.2, and equivalently in wall units, 𝑦+ 𝑐𝑤 =19.3 and 𝑦+ ℎ𝑤 =37.4, where the boundary layer thickness reads 𝛿𝐵𝐿 =Re𝜏,𝑐𝑤 𝛿𝑣=10−4with viscous length scale 𝛿𝑣=𝜇𝑐𝑤∕(𝜌𝑐𝑤 𝑢𝜏,𝑐𝑤) = 10−6and Re𝜏,𝑐𝑤 =100. 5.1. Activity of the terms As a first step, an assessment of the relative importance of each of the unclosed SFS term defined in Section 4.2 (i.e., 𝛼𝑖for 𝑖∈ [1,…,6]) is performed by comparing them with the resolved terms of the corresponding filtered equation. To this end, the DNS dataset is probed using ensemble-averaged values of 10 snapshots at different FTTs. Similarly as in the previous section, the fields are averaged along the streamwise and spanwise direction and displayed as a function of the wall-normal directions. Two different filter widths are assessed as presented in previous Section 4.2, 𝛥∕𝛥=2 and  𝛥∕𝛥=4. In order to deal with unitless numbers, each unclosed term is normalized as follows 𝛼1⋆=𝛼1𝛿 𝑢2 𝑏𝜌𝑏 , 𝛼2⋆=𝛼2𝛿 𝑢2 𝑏𝜌𝑏 , 𝛼3⋆=𝛼3𝛿 𝑢2 𝑏𝜌𝑏 , 𝛼4⋆=𝛼4𝛿 𝑢2 𝑏𝜌𝑏 , 𝛼5⋆=𝛼5𝛿 𝑢2 𝑏𝜌𝑏 , 𝛼6⋆=𝛼6 𝑢𝑏𝜌𝑏 , (49) where subscript ()𝑏corresponds to the bulk value and the resulting normalized quantity is indicated as ()⋆. This normalization is also applied to the remaining fields of each transport equation, viz. convective, pressure and viscous terms. It is important to highlight that this normalization is only used to for the sake of presenting dimensionless results, as the same scaling is applied to each term, i.e. the relative importance of each term remains unchanged. The Journal of Supercritical Fluids 207 (2024) 106191 16 M. Bernades et al. Fig. 11. PDF of the normalized eigenvector 3-principal directions orientation 𝜃∈[−𝜋, 𝜋]azimuth represented on the polar map at different 𝑥-𝑧slices corresponding to (a) 𝑦∕𝛿=0.2 and (b) 𝑦∕𝛿=1.0 and (c) 𝑦∕𝛿=1.8 for adaptive box filter with  𝛥∕𝛥=2 at top row and  𝛥∕𝛥=4 at bottom row. Fig. 12. PDF of filtered DNS and SFS model magnitude of the eigendecomposition normalized by 𝑢𝑏2on 𝑥-𝑧slices corresponding at different y-positions for (a) 𝑦∕𝛿=0.2, (b) 𝑦∕𝛿=1.0 and (c) 𝑦∕𝛿=1.8 for adaptive box filter with  𝛥∕𝛥=2. Table 6 Correlation coefficient between 𝜏𝑎 𝑖𝑗 and 𝜏𝑖𝑗 𝑎𝑆𝐹𝑆 for different models and between 𝜏𝑘𝑘 and 𝜏𝑘𝑘𝑆𝐹 𝑆 for compressible-flow trace-based SFS models at different wall-normal (𝑦∕𝛿) positions at two different filter widths (  𝛥∕𝛥). 𝑦∕𝛿 𝛥∕𝛥 𝐶(𝜏𝑎 𝑖𝑗 , 𝜏𝑖𝑗 𝑎𝑆𝐹𝑆1)𝐶(𝜏𝑎 𝑖𝑗 , 𝜏𝑎 𝑖𝑗 𝑆𝐹𝑆2)𝐶(𝜏𝑎 𝑖𝑗 , 𝜏𝑖𝑗 𝑎𝑆𝐹𝑆3)𝐶(𝜏𝑎 𝑖𝑗 , 𝜏𝑖𝑗 𝑎𝑆𝐹𝑆4)𝐶(𝜏𝑘𝑘, 𝜏𝑘𝑘𝑆𝐹𝑆1)𝐶(𝜏𝑘𝑘, 𝜏𝑘𝑘𝑆𝐹 𝑆2) 0.2 2 0.14 0.17 0.12 0.73 0.85 0.85 4 0.15 0.17 0.15 0.66 0.85 0.84 1.0 2 0.16 0.21 0.11 0.80 0.83 0.83 4 0.18 0.22 0.13 0.77 0.84 0.84 1.8 2 0.14 0.18 0.11 0.77 0.85 0.85 4 0.16 0.16 0.12 0.67 0.82 0.82 interest. Second, the anisotropic SFS stress tensor is compared to the four models presented: (i) Dynamic Smagorisnky, (ii) anisotropic minimum dissipation (AMD), (iii) WALE, and (iv) scale-similarity model. In this case, Fig. 13 depicts the PDF distributions on a barycentric map for all the models; the shapes are quite different compared with those resulting from the filtered DNS highlighted in Fig. 9. Several observations can be obtained from thse plots: (i) Dynamic Smagorinsky provide disk-like distributions near the bottom/cold wall, evolving towards sphere-like in the center and back to most probable disklike states near the top/hot wall as the model is only controlled by The Journal of Supercritical Fluids 207 (2024) 106191 17 M. Bernades et al. Fig. 13. PDF of the anisotropic tensor 𝜏𝑎 𝑖𝑗 from DNS and the SFS models represented on the barycentric map at different 𝑥-𝑧slices corresponding to left 𝑦∕𝛿=0.2, center 𝑦∕𝛿=1.0 and right column 𝑦∕𝛿=1.8 for adaptive box filter with  𝛥∕𝛥=2. First row for 𝜏𝑆𝐹𝑆1 𝑖𝑗 (dynamic Smagorinsky), second row for 𝜏𝑆𝐹𝑆2 𝑖𝑗 (AMD), third row for 𝜏𝑆𝐹𝑆3 𝑖𝑗 (WALE) and fourth row for 𝜏𝑆𝐹𝑆4 𝑖𝑗 (scale-similarity). the strain tensor, (ii) AMD results in non-determined state, near the bottom/cold wall the structures are constrained by diskand spherelike but spread towards even rod-like state elsewhere, (iii) WALE is restricted in rodand disk-like structures near the bottom/cold wall and spread towards sphere-like away from this region and (iv) Scalesimilarity is biased to 2D and 3D shapes in the bottom/cold wall and resulting to non-determined state away from this wall. Third, in terms of orientation, Figs. 14 and 15 display the eigenvector orientations for 𝜙and 𝜃, respectively. It can be seen that all models capture the 90 deg orientation of the first direction as most probable state and both WALE and Scale-Similarity results in non-determined state for second and third components, similar to the filtered DNS. Instead, the azimuthal angle results in similar (poor) performance across models, except for the scale-similarity which predicts relatively well the three orientations The Journal of Supercritical Fluids 207 (2024) 106191 18 M. Bernades et al. Fig. 14. PDF of the normalized eigenvector 3-principal directions orientation 𝜙∈[0, 𝜋]represented on the polar map for the filtered DNS and SFS models at different 𝑥-𝑧slices corresponding to left 𝑦∕𝛿=0.2, center 𝑦∕𝛿=1.0 and right column 𝑦∕𝛿=1.8 for adaptive box filter with  𝛥∕𝛥=2. First row for 𝜏𝑆𝐹𝑆1 𝑖𝑗 (dynamic Smagorinsky), second row for 𝜏𝑆𝐹𝑆2 𝑖𝑗 (AMD), third row for 𝜏𝑆𝐹𝑆3 𝑖𝑗 (WALE) and fourth row for 𝜏𝑆𝐹𝑆4 𝑖𝑗 (scale-similarity). across the wall-normal direction. The eigenvectors of the SFS models are, therefore, aligned with the strain rate tensor, and consequently introduce small-scale dissipation into the system. 6.3. SFS equation of state analysis Fig. 16 reports the results of the two SFS models considered for the equation of state, overlaid onto the filtered and LES pressure field presented in Fig. 6. Therefore, shown are (i) 𝑃(𝜌, 𝑇 ), (ii) the LES term 𝑃(𝜌,  𝑇)and the two models for 𝛼6, resulting in (iii) 𝑃(𝜌,  𝑇) + 𝛿𝑃for the Taylor expansion and (iv) 𝑃(𝜌,  𝑇) + 𝑐𝑝[𝑃(𝜌,  𝑇) − 𝑃(𝜌,   𝑇)] for the ILA model. As it can be observed, neither of the models is able to successfully overlay the filtered LES term and the 𝑃(𝜌, 𝑇 ), and apparently at larger filter widths the Taylor expansion model suffers from some low frequency oscillations, but then ILA completely fails to predict the filtered term. With the purpose of isolating this behavior, a onedimensional test was conducted at constant pressure of 2𝑃𝑐, and with temperature increasing linearly from 0.75−1.5𝑇𝑐(equivalent conditions as the reference DNS dataset) to force the fluid to operate across the pseudo-boiling line. This thermodynamic simulation was done with a uniform mesh at a resolution of 𝛥∕𝛿=0.005 across the 𝛿=1𝑚domain solving only the thermophysical fields. Fig. 17 shows the results along the longitudinal direction at two different filter widths, reporting the same terms as in the previous Figure. This test clearly highlights the models performance in a more controlled environment. As a result, the Taylor expansion model 𝑃(𝜌,  𝑇) + 𝛿𝑃struggles to accurately predict the filtered pressure field and deteriorates the correlation, in fact, its effect becomes worse at larger filter widths. Instead, the ILA model 𝑃(𝜌,  𝑇) + 𝑐𝑝[𝑃(𝜌,  𝑇) − 𝑃(𝜌,   𝑇)] collapses onto the filtered reference, although some low amplitude oscillations are present and not visible within the resolution of these figures. Regardless, this model would be preferred under transcritical conditions. This test inherently highlights the importance of using a PEP numerical framework for transcritical simulations. The pressure field based on the filtered transport quantities struggles to accurately predict the pressure field due to the strong nonlinearity even at constant pressure. Filtering, on the other hand, can amplify oscillations due to the interaction with thermodynamic nonlinearities, particularly across the pseudo-boiling line, as previously reported in [12,28,85]. 6.4. Closure expressions The activity of the terms assessment in Section 5.1 has clearly identified the unclosed terms that require modeling and those that, instead, can be neglected. Existing models are available for the terms 𝛼1(anisotropic part and trace) and 𝛼6. Therefore, for completeness, this The Journal of Supercritical Fluids 207 (2024) 106191 19 M. Bernades et al. Fig. 15. PDF of the normalized eigenvector 3-principal directions orientation 𝜃∈[−𝜋, 𝜋]azimuth represented on the polar map for the filtered DNS and SFS models at different 𝑥-𝑧slices corresponding to left 𝑦∕𝛿=0.2, center 𝑦∕𝛿=1.0 and right column 𝑦∕𝛿=1.8 for adaptive box filter with  𝛥∕𝛥=2. First row for 𝜏𝑆𝐹𝑆1 𝑖𝑗 (dynamic Smagorinsky), second row for 𝜏𝑆𝐹𝑆2 𝑖𝑗 (AMD), third row for 𝜏𝑆𝐹𝑆3 𝑖𝑗 (WALE) and fourth row for 𝜏𝑆𝐹𝑆4 𝑖𝑗 (scale-similarity). section summarizes each of the unclosed terms based on their importance, according to the results from Section 5.1; closure expressions are proposed for the important (i.e., non-negligible) quantities. To this extent, the unresolved fluxes can be modeled using different approaches; here, for simplicity, only three models are considered. First, the scalar fluxes can be described based on the gradient assumption [86], which for general non-linear terms is expressed as  𝑢𝑖𝛩−𝑢𝑖 𝛩= −𝜈𝑆𝐹 𝑆 1 𝐾𝑆𝐹 𝑆 𝜕 𝛩 𝜕𝑥𝑖 ,(55) where 𝛩is a general scalar flux, the turbulent eddy-viscosity can be estimated from the SFS models 𝜈𝑆𝐹 𝑆 , and 𝐾𝑆𝐹 𝑆 is a relevant dimensionless number of the system. This approach is often utilized for species The Journal of Supercritical Fluids 207 (2024) 106191 20 M. Bernades et al. Fig. 16. Ensemble-average of normalized unclosed terms along the wall-normal direction for equation of state for filtered DNS, LES term and both unclosed models for 𝛼6with adaptive box filter with (a)  𝛥∕𝛥=2 and (b)  𝛥∕𝛥=4. Fig. 17. One-dimensional constant pressure test with 𝑁2at 2𝑃𝑐with temperature linearly increasing from 0.75𝑇𝑐to 1.5𝑇𝑐. Results normalized to critical pressure for filtered DNS field 𝑃(𝜌, 𝑇 ), LES filtered-based 𝑃(𝜌,  𝑇)and both unclosed models for 𝛼6, the Taylor expansion 𝑃(𝜌,  𝑇) + 𝛿𝑃and ILA model 𝑃(𝜌,  𝑇) + 𝑐𝑝(𝑃(𝜌,  𝑇) − 𝑃(𝜌,   𝑇)), labeled as 𝑃𝐼𝐿𝐴 for brevity, with adaptive box filter with (a)  𝛥∕𝛥=2 and (b)  𝛥∕𝛥=4. transport, or heat/enthalpy flux closures where the Schmidt [87] and Prandtl numbers [88] of the unresolved scales are the relevant dimensionless numbers respectively. Second, the scale similarity assumption [86] whose closure expressions reads 𝐹1𝐹2−𝐹1𝐹2=1 𝐾𝛼(𝐹1𝐹2−𝐹1𝐹2),(56) where 𝐹1and 𝐹2are the filtered fields and 𝐾𝛼is a constant that depends on the DNS data. Third, Germano et al. [89] proposed an extension of this scale similarity assumption to FG-level, imposing a similarity as 𝐹1𝐹2−𝐹1𝐹2=1 𝐾𝐹𝐺 𝛼( 𝐹1𝐹2− 𝐹1 𝐹2),(57) where the expression is adjusted based on the constant 𝐾𝐹𝐺 𝛼. In the following, each unclosed term is reviewed and potential new closures are proposed. 𝛼1.The SFS stress tensor is the most relevant unclosed term in the momentum equations. The analysis performed against state-of-the-art SFS models confirms that non-eddy-viscosity type models achieve higher correlation factors with respect to eddy-viscosity-like models. In fact, neither of the models represents the same shape (in terms of eigendecomposition) of the filtered DNS, nor it possesses a similar orientation. Nonetheless, eddy-viscosity models provide the orientation of filtered DNS in terms of the first eigen-direction on the polar, and result in poor performance for azimuthal angle. In addition, the trace tensor proposed by Vreman provides higher probability to achieve 𝜏𝑘𝑘 closer to filtered DNS. 𝛼2.The unclosed term arising from viscous terms is typically neglected in incompressible or low-pressure compressible flows. However, the high-pressure results showed an increase of its importance with respect to the SFS stress tensor and its viscous filtered counterpart. Although this quantity is still relatively low, it becomes important in the streamwise and wall-normal directions and particularly in the vicinity of the walls, hence a model may be necessary to accurately predict the effect of the associated unresolved scales. Recalling that 𝛼2= ∇ ⋅(𝝈−𝜎), using the first approach of the ones mentioned above would imply that 𝛩= ∇𝜎and the closure expression would depend on the divergence of the viscous tensor. However, this approach would require a coefficient (10−6)and does not predict the viscous flux well. Instead, the scalesimilarity model turns out to be the most suitable candidate with constant defined as 𝐾𝛼2=2.0. Fig. 18(a) depicts the ensemble-average of the closure expression for stream-wise momentum. It is observed that the model captures the behavior of the filtered field near the walls, however, in the wall-normal direction it is not as accurate. Of note, the Germano approach gives similar performance to scale similarity. A posteriori analysis will confirm the suitability of this model, or in contrast, the effects of neglecting it. 𝛼3.This term can be neglected, weighing (10−6)with respect to the total flux, and a closure expression appears therefore to be unnecessary for this framework. Again, this assumption should be verified a posteriori. 𝛼4.This term achieves a similar order of magnitude as its filtered counterpart. Given the non-linearity of the three terms of the expression, a The Journal of Supercritical Fluids 207 (2024) 106191 21 M. Bernades et al. Fig. 18. Ensemble-average in wall-normal direction of SFS closure expressions for (a) 𝛼2in streamwise direction (b) 𝛼4and (c) 𝛼5for adaptive box filter with  𝛥∕𝛥=2. suitable approach appears to be the scale similarity one, which in this case would read 𝛼4𝑆𝐹 𝑆 =𝜌 𝑐2∇⋅ 𝐮−𝜌𝑐2∇⋅𝐮=1 𝐾𝛼4(𝜌 𝑐2∇⋅ 𝐮−𝜌 𝑐2∇⋅ 𝐮),(58) where in this case 𝐾𝛼4= −10 minimizes the squared errors. The robustness of this factor needs to be tested in a posteriori analyses. In Fig. 18(b) the improvement of the proposed closure expression can be identified resulting in a lower difference with respect to the unclosed trend, although this improvement is insignificant and the constant adjustment is not enough to attenuate the rapid fluctuations of this term. In this case, the similarity of the velocity fields based on Favrefiltering is causing some oscillations, instead non-density weighted filtering is applied as described in the model equation 𝛼4𝑆𝐹 𝑆 . 𝛼5.The unclosed contribution arising from viscous and Fourier terms, also related to non-linear thermodynamic variables, becomes important across the domain and a closure expression is therefore recommended. To this aim, the following expression is proposed analogously as for 𝛼4 𝛼5𝑆𝐹 𝑆 =[1 𝜌 𝛽𝑣 𝑐𝑣𝛽𝑇 (𝝈∶ ∇ ⊗𝐮− ∇ ⋅𝒒) − 1 𝜌  𝛽𝑣 𝑐𝑣 𝛽𝑇 ( 𝝈∶ ∇ ⊗ 𝐮− ∇ ⋅ 𝒒)] =1 𝐾𝛼5⎡⎢⎢⎣ 1 𝜌 𝛽𝑣 𝑐𝑣𝛽𝑇 (𝝈∶ ∇ ⊗𝐮− ∇ ⋅𝒒) − 1 𝜌  𝛽𝑣 𝑐𝑣 𝛽𝑇 ( 𝝈∶ ∇ ⊗ 𝐮− ∇ ⋅ 𝒒)⎤⎥⎥⎦ , (59) with 𝐾𝛼5=0.4. In this case, the a priori results based on DNS provides very good agreement with the filtered field, as depicted in Fig. 18(c). 𝛼6.Although the error coming from not closing this term is insignificant with respect to the filtered pressure, the ILA model is the preferred candidate. In the one-dimensional high-pressure thermophysical test, the model adjusted the unclosed subfilter term to accurately predict the filtered field. 7. Conclusions A novel LES framework has been presented based on a recently proposed discretization approach for transcritical fluids turbulence. This method attains kinetic-energy and pressure-equilibrium preservation, and is able to accomplish stable and non-dissipative scale-resolving simulations; it is therefore suitable for high-fidelity LES. Based on this framework, the set of filtered equations has been derived and the resulting subfilter terms have been identified. To perform an 𝑎 𝑝𝑟𝑖𝑜𝑟𝑖 analysis, a transcritical channel flow DNS was computed. This dataset was then filtered with two different filtered widths, corresponding to a wellresolved LES and a coarse LES. Also, the activity of the resolved and subfilter terms was analyzed to understand their relative importance, and the subfilter stress tensor was assessed in terms of magnitude, shape and orientation through eigendecomposition analysis. The set of filtered equations generates six unclosed terms. The resolved and subfilter term-by-term analysis showed that the subfilter stress tensor dominates over the viscous subfilter-term in the momentum equation. From the filtered pressure equation, three unclosed terms emerge. The terms associated with flow dilatation and diffusion are significant, whereas the SFS connected to the power of pressure forces can be neglected. Moreover, an additional unclosed term results from the non-linear equation of state. Based on an eigendecomposition analysis of the subfilter stress tensor, it can be concluded that its largest magnitude is achieved in the vicinity of the hot wall, whereas its minimum lies at the center of the channel. Here, both trace models analyzed resulted in good prediction of the mean of the trace and the variance with slight more favorable match with Vreman model; this result will be verified in a posteriori analysis as future work. In terms of shape, the largest probability showed a rod-like state near the cold wall, slightly transitioning to sphere-like at the center, although still dominated by rod-like type structures. The orientation of the polar angle is dominated by the first component, yielding 90◦with largest likelihood, the second principal direction tends to 30◦near the walls and non-determined state at the center, whereas the third orientation is slightly more favorable tilted towards 30◦at the center. Instead, the azimuthal angle stays at 0◦ for first direction near the cold wall and marginally leaned to 20◦, and 90◦and 270◦for second and third directions. The subsequent a priori analysis of the eigendecomposition on the subfilter models reported prediction discrepancies on magnitude, shape and orientation. These results could lay the foundation for consequent modeling efforts. From a correlation factor standpoint, models of the non-eddy-viscosity type showed superior performances compared to eddy-viscosity, with WALE showing slightly satisfactory performances among the eddy-viscosity models examined. On the other hand, the assessment of SFS models for the equation of state suggests that the ILA model can accurately predict the pressure field, while the Taylor expansion is more susceptible to non-liner oscillations. Both models worsen their performance at larger filter widths. Overall, the results indicate that non-eddy-viscosity-based models are suitable for the SFS stress tensor, along with the Vreman trace tensor model, while ILA model is a good candidate for the equation of state. Closure expressions based on the similarity assumptions were proposed for the unconventional terms that showed significant magnitude, i.e., 𝛼4and 𝛼5, instead 𝛼3appears to be negligible. Future work should focus on following up on this a priori analysis by implementing the selected models and computing non-dissipative LES of wall-bounded transcritical turbulent flows, expanding the results for larger filter widths and Reynolds numbers. To this aim, a posteriori analysis will be required to evaluate the suitability and accuracy of the framework by comparing the results against the DNS dataset. In particular, the dynamic assessment of firstand second-order statistics and also its performance in terms of capturing the physics-based heat transfer phenomena at the wall should be carefully examined. In addition, the characteristic local minimum behavior of the convective The Journal of Supercritical Fluids 207 (2024) 106191 22 M. Bernades et al. SFS term near the cold wall and their fluctuations in the pseudoboling region needs to be also investigated in detail. Independently, the LES implementation will require wall-model and/or wall-resolved approach which will need to be carefully developed and validated. This purported a posteriori analysis will eventually quantify the specific choices for SFS models and provide feedback for the presented closure expressions Given the current limitation of existing SFS models, more efforts are needed towards the development and implementation of more advanced physics-based expressions guided, for example, by eigendecompostion approaches, compared to the simple closure expressions proposed. In particular, the findings of this work are focused on relatively low Reynolds numbers at low-Mach-number conditions. Therefore, to achieve more universal findings, a broader range of operating conditions should be investigated. In addition, buoyancy effects need to be also considered in future works as gravity forces may become important when considering a wider range of operating conditions. CRediT authorship contribution statement Marc Bernades: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. Lluís Jofre: Writing – review & editing, Supervision, Methodology, Funding acquisition, Conceptualization. Francesco Capuano: Writing – review & editing, Supervision, Methodology, Conceptualization. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability & reproducibility The data reported in this paper was produced by an in-house software available as open-source. A MATLAB code, labeled as LargeEddy Simulation for Transcritical Turbulence (LES-TT), was used to obtain part of the numerical results presented in this work. It can be accessed at https://github.com/marc-bernades/LES-TT. LES-TT was designed to serve as a flexible tool to filter and assess DNS data and build and evaluate SFS models. The code is equipped with several comments for readability. Additionally, the repository includes (i) a description of the code, (ii) instructions and a guide for users, (iii) DNS dataset, and (iv) a filtered file based on ensemble averaging of multiple DNS snapshots. The filtered DNS dataset studied in Section 4.1, the term-by-term analysis presented in Section 5.1, the eigendecomposition analysis of Section 5.2 and the implementation and a priori analysis of the SFS models assessed in Section 6can be fully reproduced by the interested reader. The open-source Reproducible Hybrid-architecture flow solver Engineered for Academia (RHEA) [72] used for computing the DNS dataset can be accessed at: https://gitlab.com/ProjectRHEA/ flowsolverrhea. RHEA is written in C++, using object-oriented programming, utilizes YAML and HDF5 for input/output operations, and targets hybrid supercomputing architectures. RHEA was utilized to obtain the transcritical 3D channel flow data, which have been filtered to perform the a priori analysis. The test is also available in the repository and can be reproduced by the user. Data will be stored locally on clusters at Universitat Polit` cnica de Catalunya ⋅BarcelonaTech (UPC), and will be provided upon making proper arrangements with the requesters. Acknowledgments The authors gratefully acknowledge the Formació de Professorat Universitari scholarship (FPU-UPC R.D 103/2019) of the Universitat Politècnica de Catalunya ⋅BarcelonaTech (UPC) (Spain), the Serra Húnter and SGR (2021-SGR-01045) programs of the Generalitat de Catalunya (Spain), the Beatriz Galindo program (Distinguished Researcher, BGP18/00026) of the Ministerio de Educación yFormación Profesional (Spain), and the computer resources at FinisTerrae III & MareNostrum and the technical support provided by CESGA & Barcelona Supercomputing Center (RES-IM-2023-1-0005, RES-IM2023-2-0005). Funding sources This work is funded by the European Union (ERC, SCRAMBLE, 101040379). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. References [1] S.B. Pope, Turbulent Flows, first ed., Cambridge University Press, Cambridge (UK), 2000. [2] E. Garnier, N. Adams, P. Sagaut, Large Eddy Simulation for Compressible Flows, Springer Science & Business Media, 2009. [3] H. Pitsch, Large-eddy simulation of turbulent combustion, Annu. Rev. Fluid Mech. 38 (2006) 453–482. [4] N.J. Georgiadis, D.P. Rizzeta, C. Fureby, Large-eddy simulation: current capabilities, recommended practices, and future research, AIAA J. 48 (2010) 1772–1784. [5] R.O. Fox, Large-eddy-simulation tools for multiphase flows, Annu. Rev. Fluid Mech. 44 (2012) 47–76. [6] L. Jofre, J. Urzay, A characteristic length scale for density gradients in supercritical monocomponent flows near pseudoboiling, in: Annual Research Briefs, Center for Turbulence Research, Stanford University, 2020, pp. 277–282. [7] L. Jofre, J. Urzay, Transcritical diffuse-interface hydrodynamics of propellants in high-pressure combustors of chemical propulsion systems, Prog. Energy Combust. Sci. 82 (2021) 100877. [8] M. Bernades, F. Capuano, L. Jofre, Microconfined high-pressure transcritical fluid turbulence, Phys. Fluids 35 (2023) 015163. [9] M. Bernades, L. Jofre, Thermophysical analysis of microconfined turbulent flow regimes at supercritical fluid conditions in heat transfer applications, J. Heat Transfer 144 (2022) 082501. [10] L. Jofre, M. Bernades, F. Capuano, Dimensionality reduction of non-buoyant microconfined high-pressure transcritical fluid turbulence, Int. J. Heat Fluid Flow 102 (2023) 109169. [11] N. Masclans, F. Vázquez-Novoa, M. Bernades, R.M. Badia, L. Jofre, Thermodynamics-informed neural network for recovering supercritical fluid thermophysical information from turbulent velocity data, Int. J. Thermofluids 20 (2023) 100448. [12] M. Bernades, L. Jofre, F. Capuano, Kinetic-energyand pressure-equilibriumpreserving schemes for real-gas turbulence in the transcritical regime, J. Comput. Phys. 493 (2023) 112477. [13] J. Bellan, Theory, modeling and analysis of turbulent supercritical mixing, Combust. Sci. Technol. 178 (2006) 253–281. [14] L. Selle, A.N. Okongo’o, J. Bellan, K.G. Harstad, Modelling of subgrid-scale phenomena in supercritical transitional mixing layers: an a priori study, J. Fluid Mech. 593 (2007) 54–91. [15] E.S. Taşkinoğlu, J. Bellan, Subgrid-scale models and large-eddy simulation of oxygen stream disintegration and mixing with a hydrogen or helium stream at supercritical pressure, J. Fluid Mech. 679 (2011) 156–193. [16] T. Schmitt, L. Selle, B. Cuenot, T. Poinsot, Large-Eddy Simulation of Transcritical Flows, Vol. 337, 2009, pp. 528–538. [17] Y. Ren, Z. Wu, X. Meng, G. Ou, J. Kou, H. Jin, L. Guo, Large eddy simulation of water jets under transcritical and supercritical conditions, J. Supercrit. Fluids 187 (2022) 105648. [18] H. Wang, J. Zang, J. Wang, Y. Huang, Large eddy simulation of the heat transfer and unsteady pulsation of supercritical carbon dioxide in a square subchannel, Int. J. Therm. Sci. 172 (2022) 107377. [19] J. Smagorinsky, General circulation experiments with the primitive equations, Mon. Weather Rev. 91 (1963) 99–164. [20] F. Nicoud, F. Ducros, Subgrid-scale stress modelling based on the square of the velocity gradient tensor, Flow Turbul. Combust. 62 (1999) 183–200. The Journal of Supercritical Fluids 207 (2024) 106191 23 M. Bernades et al. [21] F. Ducros, V. Ferrand, F. Nicoud, C. Weber, D. Darracq, C. Gacherieu, T. Poinsot, Large-eddy simulation of the shock/turbulence interaction, J. Comput. Phys. 1252 (1999) 517–549. [22] X. Petit, G. Ribert, G. Lartigue, P. Domingo, Large-eddy simulation of supercritical fluid injection, J. Supercrit. Fluids 84 (2013) 104843. [23] B. Vreman, An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications, Phys. Fluids 16 (2004) 3670. [24] G. Borghesi, J. Bellan, A priori and a posteriori investigations for developing large eddy simulations of multi-species turbulent mixing under high-pressure conditions, Phys. Fluids 27 (2015) 035117. [25] H. Müller, C.A. Niedermeier, J. Matheis, M. Pfitzner, S. Hickel, Large-eddy simulation of nitrogen injection at transand supercritical conditions, Phys. Fluids 28 (2016) 015102. [26] S. Hickel, N.A. Adams, J.A. Domaradzki, An adaptive local deconvolution method for implicit les, J. Comput. Phys. 213 (2006) 413–436. [27] U. Unnikrishnan, J.C. Oefelein, V. Yang, Subgrid modeling of the filtered equation of state with application to real-fluid turbulent mixing at supercritical pressures, Phys. Fluids 34 (2022) 065112. [28] G. Lacaze, T. Schmitt, A. Ruiz, J. Oefelein, Comparison of energy-,pressureand enthalpy-based approaches for modeling supercritical flows, Comput. & Fluids 181 (2019) 35–56. [29] H. Terashima, M. Koshi, Approach for simulating gas-liquid-like flows under supercritical pressures using a high-order central differencing scheme, J. Comput. Phys. 231 (2012) 6907–6923. [30] M. Bernades, L. Jofre, F. Capuano, Investigation of a Novel Numerical Scheme for High-Pressure Supercritical Fluids Turbulence, Proceedings of the Summer Program 2022, Center for Turbulence Research, Stanford University, 2022, pp. 225–234. [31] N. Shima, Y. Kuya, Y. Tamaki, S. Kawai, Preventing spurious pressure oscillations in split convective form discretization for compressible flows, J. Comput. Phys. 427 (2021) 110060. [32] R. Mittal, P. Moin, Suitability of upwind-biased finite difference schemes for large-eddy simulation of turbulent flows, AIAA J. 35 (1997) 1415–1417. [33] G. Coppola, F. Capuano, L. de Luca, Discrete energy-conservation properties in the numerical simulation of the Navier-Stokes equations, Appl. Mech. Rev. 71 (2019) 010803. [34] D.Y. Peng, D.B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15 (1976) 59–64. [35] W.C. Reynolds, P. Colonna, Thermodynamics: Fundamentals and Engineering Applications, first ed., Cambridge University Press, Cambridge (UK), 2019. [36] A. Burcat, B. Ruscic, Third Millennium Ideal Gas and Condensed Phase Thermochemical Database for Combustion with Updates from Active Thermochemical Tables, Technical Report, Argonne National Laboratory, 2005. [37] T.H. Chung, L.L. Lee, K.E. Starling, Applications of kinetic gas theories and multiparameter correlation for prediction of dilute gas viscosity and thermal conductivity, Ind. Eng. Chem. Fund. 23 (1984) 8–13. [38] T.H. Chung, M. Ajlan, L.L. Lee, K.E. Starling, Generalized multiparameter correlation for nonpolar and polar fluid transport properties, Ind. Eng. Chem. Fund. 27 (1988) 671–679. [39] B.E. Poling, J.M. Prausnitz, J.P. O’Connell, Properties of Gases and Liquids, fifth ed., McGraw Hill, New York (USA), 2001. [40] P.J. Linstrom, W.G. Mallard, Thermophysical properties of fluid systems, in: NIST ChEmistry Webbook (SRD 69), 2021. [41] G. Coppola, F. Capuano, S. Pirozzoli, L. de Luca, Numerically stable formulations of convective terms for turbulent compressible flows, J. Comput. Phys. 382 (2019) 86–104. [42] M. Bernades, F. Capuano, F.X. Trias, L. Jofre, Energy-preserving stable computations of high-pressure supercritical fluids turbulence, in: 8th European Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS, 2022, pp. 1–12. [43] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev. 43 (2001) 89–112. [44] P.E. Lapenna, R. Lamioni, P.P. Ciottoli, F. Creta, Low-Mach number simulations of transcritical flows, AIAA J. 128 (2018) 0346. [45] A. Leonard, Energy cascade in large-eddy simulations of turbulent fluid flows, Adv. Geophys. 18 (1974) 237–248. [46] A. Favre, Turbulence: space–time statistical properties and behavior in supersonic flows, Phys. Fluids 26 (1983) 2851–2863. [47] B. Vreman, Direct and Large-Eddy Simulation of the Compressible Turbulent Mixing Layer, Universiteit Twente Enschede, 1995. [48] B. Geurts, B. Vreman, H. Kuerten, Comparison of dns and les of transitional and turbulent compressible flow: flat plate and mixing layer, in: Application of Direct and Large Eddy Simulation to Transition and Turbulence, AGARD Conference Proceedings, AGARD, France, 1994, pp. 5/1–5/14, 74th Fluid Dynamics Symposium 1994 ; Conference date: 18-04-1994 Through 21-04-1994. [49] S. Ghosal, P. Moin, The basic equations for the large eddy simulation of turbulent flows in complex geometry, J. Comput. Phys. 118 (1995) 24–37. [50] O.V. Vasilyev, T.S. Lund, P. Moin, A general class of commutative filters for LES in complex geometries, J. Comput. Phys. 146 (1998) 82–104. [51] A.L. Marsden, O.V. Vasilyev, P. Moin, Construction of commutative filters for LES on unstructured meshes, J. Comput. Phys. 175 (2002) 584–603. [52] T.A. Zang, R.B. Dahlburg, J.P. Dahlburg, Direct and large-eddy simulations of three-dimensional compressible navier–stokes turbulence, Phys. Fluids 4 (1992) 127. [53] A. Yoshizawa, Statistical theory for compressible turbulent shear flows, with the application to subgrid modeling, Phys. Fluids 29 (1986) 2152–2164. [54] B. Vreman, B. Geurts, H. Kuerten, On the formulation of the dynamic mixed subgrid-scale model, Phys. Fluids 6 (1994) 4057–4059. [55] R.S. Rogallo, P. Moin, Numerical simulation of turbulent flow, Annu. Rev. Fluid Mech. 16 (1984) 2150. [56] M. Germano, Turbulence: the filtering approach, J. Fluid Mech. 238 (1992) 325–336. [57] Y. Zang, R.L. Street, J.P. Koseff, A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows, Phys. Fluids 5 (1993) 3186–3196. [58] B. Vreman, B. Geurts, H. Kuerten, Large eddy simulation of the temporal mixing layer using the clark model, Theor. Comput. Fluid Dyn. 8 (1996) 309–324. [59] R.W.C.P. Verstappen, S. Bose, J. Lee, H. Choi, P. Moin, A dynamic eddy-viscosity model based on the invariants of the rate-of-strain, in: Proceedings of the Summer Program 2010, Center for Turbulence Research, Stanford University, 2010, pp. 183–192. [60] R.W.C.P. Verstappen, W. Rozema, H.J. Baea, Numerical Scale Separation in Large-Eddy Simulation, Proceedings of the Summer Program 2014, Center for Turbulence Research, Stanford University, 2014, pp. 417–426. [61] W. Rozema, H.J. Bae, P. Moin, R. Verstappen, Minimum-dissipation models for large-eddy simulation, Phys. Fluids 27 (2015) 085107. [62] M. Abkar, P. Moin, Les of the convective boundary layer: A minimum-dissipation modeling approach, in: Annual Research Briefs, Center for Turbulence Research, Stanford University, 2016, pp. 87–96. [63] M.H. Silvis, R.A. Remmerswaal, R. Verstappen, Physical consistency of subgridscale models for large-eddy simulation of incompressible turbulent flows, Phys. Fluids 29 (2017) 015105. [64] H.B. Toda, O. Cabrit, G. Balarac, s. Bose, J. Lee, H. Choi, F. Nicoud, A subgrid-scale model based on singular values for les in complex geometriess, in: Proceedings of the Summer Program 2010, Center for Turbulence Research, Stanford University, 2010, pp. 193–202. [65] M. Rieth, F. Proch, O. Stein, M. Pettit, A. Kempf, Comparison of the sigma and smagorinsky les models for grid generated turbulence and a channel flow, Comput. & Fluids 99 (2014) 172–181. [66] J. Bardina, J.H. Ferziger, W.C. Reynolds, Improved Turbulence Models Based on LES of Homogeneous Incompressible Turbulent Flows, Stanford University, 1983. [67] S. Liu, C. Meneveau, J. Katz, On the properties of similarity subgrid-scale models as deduced from measurements in a turbulent jet, J. Fluid Mech. 275 (1994) 83–119. [68] L. Jofre, S.P. Domino, G. Iaccarino, A framework for characterizing structural uncertainty in large-eddy simulation closures, Flow Turbul. Combust. 100 (2018) 341–363. [69] P. Moin, K. Squires, W. CAbot, S. Lee, A dynamic subgridscale model for compressible turbulence and scalar transport, Phys. Fluids 3 (1991) 2746. [70] B. Vreman, B. Geurts, H. Kuerten, Realizability conditions for the turbulent stress tensor in large-eddy simulation, J. Fluid Mech. 278 (1994) 351–362. [71] L. Selle, G. Ribert, Modeling requirements for large-eddy simulation of turbulent flows under supercritical thermodynamic conditions, in: Proceedings of the Summer Program 2008, Center for Turbulence Research, Stanford University, 2008, pp. 195–207. [72] L. Jofre, A. Abdellatif, G. Oyarzun, RHEA - an open-source Reproducible hybridarchitecture flow solver engineered for academia, J. Open Source Softw. 8 (2023) 4637. [73] A. Abdellatif, J. Ventosa-Molina, J. Grau, R. Torres, L. Jofre, Artificial compressibility method for high-pressure transcritical fluids at low Mach numbers, Comput. & Fluids 270 (2023) 106163. [74] P. Sagaut, R. Grohens, Discrete filters for large eddy simulation, Internat. J. Numer. Methods Fluids 31 (1999) 1195–1220. [75] A.B. Vidal, O. Lehmkuhl, F.X. Trias, C.D. Pérez-Segarra, On the properties of discrete spatial filters for CFD, J. Comput. Phys. 326 (2016) 474–498. [76] S.T. Bose, G.I. Park, Wall-modeled large-eddy simulation for complex turbulent flows, Annu. Rev. Fluid Mech. 50 (2018) 535–561. [77] J. Larsson, S. Kawai, J. Boadart, I. Bermejo-Moreno, Large eddy simulation with modeled wall-stress: recent progress and future directions, Mech. Eng. Rev. 3 (2015) 15–00418–15–00418. [78] L. Jofre, S.P. Domino, G. Iaccarino, Eigensensitivity analysis of subgrid-scale stresses in large-eddy simulation of a turbulent axisymmetric jet, Int. J. Heat Fluid Flow 77 (2019) 314–335. [79] B. Tao, J. Katz, C. Meneveau, Statistical geometry of subgrid-scale stresses determined from holographic velocimetry measurements, J. Fluid Mech. 457 (2002) 35–78. [80] K. Horiuti, Roles of non-aligned eigenvectors of strain-rate and subgrid-scale stress tensors in turbulence generation, J. Fluid Mech. 491 (2003) 65–100. [81] F. Dabbagh, F.X. Trias, A. Gorobets, A. Oliva, On the evolution of flow topology in turbulent Rayleigh-Bénard convection, Phys. Fluids 28 (2016) 115105. [82] F. Dabbagh, F.X. Trias, A. Gorobets, A. Oliva, A priori study of subgrid-scale features in turbulent Rayleigh-Bénard convection, Phys. Fluids 29 (2017) 115109. The Journal of Supercritical Fluids 207 (2024) 106191 24 M. Bernades et al. [83] S. Banerjee, R. Krahl, F. Durst, C. Zenger, Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches, J. Turbul. 8 (2007) N32. [84] M. Bernades, F. Capuano, L. Jofre, Flow Physics Characterization of Microconfined High-Pressure Transcritical Turbulence, Proceedings of the Summer Program 2022, Center for Turbulence Research, Stanford University, 2022, pp. 215–224. [85] M. Bernades, L. Jofre, F. Capuano, Non-dissipative large-eddy simulation of highpressure transcritical turbulent flows: formulation and a priori analysis, in: 14th International ERCOFTAC Symposium on Engineering Turbulence Modelling and Measurements, ERCOFTAC,2023, pp. 546–551. [86] T. Poinsot, D. Veynante, Theoretical and Numerical Combustion, second ed., RT Edwards, Inc., 2005. [87] S. Guo, K. Chen, E. Tsotsas, F. Shang, Z. Ge, H. Jin, Y. Chen, L. Guo, A largeeddy simulation study of the transcritical mixing process in coaxial jet flow under supercritical condition, J. Supercrit. Fluids 203 (2023) 106080. [88] J. Xie, G. Xie, Assessment on heat transfer deterioration to supercritical carbon dioxide in upward flows via large eddy simulation, Int. J. Heat Fluid Flow 95 (2022) 108954. [89] M. Germano, A. Maffio, S. Sello, G. Mariotti, in: J.-P. Chollet (Ed.), On the Extension of the Dynamic Modelling Procedure to Turbulent Reacting Flows, Dordrecht, 1997, pp. 291–300.