## Abstract

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.

## 1 Introduction

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 [13]. 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. [48]. 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 [47,913].

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. [35,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 [16]. 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 [1723] 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 [2428]. 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. [29] 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 [30]. 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.

We consider general transport variables, varying in space x = xi(i = 1, 2, 3) and time t: fluid density ρ(x, t), velocity vector u = ui(x, t) along xi direction, specific internal energy e(x, t), pressure p(x, t), and mass fractions $Y=Yα(x,t)(α=1,2,…,Ns)$, where Ns is the number of chemical species. For convenience, we present the scalar variables in a common form as the scalar array $ϕ≡[ϕ1,ϕ2,…,ϕNs,ϕσ]=$$[Y1,Y1,…,YNs,e]$, where σ = Ns + 1. The transport variables satisfy the compressible form of equations for conservation of mass, momentum, energy, and species mass fractions for a binary mixture of Newtonian fluids
$∂ρ∂t+∂ρui∂xi=0$
(1)
$∂ρui∂t+∂ρuiuj∂xj=−∂p∂xi+∂τij∂xj$
(2)
$∂ρϕα∂t+∂ρuiϕα∂xi=−∂Jiα∂xi+ρSα$
(3)
where $Jiα$ represents the molecular mass and heat diffusion flux vectors for α = 1,…,Ns and α = σ, respectively; $Sα$ denotes the chemical reaction source term for α = 1,…,Ns. The present study concerns binary mixture of two passive species; we thus set Ns = 2 and $Sα=0$ for α = 1, 2. The source term in the energy equation is $Sσ=1/ρ(τij(∂ui/∂xj)−p(∂ui/∂xi))$. For a Newtonian fluid, the viscous stress tensor τij is represented as τij = μ((∂ui/∂xj) + (∂uj/∂xi) − (2/3)(∂uk/∂xk)δij) where μ and δij are the fluid dynamic viscosity and Kronecker delta, respectively. The viscosity is calculated using the power-law model similar to that in Ref. [8].

### 2.2 Equation of State and Thermodynamic Properties.

The conservation equations are closed by the cubic Peng–Robinson (PR) equation of state
$p=RuTv¯−Bm−Amv¯2+2Bmv¯−Bm2$
(4)
where $v¯$ and Ru denote the molar specific volume and universal gas constant, respectively. The two mixture parameters are defined as
$Am=∑α∑βXαXβAαβ,Bm=∑αXαBα$
(5)
in which Xα denotes the mole fraction of species α. These parameters incorporate the mixing rules, AαβAαβ(T) and Bα as given by Refs. [8,31]. Real fluid thermodynamic properties are accounted for using the departure function which describes the deviation of a real fluid from the ideal gas state: $λ¯−λo¯$, where $λ¯$ is a thermodynamics property per unit mole of the mixture at a general (non-ideal) fluid state with pressure p and temperature T; $λo¯$ is the value of $λ¯$ at the same temperature and composition but at the ideal gas state with pressure p°. From $g¯=a¯+pv¯$, where $g¯$ and $a¯$ are the Gibbs and Helmholtz free energy per unit mole, respectively, we can write
$g¯−go¯=(a¯−ao¯)+pv¯−povo¯$
(6)
At ideal gas pressure po, the ideal gas law $povo¯=RuT$ holds. Therefore, the departure function $g¯−go¯$ can be written as
$g¯−go¯=(a¯−ao¯)+pv¯−RuT$
(7)
For a fluid mixture with fixed composition and constant T, we can write
$da¯=−s¯dT−pdv¯+∑αμαdXα=−pdv¯$
(8)
where $s¯$ and μα are the entropy per unit mole of the mixture and the chemical potential, respectively. Integrating this equation gives the departure function $a¯−ao¯$ as
$a¯−ao¯=−∫vo¯v¯pdv′$
(9)
The Gibbs free energy $go¯$ at the ideal gas pressure po is expressed as
$go¯=∑αXα(gαo¯+RuTlnXα)$
(10)
Substituting Eqs. (9) and (10) into Eq. (7), we obtain
$g¯(T,p,X)=−∫vo¯v¯p(v′,T,X)dv′+pv¯−RuT+∑αXα(go¯α+RuTlnXα)$
(11)
where $X≡Xα(x,t)(α=1,2,…,Ns)$ and p(v′, T, X) is obtained from Eq. (4). The enthalpy per unit mole of the mixture is calculated using $h¯=g¯−T(∂g¯/∂T)p,X$ and Eq. (11) as
$h¯(T,p,X)=−∫vo¯v¯pdv′+pv¯−RuT+go¯−T∂∂T(−∫vo¯v¯pdv′+pv¯−RuT+go¯)p,X$
(12)
Using $ho¯=go¯−T(∂go¯/∂T)p,X$, Eq. (12) is simplified as
$h¯(T,p,X)=ho¯−∫vo¯v¯pdv′+pv¯−RuT−T∂∂T(−∫vo¯v¯pdv′+pv¯−RuT)p,X$
(13)
The last term on the right-hand side (RHS) can be expanded as
$T∂∂T(∫vo¯v¯pdv′)p,X=T∫vo¯v¯∂p∂Tdv′+pT(∂v¯∂T)p,X−poT(∂vo¯∂T)p,X$
(14)
Substituting Eqs. (14) and (4) into Eq. (13) results in
$h¯=ho¯+pv¯−RuT+(Am−T∂Am∂T)K1$
(15)
where K1 can be calculated using the following relation [8]:
$K1=122Bmln[v¯+(1−2)Bmv¯+(1+2)Bm]$
(16)
Using $e¯=h¯−pv¯$, the internal energy per unit mole of the mixture is
$e¯=eo¯+(Am−T∂Am∂T)K1$
(17)
where $eo¯$ is molar internal energy at the ideal gas state. This equation can be used to calculate molar specific heat at constant volume as
$c¯v≡(∂e¯∂T)v¯,X=c¯vo−T∂2Am∂T2K1$
(18)
where $c¯vo$ is the corresponding ideal gas molar specific heat. To obtain the molar specific heat at constant pressure $c¯p$ from Eq. (18), we first consider the following thermodynamic relations:
$ds¯=(∂s¯∂T)v¯,XdT+(∂s¯∂v¯)T,Xdv¯$
(19)
$dv¯=(∂v¯∂T)p,XdT+(∂v¯∂p)T,Xdp$
(20)
Combining these two equations yields
$(∂s¯∂T)p,X=(∂s¯∂T)v¯,X+(∂s¯∂v¯)T,X(∂v¯∂T)p,X$
(21)
Using $c¯p≡T(∂s¯/∂T)p,X$, $c¯v≡T(∂s¯/∂T)v¯,X$, and Eq. (21), we obtain
$c¯p−c¯v=T(∂s¯∂v¯)T,X(∂v¯∂T)p,X$
(22)
This equation can be further simplified using Maxwell’s relation $(∂s¯/∂v¯)T,X=(∂p/∂T)v¯,X$ and $dp=(∂p/∂T)v¯,XdT+(∂p/∂v¯)T,Xdv¯$. We thus have
$c¯p−c¯v=−T(∂p/∂T)v¯,X2(∂p/∂v¯)T,X$
(23)
Using this equation along with Eq. (18), the molar specific heat at constant pressure $c¯p=(∂h¯/∂T)p,X$ can be obtained as
$c¯p=c¯po−T(∂p/∂T)v¯,X2(∂p/∂v¯)T,X−Ru−T∂2Am∂T2K1$
(24)
where $c¯po$ is the corresponding ideal gas molar specific heat. In this equation, $(∂p/∂T)v¯,X$ and $(∂p/∂v¯)T,X$ are calculated using Eq. (4). The speed of sound as is calculated as $as≡(∂p/∂ρ)s=1/ρκs,$ where κs is the isentropic compressibility expressed as [8]: $κs=κT−v¯Tαv2/c¯p$. The expansivity αv and the isothermal compressibility κT in this equation are obtained from
$αv=−1v¯(∂p/∂T)v¯,X(∂p/∂v¯)T,X,κT=−1v¯1(∂p/∂v¯)T,X$
(25)
The compression factor is Z = pM/ρRuT where $M=∑α=1NsXαMα$ is the mixture molecular weight and $Mα$ denotes the molecular weight for species α.

### 2.3 Heat and Mass Transport.

The heat and mass diffusion fluxes for binary mixtures of real fluids are expressed as linear combinations of terms corresponding to gradients of temperature, chemical species concentrations, and pressure as discussed in Refs. [8,32]. The heat diffusion is described as the Irwing–Kirkwood form of the heat flux
$Jiσ=−[λIK∂T∂xi+αIKRuT(M∏η=1NsMη)Ji′1]$
(26)
The mass diffusion flux for methane is expressed as
$Ji1=−[Ji′1+αBK(∏η=1NsYη)ρDT∂T∂xi]$
(27)
where
$Ji′1=ρD[αD∂Y1∂xi+∏η=1NsYηRuT(∏η=1NsMηM)(V1¯M1−V2¯M2)∂p∂xi]$
(28)
In these equations, indices α = 1, 2 refer to methane and carbon dioxide, respectively. D is the binary diffusion coefficient, αD is the mass diffusion factor, and λIK is the thermal conductivity. The coefficients αIK and αBK are the thermal diffusion factors corresponding to the Irwing–Kirkwood and Bearman–Kirkwood forms of the heat flux, respectively, and $Vα¯=∂v¯/∂Xα$ denotes the partial molar volume. These parameters are calculated as discussed in Ref. [8].

### 2.4 Large Eddy Simulation Equations.

In LES, we consider the spatial filtering operation [22] $⟨Q(x,t)⟩=∫−∞+∞Q(x′,t)G(x′−x)dx′,$ where $G$ denotes the filter function of characteristic width Δ, and 〈Q(x, t)〉 represents the filtered value of the transport variable Q(x, t). In reacting flows, it is convenient to consider Favre filtered quantities as 〈Q(x, t)〉L = 〈ρQ〉/〈ρ〉. Applying the filtering operation to conservation equations, Eqs. (1)(3), results in Eqs. (29)(31) which are the transport equations for the filtered variables. For convenience in these equations we define θ ≡ [ρ, u, ϕ] and 〈θ〉 ≡ [〈ρ〉, 〈uL, 〈ϕL]. The second-order SGS correlations are defined as τL(a, b) ≡ 〈abL − 〈aLbL.
$∂⟨ρ⟩∂t+∂⟨ρ⟩⟨ui⟩L∂xi=0$
(29)
$∂⟨ρ⟩⟨ui⟩L∂t+∂⟨ρ⟩⟨ui⟩L⟨uj⟩L∂xj=−∂p(⟨θ⟩)∂xi+∂τij(⟨θ⟩)∂xj−∂⟨ρ⟩τL(ui,uj)∂xj−∂∂xi[⟨p(θ)⟩−p(⟨θ⟩)]+∂∂xj[⟨τij(θ)⟩−τij(⟨θ⟩)](30)$
(30)
$∂⟨ρ⟩⟨ϕα⟩L∂t+∂⟨ρ⟩⟨ui⟩L⟨ϕα⟩L∂xi=−∂Jiα(⟨θ⟩)∂xi−∂⟨ρ⟩τL(ui,ϕα)∂xi−∂∂xi[⟨Jiα(θ)⟩−Jiα(⟨θ⟩)]+⟨ρ⟩⟨Sα⟩L$
(31)
The unclosed SGS terms in these equations include the second-order correlations corresponding to SGS stress τL(ui, uj) in Eq. (30) and the SGS scalar fluxes of heat and mass $τL(ui,ϕα)$ in Eq. (31). The dynamic SGS model [33] is used to close these terms. Besides the second-order SGS moments, in Eq. (31), the source terms due to chemical reaction and that in the energy equation also appear as unclosed; these terms are accounted for using the FDF as explained in the next section. Aside from the aforementioned unclosed SGS effects, additional SGS terms appear in high-p flows, shown in square brackets in Eqs. (30) and (31), due to the nonlinear nature of the real fluid equation of state and generalized diffusion expressions. These high-p SGS terms are neglected similar to several other LES studies [15,34,35]. Among these terms, however, the dominant effects have been subject to SGS modeling in a few investigations [4,5,7]. These are the SGS pressure [〈p(θ)〉 − p(〈θ〉)] in Eq. (30) and the SGS heat and mass diffusion fluxes $[⟨Jiα(θ)⟩−Jiα(⟨θ⟩)]$ in Eq. (31). For discussions regarding the importance of these terms, we refer to Refs. [4,5,7,9,10,19,36].

### 2.5 The Filtered Density Function Transport Equation.

The FDF formulation considered in this study is that of the scalar filtered mass density function [25]
$FL(ψ,x;t)=∫−∞+∞ρ(x′,t)ζ(ψ,ϕ(x′,t))G(x′−x)dx′$
(32)
where ψ is the sample space variable vector corresponding to the scalar array ϕ. This equation defines the FDF as the mass-averaged spatially filtered value of fine-grained density $ζ(ψ,ϕ(x,t))≡∏α=1σδ(ψα−ϕα(x,t))$ where δ denotes the Dirac delta function. Using the scalar FDF, the filtered value of any function Q(x, t) can be obtained as $⟨ρQ(x,t)⟩=∫−∞+∞⟨Q|ψ⟩FL(ψ,x;t)dψ,$ where denotes the conditional filtered quantities [24,25]. As detailed in Ref. [25], by taking the time derivative of the fine-grained density function, we obtain the exact scalar FDF transport equation
$∂FL∂t+∂[⟨ui|ψ⟩FL]∂xi=∑α=1σ∂∂ψα[⟨1ρ∂Jiα∂xi|ψ⟩FL]−∑α=1σ∂∂ψα[⟨Sα|ψ⟩FL]$
(33)
In this equation, the chemical reaction source term (the last term on the RHS for α = 1,…, Ns) appears in a closed form but all other conditional filtered quantities are unclosed. The source term in the energy equation $Sσ$ is included in a similar fashion as in Ref. [30]. As proposed in Ref. [25], the second term on the left-hand side is modeled using a gradient diffusion model. The first term on the RHS corresponds to heat and mass diffusion which is divided into resolved and SGS parts. The former is accounted for explicitly while the latter is modeled using the linear mean-square estimation SGS mixing model [37,38]. Employing these models, the modeled FDF transport equation is expressed as
$∂FL∂t+∂[⟨ui⟩LFL]∂xi=∂∂xi[γt∂(FL/⟨ρ⟩)∂xi]+∑α=1σ∂∂ψα[1⟨ρ⟩∂⟨Jiα⟩∂xiFL]+∑α=1σ∂∂ψα[Cϕω(ψα−⟨ϕα⟩L)−⟨Sα|ψ⟩]FL$
(34)
where γt = 〈ρνt/Sct is the SGS diffusion coefficient with νt and Sct denoting the SGS kinematic viscosity and Schmidt number, respectively. The former is determined by the SGS model mentioned in Sec. 2.4 and the latter is set to Sct = 1 for all scalars following the previous study [25]. The parameter ω is the turbulent mixing frequency within the SGS modeled as ω = (γ + γt)/(〈ρ〉Δ2) [24,25] where γ denotes the molecular diffusion coefficient calculated as γ = μ/Sc, i.e., similar diffusivity is considered for all scalars within the SGS; the Schmidt number is expressed as Sc = μ/ρDαD [8]. The model constant $Cϕ$ is set to $Cϕ=1$ following the previous work [25]. In this equation, the first term on the RHS represents the turbulent diffusion of the FDF in the physical space. In the original scalar FDF formulation [24,25], this term also includes the molecular diffusion with the same diffusivity for all scalars. In order to employ the generalized diffusion models, the molecular transport is eliminated from this term and added as the second term on the RHS. This formulation is proposed in Ref. [29] along with Fickian diffusion to introduce the differential diffusion effect into the FDF framework. In this study, a similar approach is taken together with the generalized forms of heat and mass diffusion, described in Sec. 2.3, which allows more realistic representation of molecular transport in real fluids.
Moments of Eq. (34) provide the transport equations implied by the modeled FDF formulation. Comparing these transport equations with the exact ones, derived by directly filtering the conservation equations (Eqs. (1)(3)), reveals the implied closure by the FDF. This is demonstrated here for the first two moments. We take the first scalar moments of Eq. (34) resulting in
$∂⟨ρ⟩⟨ϕα⟩L∂t+∂⟨ρ⟩⟨ui⟩L⟨ϕα⟩L∂xi=−∂⟨Jiα⟩∂xi+∂∂xi(γt∂⟨ϕα⟩L∂xi)+⟨ρ⟩⟨Sα⟩L$
(35)
This equation is identical to Eq. (31) with the SGS heat and scalar fluxes represented by a gradient diffusion model. The filtered chemical reaction source term results from the FDF and thus appears in a closed form. The source term in the energy equation includes the pressure dilatation and viscous dissipation terms as described in Sec. 2.1. For the second moment, we consider the second-order SGS covariance $τL(ϕα,ϕβ)$ between any two scalars governed by the following exact equation derived using the filtering operation:
$∂⟨ρ⟩τL(ϕα,ϕβ)∂t+∂⟨ρ⟩⟨ui⟩LτL(ϕα,ϕβ)∂xi=−∂⟨ρ⟩τL(ui,ϕα,ϕβ)∂xi−⟨ρ⟩(τL(ϕα,ui)∂⟨ϕβ⟩L∂xi+τL(ϕβ,ui)∂⟨ϕα⟩L∂xi)−[⟨ϕα∂Jiβ∂xi⟩−⟨ϕα⟩L∂⟨Jiβ⟩∂xi]−[⟨ϕβ∂Jiα∂xi⟩−⟨ϕβ⟩L∂⟨Jiα⟩∂xi]+⟨ρ⟩(τL(ϕα,Sβ)+τL(ϕβ,Sα))$
(36)
where the third-order correlation is τL (a, b, c) = 〈a b cL − 〈aLτL (b, c) − 〈bLτL (a, c) − 〈cLτL (a, b) − 〈aLbLcL. A similar transport equation results from the FDF by taking the second moment of Eq. (34) as
$∂⟨ρ⟩τL(ϕα,ϕβ)∂t+∂⟨ρ⟩⟨ui⟩LτL(ϕα,ϕβ)∂xi=∂∂xi(γt∂τL(ϕα,ϕβ)∂xi)+2γt∂⟨ϕα⟩L∂xi∂⟨ϕβ⟩L∂xi−2Cϕω⟨ρ⟩τL(ϕα,ϕβ)+⟨ρ⟩[τL(ϕα,Sβ)+τL(ϕβ,Sα)]$
(37)
This equation is similar to Eq. (36) with the third-order correlation and the SGS heat and mass fluxes (the first three terms on the RHS of Eq. (36)) modeled using the gradient diffusion closure. This results in the first two terms on the RHS of Eq. (37). The correlations involving molecular heat and mass diffusion in Eq. (36) (the terms in the square brackets) contain both molecular transport and scalar dissipation effects. The latter is modeled by the third term on the RHS of Eq. (37). However, there are no terms corresponding to the former in Eq. (37), as also indicated in Ref. [29]. The correlations involving source terms on the last two terms are obtained from the FDF. This comparison shows the implied closure by the FDF and that the transport equations implied by the FDF are consistent with those of the original LES.

### 2.6 Monte Carlo Solution of the Filtered Density Function.

To solve the FDF transport equation, a system of SDEs representing the FDF is considered [24,25]. These SDEs are solved by a Lagrangian Monte Carlo (MC) method in which the MC particles undergo motion in physical space due to filtered velocity and SGS diffusion. The thermochemical state of particles is changed due to SGS mixing, molecular diffusion, as well as source terms, including the chemical reaction effects. These are described by the following system of SDEs:
$dxi+=[⟨ui⟩L+1⟨ρ⟩∂γt∂xi]dt+2γt⟨ρ⟩dWi$
(38)
$dϕα+=−Cϕω[ϕα+−⟨ϕα⟩L]dt−1⟨ρ⟩∂⟨Jiα⟩∂xidt+Sαdt$
(39)
where α = 1,…, σ. Variables $xi+$ and $ϕα+$ denote the stochastic processes for Lagrangian position and scalars, respectively. The term Wi denotes the Wiener process. Compared to the original scalar FDF formulation [24,25], the molecular diffusion coefficient is removed from Eq. (38) and instead the resolved molecular transport is included in Eq. (39) as a drift term to incorporate the generalized diffusion models. A similar approach is used in Refs. [29,39] along with Fickian diffusion. The source term in the internal energy SDE includes the filtered form of pressure dilatation and viscous dissipation effects in a similar way as in Ref. [30]. The Fokker–Planck equation corresponding to this set of SDEs is identical to the modeled FDF transport equation (Eq. (34)). Therefore, the joint scalar PDF resulting from the SDEs corresponds to the FDF at the level of one-point, one-time statistics [25,28].

## 3 Numerical Solution Procedure

Numerical implementation of the FDF consists of a hybrid finite difference (FD)/MC method as considered in previous work [2428]. 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 $xi(n)(t)$ and scalars $ϕα(n)(t)$, where n = 1,…, Np. These particles are evolved in time by integrating the set of SDEs discretized using the Euler–Maruyamma discretization [40]. 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 [8]. 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 ≤ xL, −(L/2) ≤ yL/2, 0 ≤ zL 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 [2427], 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 [44] 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 $⟨⟩L¯$, respectively.

## 5 Results

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.

Fig. 1
Fig. 1
Close modal

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) [49] 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 $10%$, 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.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal

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 xy plane along with secondary structures on the yz 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 [7] 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 xy and yz 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.

Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal
Fig. 7
Fig. 7
Close modal

