The scalar filtered density function (FDF) methodology is extended and employed for large eddy simulation (LES) of turbulent flows under supercritical condition. To describe real fluid behavior, the extended methodology incorporates the generalized heat and mass diffusion models along with real fluid thermodynamic relations which are derived using the cubic Peng–Robinson equation of state. These models are implemented within the stochastic differential equations comprising the scalar FDF transport. Simulations are conducted of a temporally developing mixing layer under supercritical condition and the results are assessed by comparing with data generated by direct numerical simulation (DNS) of the same layer. The consistency of the proposed FDF methodology is assessed. The LES-FDF predictions are shown to agree favorably with the DNS data and exhibit several key features pertaining to supercritical turbulent flows.
High-pressure (high-p) turbulent mixing and reacting flows are encountered in many systems such as diesel engines, rocket propulsion, gas turbines, as well as heat exchangers. These flows are also relevant to supercritical CO2 power cycles which have gained much attention recently because of their higher efficiency, lower capital costs, and more compact components [1–3]. In some of these systems, pressure and temperature are large enough that the working fluid is in supercritical state wherein no phase change occurs but the fluid exhibits large variations in thermodynamic properties such as density and specific heat. Despite important applications of such flows, there is still limited understanding of their turbulent behavior integrated with high-p flow physics. Modeling and simulation of turbulent flows of supercritical nature has been the subject of relatively few investigations, e.g., Refs. [4–8]. Previous studies reveal the intricate physics associated with these flows, modeling of which requires consideration of generalized diffusion (including multicomponent, cross and differential diffusion), real fluid thermodynamic and transport properties, and chemical kinetics together with models describing their interactions with turbulence. It has been shown that these are the key elements to take into account to ensure reliable simulation of high-p flows [4–7,9–13].
Large eddy simulation (LES) has proven quite effective in simulation of turbulent flows. Within the past several decades, there has been a significant progress in application of LES in studying turbulent (reacting) flows, predominantly at atmospheric and moderately elevated pressure conditions. However, the literature pertaining to the LES study of supercritical turbulent flows has been limited; see, e.g., Refs. [3–5,7,14,15]. A challenging issue in LES is accurate representation of the subgrid scale (SGS) quantities. This becomes even more difficult in high-p flows as the nonlinearity associated with high-p effects gives rise to additional unclosed SGS terms; this is even the case for the filtered real fluid equations of state alone . Besides, modeling the turbulence-chemistry interactions poses major challenges to LES of turbulent reacting flows. For this purpose, the filtered density function (FDF) methodology has evolved as one of the promising strategies [17–23] due to its capacity to provide a closed form for the chemical reaction effects. The FDF essentially presents the probability density function (PDF) of the SGS variables. A rigorous means of obtaining the FDF is via its transport equation [24–28]. The simplest form of the transported FDF formulations is the scalar FDF [24,25] which considers the FDF of scalar quantities. This form of the FDF, which is the subject of the present study, is extended in Ref.  to account for differential diffusion. This is accomplished by including the molecular transport of heat and mass as separate drift terms in the scalar stochastic differential equations (SDEs) within the FDF formulation. Another related development in FDF modeling is the extension of the scalar FDF to compressible flows by including the pressure and viscous dissipation effects, as appear in the energy equation, in the modeled FDF transport equation . All previous FDF formulations however involve ideal gas thermodynamic relations together with simplified diffusion described by Fick’s law. To employ the FDF for LES of high-p turbulent reacting flows, in particular those occurring under supercritical condition, these limitations need to be alleviated. The present work focuses on several modeling aspects towards this objective. This includes incorporating real fluid thermodynamic relations and general description of heat and mass diffusion for binary mixtures within the FDF framework. The accuracy of the proposed methodology is assessed by performing LES of a temporally developing mixing layer involving transport of passive scalars at supercritical pressure. The consistency of the FDF is demonstrated and its predictions are compared with those of direct numerical simulation (DNS) of the same flow.
2 Governing Equations
In this section, we present the equations governing high-p compressible flows with inclusion of real fluid thermodynamic relations and general heat and mass diffusion models for binary mixtures. Following the presentation of these models, we discuss their incorporation in the FDF formulation.
2.1 Conservation Equations.
2.2 Equation of State and Thermodynamic Properties.
2.3 Heat and Mass Transport.
2.4 Large Eddy Simulation Equations.
2.5 The Filtered Density Function Transport Equation.
2.6 Monte Carlo Solution of the Filtered Density Function.
3 Numerical Solution Procedure
Numerical implementation of the FDF consists of a hybrid finite difference (FD)/MC method as considered in previous work [24–28]. The MC part involves solution of the SDEs, Eqs. (38) and (39), using a Lagrangian MC method. In this method, the FDF is represented by an ensemble of Np number of particles, which carry information pertaining to the FDF variables: position and scalars , where n = 1,…, Np. These particles are evolved in time by integrating the set of SDEs discretized using the Euler–Maruyamma discretization . The MC particles are distributed randomly and experience free movement within the flow field in a Lagrangian fashion. The statistical information, e.g., filtered quantities, at any point is obtained by considering an ensemble of particles residing within an ensemble domain centered around the point. To obtain reliable statistics with minimal numerical diffusion, it is desired to minimize the size of the ensemble domain and maximize the number of MC particles causing convergence of statistics to the corresponding filtered quantities [24,27]. The FD part involves FD solution of Eqs. (29)–(31) on equally spaced Cartesian grid points using a compact parameter numerical scheme [41,42]. The transfer of information from FD grid points to MC particles is done by interpolation. To reduce the computational cost, non-uniform weights for the particles are considered [25,27,28]. The sum of weights within the ensemble domain corresponds to filtered fluid density at any instant, and the Favre filtered values are constructed from the weighted average values [25,27,28]. To ensure reliable statistics, distribution of MC particles throughout the flow field is monitored instantaneously to ensure their uniform distribution at all times. Message passing interface library is used to parallelize both FD and MC solvers in all three coordinate directions.
4 Flow Configuration and Specifications
Simulations are conducted of a three-dimensional (3D) temporal mixing layer, involving transport of passive scalars. The flow configuration consists of two parallel streams on top and bottom flowing in opposite directions, carrying carbon dioxide and methane, respectively. The unfiltered streamwise velocity, temperature, and mixture fraction fields are initialized using hyperbolic tangent profiles with freestream conditions: u = 104.5 m s−1 and T = 1000 K on the top and u = −42.8 m s−1 and T = 300 K on the bottom. Pressure values of up to 150 atm are considered. The streams are in the supercritical state when their temperature and pressure are above the critical values (304.1 K and 72.8 atm for carbon dioxide and 190.6 K and 45.4 atm for methane). The initial pressure is uniform and the initial density and internal energy fields are determined using real fluid thermodynamic relations. The initial filtered quantities are obtained from the unfiltered ones by applying the LES filtering operation as explained below. Variables are nondimensionalized using reference values. The reference length Lo is set to the initial vorticity thickness of 0.00686 m . The reference velocity is defined as the velocity difference across the layer, Uo = 147.3 m s−1. Time is nondimensionalized using Lo/Uo. The reference temperature (650 K) and density (98.3 kg m−3) are defined as the average of top and bottom stream values. The Reynolds number based on reference values is Re = 200.
The computational domain is a cubic box, 0 ≤ x ≤ L, −(L/2) ≤ y ≤ L/2, 0 ≤ z ≤ L where the coordinates x, y, and z denote the nondimensional streamwise, the cross-stream, and the spanwise directions, respectively. The nondimensional length L is specified as L = Nrλu/Lo, in which Nr is the desired number of successive vortex pairings and λu is the wavelength of the most unstable mode corresponding to the mean streamwise velocity profile imposed at the initial time [8,43]. LES is performed on 33 × 33 × 33 grid points in x, y, and z directions, respectively, and the results are compared with those of DNS of the same layer with the grid resolution of 193 × 193 × 193. The LES filter size is twice as large as the grid spacing in each direction. The ensemble domain size for calculation of filtered quantities from the MC solver is equal to the filter size. The ensemble averaging is performed instantaneously using approximately 64 MC particles at each grid point. As per results of previous work [24–27], this number is sufficient to provide reliable statistics here. The total number of MC particles within the domain is approximately 0.6 million at all times. The treatment of MC particles at the boundaries include periodic movements in x and z directions as well as mirror-reflections in y direction to represent the zero-gradient boundary condition for scalar variables. For comparison, the DNS data are filtered via a top-hat filter. The FD boundary conditions include periodic condition in the streamwise and spanwise directions as well as the Navier–Stokes characteristic boundary condition  on top and bottom boundaries in the cross-stream direction. The formation of large-scale flow structures are expedited using eigenfunction-based initial 3D perturbations [43,45,46] with a random phase shift between the 3D modes. The Reynolds-averaged statistics are constructed from the instantaneous data by spatial averaging over the homogeneous directions (x and z). In the presentation below, the Reynolds-averaged filtered and Favre filtered variables are denoted as and , respectively.
The temporal mixing layer case considered exhibits two successive vortex pairings and 3D flow structures [47,48] as illustrated in Fig. 1. This figure displays the iso-surfaces of CO2 mass fraction at 150 atm predicted by DNS at the nondimensional time of t = 50. By this time, the flow has gone through vortex pairings and as a result, large-scale spanwise rollers have formed along with secondary flow structures in streamwise planes, as shown in Fig. 1.
Fluids under supercritical state exhibit large variations in thermodynamic properties with small changes in temperature and pressure. This is shown in Fig. 2 for CO2. In the absence of phase change, large variation in density occurs when temperature and pressure exceed their critical values. Heat capacity also shows a large peak near the critical point. In this figure, predictions obtained using the PR equation of state show good agreement with data from the National Institute of Standards and Technology (NIST)  at 60 and 100 atm. For comparison, ideal gas equation of state results are also shown. As expected, significant deviation from ideal gas behavior exist, especially near the critical point. It is notable that this is even the case for 60 atm which is still at subcritical state, suggesting that the ideal gas assumption may be unsuitable for flow simulations at such elevated pressures. This is also evident in LES results in Fig. 3 which shows the filtered compression factor for different pressures obtained using the PR equation of state. Compression factor values different from one indicate departure from ideal gas behavior. While at 30 atm, its deviation from one is less than , at 60 atm the difference is too large for the ideal gas law to be sufficient. As shown, by increasing the pressure, the range of compression factor values further increases which necessitates consideration of real fluid thermodynamic models.
Turbulent mixing under supercritical condition includes important features as observed in the present DNS results. The temperature field predicted by DNS in Fig. 4 depicts rolling up of two neighboring vortices on the x–y plane along with secondary structures on the y–z plane. As mixing of the two fluids occurs, the locations where temperature is near its critical value are at the bottom of the layer and the braid regions between the vortices. This causes significant density variations because of the way density varies with temperature near the critical point (as viewed in Fig. 2 for CO2). Supercritical mixing of liquid-like CH4 and CO2 creates finger-like structures  which are clearly evident in Fig. 5. Sudden change in density due to supercritical state of the fluids causes local large density gradient magnitude regions as shown in Fig. 6. Existence of local peaks in density is one of the notable aspects of supercritical flows [7,50]. As shown on both x–y and y–z planes in Fig. 6, such peaks occur at the bottom of the shear layer and in the braid regions following the temperature field. Such features are also observed in LES results. In Fig. 7, the filtered density gradient magnitude field obtained from LES and filtered DNS data are compared at t = 50. As shown, both magnitude and location of large density gradient regions are predicted consistently by LES.
To ensure accuracy of both FD and MC solvers of the FDF, we evaluate consistency of their results, as performed in previous work [24–28]. In the scalar FDF, all scalar fields predicted by the solvers are compared instantaneously to establish this consistency. This is shown in Fig. 8 for the filtered mass fraction of CO2 and temperature at t = 50. As shown, there is a close agreement between these variables obtained from the FD and MC solvers. Similar agreement is obtained for all scalars at all times, which verifies the consistency of the FDF solvers. Based on this agreement, we can also conclude that the thermodynamic state of the system is accurately represented at the MC particle level.
The accuracy of the extended FDF is assessed by comparing its results with the filtered DNS data at 150 atm. Figure 9 shows the Reynolds-averaged filtered fields of the streamwise velocity, temperature, density, and compression factor. As shown, LES results generally show favorable agreement with the DNS data. The thickness of the layer is shown to be slightly underpredicted which is attributed to the underprediction of SGS dissipation by the dynamic SGS model in this flow. The real fluid behavior is evident from the compression factor profile exhibiting values considerably below one which justifies the use of real fluid models. Similar agreement is obtained for other variables which indicates the accuracy of the extended FDF methodology for LES of supercritical turbulent flows.
6 Concluding Remarks
The scalar FDF methodology is extended and applied for LES of turbulent flows under supercritical condition. This is accomplished by incorporating models to represent real fluid behavior within the LES transport equations and the system of SDEs comprising the modeled FDF transport. These models include the generalized description of heat and mass diffusion for binary mixture of fluids and real fluid thermodynamic relations. The modified FDF formulation allows consideration of separate molecular transport for heat and mass fractions. The solution procedure is based on a hybrid Eulerian, Lagrangian MC method in which a MC particle method is used to solve the set of SDEs. It is shown that the thermodynamic state of the system is represented accurately at the MC particle level. The extended FDF is applied for LES of a temporally developing mixing layer under supercritical condition involving transport of passive scalars. The results are assessed against the data generated by DNS of the same layer. The consistency of the extended FDF formulation is demonstrated. The agreement between the LES-FDF and DNS results is reasonable and the LES-FDF is shown to be capable of capturing the supercritical flow features as observed in the DNS.
The authors are indebted to Dr. Mehdi Safari for useful discussions and contributions to this study. This work is supported in part by Department of Energy under Grant No. DE-SC0017097. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.
Conflict of Interest
There are no conflicts of interest.