To ensure accuracy of both FD and MC solvers of the FDF, we evaluate consistency of their results, as performed in previous work [2428]. 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.

Fig. 8
Fig. 8
Close modal

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.

Fig. 9
Fig. 9
Close modal

## 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.

## Acknowledgment

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.

## References

1.
Wright
,
S. A.
,
,
R. F.
,
Vernon
,
M. E.
,
Rochau
,
G. E.
, and
Pickard
,
P. S.
,
2010
, “
Operation and Analysis of a Supercritical CO2 Brayton Cycle
,”
Sandia National Laboratories
,
Livermore, CA
.
2.
Vesely
,
L.
,
Manikantachari
,
K. R. V.
,
Vasu
,
S.
,
Kapat
,
J.
,
Dostal
,
V.
, and
Martin
,
S.
,
2019
, “
Effect of Impurities on Compressor and Cooler in Supercritical CO2 Cycles
,”
ASME J. Energy Resour. Technol.
,
141
(
1
), p.
012003
.
3.
Strakey
,
P. A.
,
2019
, “
Oxy-Combustion Modeling for Direct-Fired Supercritical CO2 Power Cycles
,”
ASME J. Energy Resour. Technol.
,
141
(
7
), p.
070706
.
4.
Taşkinoğlu
,
E. S.
, and
Bellan
,
J.
,
2011
, “
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
, pp.
156
193
.
5.
Borghesi
,
G.
, and
Bellan
,
J.
,
2015
, “
A Priori and a Posteriori Investigations for Developing Large Eddy Simulations of Multi-Species Turbulent Mixing Under High-Pressure Conditions
,”
Phys. Fluids
,
27
(
3
), p.
035117-1
035117-35
.
6.
Masi
,
E.
, and
Bellan
,
J.
,
2011
, “
The Subgrid-Scale Scalar Variance Under Supercritical Pressure Conditions
,”
Phys. Fluids
,
23
(
8
), p.
085101
, 1–22.
7.
Selle
,
L. C.
,
Okong’o
,
N. A.
,
Bellan
,
J.
, and
,
K. G.
,
2007
, “
Modelling of Subgrid-scale Phenomena in Supercritical Transitional Mixing Layers: An a Priori Study
,”
J. Fluid Mech.
,
593
, pp.
57
91
.
8.
Miller
,
R. S.
,
,
K. G.
, and
Bellan
,
J.
,
2001
, “
Direct Numerical Simulations of Supercritical Fluid Mixing Layers Applied to Heptane–Nitrogen
,”
J. Fluid Mech.
,
436
, pp.
1
39
.
9.
Foster
,
J.
, and
Miller
,
R. S.
,
2010
, “Fundamentals of High Pressure Combustion,”
High Pressure Processes in Chemical Engineering
,
M.
Lackner
, ed.,
ProcessEng Engineering GmbH
,
Vienna
, pp.
53
75
.
10.
Foster
,
J.
, and
Miller
,
R. S.
,
2012
, “
A Priori Analysis of Subgrid Mass Diffusion Vectors in High Pressure Turbulent Hydrogen/Oxygen Reacting Shear Layer Flames
,”
Phys. Fluids
,
24
(
7
), p.
075114
, 1–26.
11.
Guven
,
U.
, and
Ribert
,
G.
,
2019
, “
Impact of Non-Ideal Transport Modeling on Supercritical Flow Simulation
,”
Proc. Combust. Inst.
,
37
(
3
), pp.
3255
3262
.
12.
Park
,
S.
,
Urso
,
J.
,
Manikantachari
,
K.
,
,
A.
,
Zambon
,
A.
, and
Vasu
,
S.
,
2020
, “
Measurements of Density and Sound Speed in Mixtures Relevant to Supercritical CO2 Cycles
,”
ASME J. Energy Res. Technol.
,
142
(
10
), p.
102105
.
13.
Manikantachari
,
K.
,
Rahman
,
R.
,
Martin
,
S.
,
Velez
,
C.
, and
Vasum
,
S.
,
2021
, “
Influence of Equation-of-States on Supercritical CO2 Combustion Mixtures
,”
ASME J. Energy Res. Technol.
,
143
(
6
), p.
062106
.
14.
Schmitt
,
T.
,
Méry
,
Y.
,
Boileau
,
M.
, and
Candel
,
S.
,
2011
, “
Large-Eddy Simulation of Oxygen/Methane Flames Under Transcritical Conditions
,”
Proc. Combust. Inst.
,
33
(
1
), pp.
1383
1390
.
15.
Schmitt
,
T.
,
Selle
,
L.
,
Cuenot
,
B.
, and
Poinsot
,
T.
,
2009
, “
Large-Eddy Simulation of Transcritical Flows
,”
Comp. Rendus Mecanique
,
337
(
6
), pp.
528
538
.
16.
Ribert
,
G.
,
Petit
,
X.
, and
Domingo
,
P.
,
2017
, “
High-Pressure Methane-Oxygen Flames. Analysis of Sub-Grid Scale Contributions in Filtered Equations of State
,”
J. Supercrit. Fluids
,
121
, pp.
78
88
.
17.
Pope
,
S. B.
,
2013
, “
Small Scales, Many Species and the Manifold Challenges of Turbulent Combustion
,”
Proc. Combust. Inst.
,
34
(
1
), pp.
1
31
.
18.
Haworth
,
D. C.
,
2010
, “
Progress in Probability Density Function Methods for Turbulent Reacting Flows
,”
Prog. Energy Combust. Sci.
,
36
(
2
), pp.
168
259
.
19.
Miller
,
R. S.
, and
Foster
,
J. W.
,
2016
, “
Survey of Turbulent Combustion Models for Large-Eddy Simulations of Propulsive Flowfields
,”
AIAA J.
,
54
(
10
), pp.
2930
2946
.
20.
Ansari
,
N.
,
Jaberi
,
F. A.
,
Sheikhi
,
M. R. H.
, and
Givi
,
P.
,
2011
, “Filtered Density Function as a Modern CFD Tool,”
Engineering Applications of CFD
,
R. S.
Maher
, ed., Vol.
1
,
International Energy and Environment Foundation
,
Fluid Mechanics and Its Applications
Al-Najaf
, pp.
1
22
, Chapter 1.
21.
Givi
,
P.
,
2006
, “
Filtered Density Function for Subgrid Scale Modeling of Turbulent Combustion
,”
AIAA J.
,
44
(
1
), pp.
16
23
.
22.
Pope
,
S. B.
,
2000
,
Turbulent Flows
,
Cambridge University Press
,
Cambridge
.
23.
Sammak
,
S.
,
Ren
,
Z.
, and
Givi
,
P.
,
2020
, “Modern Developments in Filtered Density Function,”
Modeling and Simulation of Turbulent Mixing and Reaction
,
D.
Livescu
,
A.
Nouri
,
F.
Battaglia
, and
P.
Givi
, eds.,
Heat and Mass Transfer
,
Springer
,
Singapore
, pp.
181
200
.
24.
Colucci
,
P. J.
,
Jaberi
,
F. A.
,
Givi
,
P.
, and
Pope
,
S. B.
,
1998
, “
Filtered Density Function for Large Eddy Simulation of Turbulent Reacting Flows
,”
Phys. Fluids
,
10
(
2
), pp.
499
515
.
25.
Jaberi
,
F. A.
,
Colucci
,
P. J.
,
James
,
S.
,
Givi
,
P.
, and
Pope
,
S. B.
,
1999
, “
Filtered Mass Density Function for Large Eddy Simulation of Turbulent Reacting Flows
,”
J. Fluid Mech.
,
401
, pp.
85
121
.
26.
Gicquel
,
L. Y. M.
,
Givi
,
P.
,
Jaberi
,
F. A.
, and
Pope
,
S. B.
,
2002
, “
Velocity Filtered Density Function for Large Eddy Simulation of Turbulent Flows
,”
Phys. Fluids
,
14
(
3
), pp.
1196
1213
.
27.
Sheikhi
,
M. R. H.
,
Drozda
,
T. G.
,
Givi
,
P.
, and
Pope
,
S. B.
,
2003
, “
Velocity-Scalar Filtered Density Function for Large Eddy Simulation of Turbulent Flows
,”
Phys. Fluids
,
15
(
8
), pp.
2321
2337
.
28.
Sheikhi
,
M. R. H.
,
Givi
,
P.
, and
Pope
,
S. B.
,
2007
, “
Velocity-Scalar Filtered Mass Density Function for Large Eddy Simulation of Turbulent Reacting Flows
,”
Phys. Fluids
,
19
(
9
), p.
095106
.
29.
McDermott
,
R.
, and
Pope
,
S. B.
,
2007
, “
A Particle Formulation for Treating Differential Diffusion in Filtered Density Function Methods
,”
J. Comput. Phys.
,
226
(
1
), pp.
947
993
.
30.
,
A.
,
Li
,
Z.
, and
Jaberi
,
F. A.
,
2011
, “
Compressible Scalar Filtered Mass Density Function Model for High-Speed Turbulent Flows
,”
AIAA J.
,
49
(
10
), pp.
2130
2143
.
31.
,
K. G.
,
Miller
,
R. S.
, and
Bellan
,
J.
,
1997
, “
Efficient High Pressure State Equations
,”
AIChE J.
,
43
(
6
), pp.
1605
1610
.
32.
Okong’O
,
N.
, and
Bellan
,
J.
,
2002
, “
Direct Numerical Simulation of a Transitional Supercritical Binary Mixing Layer: Heptane and Nitrogen
,”
J. Fluid Mech.
,
464
, pp.
1
34
.
33.
Moin
,
P.
,
Squires
,
W.
,
Cabot
,
W. H.
, and
Lee
,
S.
,
1991
, “
A Dynamic Subgrid-Scale Model for Compressible Turbulence and Scalar Transport
,”
Phys. Fluids A
,
3
(
11
), pp.
2746
2757
.
34.
Selle
,
L.
, and
Schmitt
,
T.
,
2010
, “
Large-Eddy Simulation of Single-Species Flows Under Supercritical Thermodynamic Conditions
,”
Combust. Sci. Technol.
,
182
(
4–6
), pp.
392
404
.
35.
Müller
,
H.
,
Niedermeier
,
C. A.
,
Matheis
,
J.
,
Pfitzner
,
M.
, and
Hickel
,
S.
,
2016
, “
Large-Eddy Simulation of Nitrogen Injection at Trans- and Supercritical Conditions
,”
Phys. Fluids
,
28
(
1
), p.
015102
.
36.
Ma
,
Z.
,
Korucu
,
A.
, and
Miller
,
R. S.
,
2015
, “
A Priori Analysis of Subgrid Scale Pressure and Heat Flux in High Pressure Mixing and Reacting Shear Layers
,”
Combust. Theor. Model.
,
19
(
6
), pp.
807
832
.
37.
Dopazo
,
C.
, and
O’Brien
,
E. E.
,
1976
, “
Statistical Treatment of Non-Isothermal Chemical Reactions in Turbulence
,”
Combust. Sci. Technol.
,
13
(
1–6
), pp.
99
112
.
38.
O’Brien
,
E. E.
,
1980
, “The Probability Density Function (PDF) Approach to Reacting Turbulent Flows,”
Turbulent Reacting Flows
,
P. A.
Libby
, and
F. A.
Williams
, eds., Vol.
44
,
Topics in Applied Physics
,
Springer-Verlag
,
Heidelberg
, pp.
185
218
, Chapter 5.
39.
Turkeri
,
H.
,
Zhao
,
X.
,
Pope
,
S. B.
, and
,
M.
,
2019
, “
Large Eddy Simulation/Probability Density Function Simulations of the Cambridge Turbulent Stratified Flame Series
,”
Combust. Flame
,
199
, pp.
24
45
.
40.
Kloeden
,
P. E.
,
Platen
,
E.
, and
Schurz
,
H.
,
1997
,
Numerical Solution of Stochastic Differential Equations Through Computer Experiments
, 2nd ed.,
Springer-Verlag
,
New York
.
41.
Carpenter
,
M. H.
,
1990
, “A High-Order Compact Numerical Algorithm for Supersonic Flows,”
Twelfth International Conference on Numerical Methods in Fluid Dynamics
, Vol.
371
, Lecture Notes in Physics,
K. W.
Morton
, ed.,
Springer-Verlag
,
New York
, pp.
254
258
.
42.
Kennedy
,
C. A.
, and
Carpenter
,
M. H.
,
1994
, “
Several New Numerical Methods for Compressible Shear-Layer Simulations
,”
Appl. Numer. Math.
,
14
(
4
), pp.
397
433
.
43.
Vreman
,
B.
,
Geurts
,
B.
, and
Kuerten
,
H.
,
1997
, “
Large-Eddy Simulation of the Turbulent Mixing Layer
,”
J. Fluid Mech.
,
339
, pp.
357
390
.
44.
Poinsot
,
T.
, and
Veynante
,
D.
,
2011
,
Theoretical and Numerical Combustion
, 3rd ed.,
R. T. Edwards Inc.
,
.
45.
Metcalfe
,
R. W.
,
Orszag
,
S. A.
,
Brachet
,
M. E.
,
Menon
,
S.
, and
Riley
,
J. J.
,
1987
, “
Secondary Instabilities of a Temporally Growing Mixing Layer
,”
J. Fluid Mech.
,
184
, pp.
207
243
.
46.
Lin
,
S. J.
, and
Corcos
,
G. M.
,
1984
, “
The Mixing Layer: Deterministic Models of a Turbulent Flow. Part 3. The Effect of Plane Strain on the Dynamics of Streamwise Vortices
,”
J. Fluid Mech.
,
141
, pp.
139
178
.
47.
Sandham
,
N. D.
, and
Reynolds
,
W. C.
,
1991
, “
Three-dimensional Simulations of Large Eddies in the Compressible Mixing Layer
,”
J. Fluid Mech.
,
224
, pp.
133
158
.
48.
Moser
,
R. D.
, and
Rogers
,
M. M.
,
1993
, “
Spanwise Scale Selection in Plane Mixing Layers
,”
J. Fluid Mech.
,
247
, pp.
321
337
.
49.
NIST Chemistry WebBook
,
2021
, http://webbook.nist.gov/chemistry
50.
Bellan
,
J.
,
2017
, “
Direct Numerical Simulation of a High-Pressure Turbulent Reacting Temporal Mixing Layer
,”
Combust. Flame
,
176
, pp.
245
262
.