Partially-Averaged Navier-Stokes Closure Modeling for Variable-Density Turbulent Flow
F.S. Pereira, F.F. Grinstein, D.M. Israel, R. Rauenzahn, S.S. Girimaji
PPartially-Averaged Navier-Stokes Closure Modeling for Variable-Density TurbulentFlow
F.S. Pereira , ∗ F.F. Grinstein , D.M. Israel , R. Rauenzahn , and S.S. Girimaji Los Alamos National Laboratory, X Division, Los Alamos, New Mexico, USA and Texas A&M University, Ocean Engineering, College Station, Texas, USA
We extend the framework of the Partially-Averaged Navier-Stokes (PANS) equations to variable-density flow, i.e., multi-material and/or compressible mixing problems with density variations andproduction of turbulence kinetic energy by both shear and buoyancy mechanisms. This includes a-priori testing to develop guidelines toward the efficient selection of the parameters controlling thephysical resolution (range of resolved scales) of PANS. The method is then evaluated predictingtwo archetypal test-cases involving transient turbulence, hydrodynamic instabilities, and coherentstructures: the Taylor-Green vortex (TGV) at Reynolds number Re = 3000, and the Rayleigh-Taylor (RT) flow at Atwood number 0 . max ≈ I. INTRODUCTION
The numerical prediction of variable-density (multi-material and/or compressible) flows is crucial to numer-ous applications of fundamental and applied fluid me-chanics – e.g., scramjets, mixing problems, oceanography,supernova explosions, and inertial confinement fusion. Inaddition to the shear production mechanism of constantdensity turbulence, variable-density flows include baro-clinic production due to local mean density and pressuregradients. These flows are inherently transient and alsoinclude complex physics involving dilatation effects, char-acteristic hydrodynamic instabilities and coherent struc-tures [1–6], and interactions between material and veloc-ity fields. All these aspects make modeling and simula-tion of variable-density flows rife with challenges.Direct numerical simulation (DNS) is the ideal optionfor prediction of any flow, since it is an exact solutionto the governing physical equations. Large-eddy simu-lation (LES) [7], since it still resolves most of the tur-bulence spectrum, is expected to lead to accurate repre-sentations of the flow dynamics. Yet, such high-fidelityscale-resolving simulations (SRS) come at with computa-tional expense that may be too prohibitive for practicalapplications. Also, selecting proper initial and boundaryconditions for DNS and LES of variable-density problemsis difficult.The Reynolds-averaged Navier-Stokes (RANS) equa-tions [8–12] are an alternative formulation to simulatepractical flows of variable-density. In contrast to DNSand LES, RANS relies on a fully statistical descriptionof turbulence, in which all turbulence scales are modeledthrough a constitutive relationship named turbulence clo-sure [13–16]. This modeling strategy entails computa- ∗ F.S. Pereira: [email protected] tion of the mean and coherent fields [17, 18]) and, con-sequently, significantly reduces the computations’ cost.However, RANS closures are usually inaccurate predict-ing many problems of interest in this study.The caveats of the former modeling strategies drovethe emergence of a new paradigm of SRS methods spe-cially designed to predict complex flows efficiently. Thisclass of models named bridging methods was proposedby Germano [19, 20] and Speziale [21] to seamlessly op-erate between DNS and RANS, and only resolve the flowscales not amenable to modeling. The remaining scalescan be accurately represented through an adequate clo-sure [22]. This strategy is responsible for the potentialefficiency (accuracy vs. cost) of bridging methods, andintroduces the idea of accuracy-on-demand. Very Large-Eddy Simulation (VLES) [21], Limited Numerical Scales(LNS) [23], Flow Simulation Methodology (FSM) [24],Partially-Integrated Transport Model (PITM) [25, 26],and Partially-Averaged Navier-Stokes (PANS) equations[27] are examples of bridging formulations.Despite being widely used in many scientific areas,bridging models are still not common in the variable-density flow community. There are three main factorscontributing to this outcome: i ) complexity - bridging closures are typically basedon one-point RANS models, which are calibratedfor total turbulent quantities. However, sincebridging methods can operate at any range ofresolved scales, the inability to reliably estimateRANS variables may lead to calibration deficienciesof the closure [28]. This is expected to be particu-larly relevant for low physical resolutions, transientand transitional flows, and second-moment closures[13–15, 29, 30]. Further, the RANS closures forvariable-density flow possess terms to account fordensity variations that can be difficult to extend tothe bridging framework. a r X i v : . [ phy s i c s . f l u - dyn ] F e b ii ) physical resolution - bridging models normallyrequire a physical resolution parameter which de-termines the range of resolved scales and, conse-quently, the computational efficiency. Whereas ex-cessively large physical resolutions increase the sim-ulations cost unnecessarily, low values of such a pa-rameter can lead to inaccurate computations. Thismakes the selection of the physical resolution cru-cial to the simulations accuracy. A closure forvariable-density flow will rely on multiple modelevolution equations [14, 15]. This raises the ques-tion, how do different modeled turbulence quanti-ties behave with the physical resolution. I.e., whatis the fraction of each dependent turbulence quan-tity being modeled for a given range of resolvedscales. VLES, LNS, and FSM formulations use apragmatic approach where the magnitude of theReynolds-stress tensor is scaled by a given factor.Despite being simple, this strategy does not yieldthe correct fixed point behavior for the closure sys-tem [31]. PANS and PITM, on the other hand,do not experience the these issues since they useparameters to define the modeled-to-total ratio ofeach turbulence dependent quantity. iii ) commutation errors - SRS formulations arebased upon the scale-invariance property of theNavier-Stokes equations [19]. This property is re-sponsible for the formal similarity between the fil-tered Navier-Stokes and RANS equations, and re-quires the implicit model filter operator to com-mute with temporal and spatial differentiation [19].Most SRS computations are conducted with a spa-tially varying physical resolution to optimize theuse of the grid. Even for conventional LES models,this raises potential modeling issues (as shown inHamba [32]). For bridging models, such as PANS,spatial variation of the resolution control param-eters lead to additional commutation terms (Giri-maji and Wallin [33]), that can reach the magni-tude of the convective terms of the filtered Navier-Stokes equations. This is expected to be relevantfor variable-density turbulent flows due to theirtransient and transitional nature.This work proposes a new PANS bridging model forvariable-density turbulent flow. Specifically, the PANSframework of Girimaji [27] and Suman and Girimaji [34]is extended to variable-density flow by deriving the PANSversion of the six-equation BHR-LEVM (linear eddy vis-cosity model) closure [29, 35]. We investigate the effectof the range of resolved scales on PANS turbulent de-pendent quantities through a-priori testing. The accu-racy of PANS BHR-LEVM is evaluated using two bench-mark problems: the Taylor-Green vortex (TGV) [36] atReynolds number (Re) 3000 and initial Mach number(Ma) 0 .
28, and the Rayleigh-Taylor flow [4, 37] at At-wood number (At) 0 .
50, (Re) max ≈ < . II. GOVERNING EQUATIONS
The partially-averaged Navier-Stokes (PANS) equa-tions are based upon the scale-invariance property of theNavier-Stokes equations. This property has been demon-strated by Germano [19] for incompressible flow, and ex-tended to compressible flow by Suman and Girimaji [34].To derive the PANS equations for variable-density flow(multi-material and/or compressible), let us start by con-sidering a general linear and constant preserving filteringoperator (cid:104) · (cid:105) , (cid:104) Φ + Φ (cid:105) = (cid:104) Φ (cid:105) + (cid:104) Φ (cid:105) , (1) (cid:104) α Φ (cid:105) = α (cid:104) Φ (cid:105) , (2)where Φ is a generic variable, and α is a constant. Thisfilter commutes with spatial and temporal differentiationso that (cid:28) ∂ Φ ∂x i (cid:29) = ∂ (cid:104) Φ (cid:105) ∂x i , (cid:28) ∂ Φ ∂t (cid:29) = ∂ (cid:104) Φ (cid:105) ∂t , (3)and decomposes any instantaneous flow quantity Φ intoa filtered (resolved), (cid:104) Φ (cid:105) , and modeled (unresolved), φ ,component, Φ ≡ (cid:104) Φ (cid:105) + φ . (4)This decomposition can be extended to variable-densityflow through the concept of Favre-averaging [9–12],Φ ≡ { Φ } + φ ∗ , (5)where { Φ } and φ ∗ are the density-weighted filtered, { Φ } ≡ (cid:104) ρ Φ (cid:105) / (cid:104) ρ (cid:105) , and modeled fluctuating, φ ∗ = Φ − { Φ } ,components of Φ. In the limit of all turbulence scales be-ing modeled, the former decompositions are equivalentto Reynolds- [8] and Favre-averaging [9–12],Φ = Φ + φ (cid:48) , (6)Φ = ˜Φ + φ (cid:48)(cid:48) , (7)Φ and φ (cid:48) being the (time, ensemble or spatial) averagedand turbulent components of Φ, whereas ˜Φ and φ (cid:48)(cid:48) thedensity-weighted averaged and fluctuating counterpartsof Φ.The application of such filtering operators to theconservation equations for mass, momentum, total en-ergy, and fluid species [40, 41] leads to the filtered orpartially-averaged form of the Navier-Stokes equationsfor variable-density flow [19, 27, 34, 42], ∂ (cid:104) ρ (cid:105) ∂t + ∂ ( (cid:104) ρ (cid:105){ V i } ) ∂x i = 0 , (8) ∂ ( (cid:104) ρ (cid:105){ V i } ) ∂t + ∂ ( (cid:104) ρ (cid:105){ V j }{ V i } ) ∂x j = − ∂ (cid:104) P (cid:105) ∂x i + ∂ (cid:104) σ ij (cid:105) ∂x j + ∂ (cid:0) (cid:104) ρ (cid:105) τ ( V i , V j ) (cid:1) ∂x j + (cid:104) ρ (cid:105) g i , (9) ∂ ( (cid:104) ρ (cid:105){ E } ) ∂t + ∂ ( (cid:104) ρ (cid:105){ E }{ V j } ) ∂x j = − ∂ (cid:0) (cid:104) ρ (cid:105) τ ( V j , E ) (cid:1) ∂x j − ∂ ( { V j }(cid:104) P (cid:105) ) ∂x j − ∂τ ( V j , P ) ∂x j + ∂ ( { V i }(cid:104) σ ij (cid:105) ) ∂x j + ∂τ ( V i , σ ij ) ∂x j − ∂ (cid:104) q cj (cid:105) ∂x j − ∂ (cid:104) q hj (cid:105) ∂x j , (10) ∂ ( (cid:104) ρ (cid:105){ c n } ) ∂t + ∂ ( (cid:104) ρ (cid:105){ c n }{ V j } ) ∂x j = − ∂ (cid:104) J nj (cid:105) ∂x j . (11)Here, t is the time, x i are the coordinates of a Cartesiansystem, ρ is the fluid density, V i are the Cartesian velocitycomponents, P is the pressure, σ ij is the viscous-stresstensor assuming Newtonian fluid, (cid:104) σ ij (cid:105) = 2 µ (cid:18) { S ij } − ∂ { V k } ∂x k (cid:19) , (12) { S ij } is the resolved strain-rate tensor, { S ij } = 12 (cid:18) ∂ { V i } ∂x j + ∂ { V j } ∂x i (cid:19) , (13) µ is the fluid’s dynamic viscosity, g i is the gravitationalacceleration vector, E = V i + e is the total energy of the fluid, e is the internal energy, c n is the mass concen-tration of material n , q c is the conductive heat flux, q d is the interdiffusional enthalpy flux, and J n is the massfraction diffusivity flux of material n . Also, τ (Φ i , Φ j )and τ (Φ i , Φ j ) are generalized central second momentswhich account for the effect of the modeled turbulence inthe resolved flow field. Expressing the PANS equations interms of generalized central second moments [19] guaran-tees scale-invariance. Such tensors are formally definedas [19, 34] τ (Φ i , Φ j ) ≡ { Φ i Φ j } − { Φ i }{ Φ j } , (14) τ (Φ i , Φ j ) ≡ (cid:104) Φ i Φ j (cid:105) − { Φ i }(cid:104) Φ j (cid:105) . (15)In equations 8 to 11, the pressure is calculated assuminga thermally perfect gas ( P = ρRT ). Thus, its resolvedcomponent is given by [34], (cid:104) P (cid:105) = ( γ − (cid:104) ρ (cid:105) (cid:18) { E } − { V k }{ V k } − k u (cid:19) , (16)where T is the temperature, γ is the ratio between spe-cific heats, and k u is the unresolved or modeled specificturbulence kinetic energy.The generalized central second-moments and fluxespresent in the PANS equations need modeling to closethe resultant system of equations. In the present work,this is accomplished through the Boussinesq approxima-tion [43], τ ( V i , V j ) = 2 ν u { S ij } − k u δ ij , (17)and the relationships given in Besnard et al. [13], Sumanand Girimaji [34], Stalsberg-Zarling and Gore [29], andSchwarzkopf et al. [30]. These lead to the partially-averaged form of the energy and fluid species equations, ∂ ( (cid:104) ρ (cid:105){ E } ) ∂t + ∂ ( (cid:104) ρ (cid:105){ E }{ V j } ) ∂x j = − ∂ ( { V j }(cid:104) P (cid:105) ) ∂x j + ∂ ( { V i }(cid:104) σ ij (cid:105) ) ∂x j + ∂ (cid:0) (cid:104) ρ (cid:105){ V i } τ ( V i , V j ) (cid:1) ∂x j + ∂∂x j (cid:20)(cid:18) µ + µ u σ k (cid:19) ∂k u ∂x j (cid:21) − ∂∂x j (cid:20) c p (cid:18) µ Pr + µ u Pr t (cid:19) ∂ (cid:104) T (cid:105) ∂x j (cid:21) − ∂∂x j (cid:34) n t (cid:88) n =1 h n (cid:104) J nj (cid:105) (cid:35) , (18) ∂ ( (cid:104) ρ (cid:105){ c n } ) ∂t + ∂ ( (cid:104) ρ (cid:105){ c n }{ V j } ) ∂x j = − ∂ (cid:104) J nj (cid:105) ∂x j = − ∂∂x j (cid:20) (cid:104) ρ (cid:105) (cid:18) D + µ u σ c (cid:19) ∂ { c n } ∂x j (cid:21) , (19)where ν u = µ u / (cid:104) ρ (cid:105) is the kinematic turbulent viscosityof the unresolved scales, σ k and σ c are turbulent diffu-sion coefficients, c v is the constant specific heat (idealgas is assumed), Pr is the Prandtl number, Pr t is theturbulent Prandtl number defined as Pr t = c v ν u /κ , κ is the effective thermal conductivity, and h n is the en-thalpy of material n . Throughout this manuscript, allmodeled/unresolved PANS and RANS turbulence quan-tities are denoted by the subscripts u and t , respectively.The relationships above create two additional turbu-lence quantities, k u and ν u , that need modeling. This isaccomplished through the BHR-LEVM [29, 35] closuremodel which is now derived for PANS. II.1. PANS BHR-LEVM closure
Most PANS closures are based on one-point, linear tur-bulent viscosity RANS closures. This modeling strategyis chosen to balance sufficient complexity to accuratelyoperate at any degree of physical resolution with scale-aware minor modifications [44], without the loss of ro-bustness observed for full Reynolds-stress closures. Forvariable-density flow we choose the BHR model originalproposed by [13], in k t − S t linear turbulent viscosity formfound in [29, 35]. The model requires transport equationsfor six turbulent dependent variables: the turbulence ki-netic energy, k t , turbulence dissipation length-scale, S t ,velocity mass flux, a i t = ρ (cid:48) v (cid:48) i /ρ , and density-specific vol-ume correlation, b t = − ρ (cid:48) (1 /ρ ) (cid:48) . For the b t equation, weuse the newer formulation of [45]. The subscript t indi-cates the total turbulence quantity, that is, the quantitypredicted by the RANS model that includes the action ofall the turbulent scales of motion. The subscript u willbe used for partial-averaged quantities, including onlythe unresolved portion of the turbulent scales.The RANS BHR-LEVM model used in this work cal-culates the total kinematic turbulent viscosity as, ν t = µ t ρ = c µ S t (cid:112) k t , (20)where c µ is a coefficient given in table I, and S t is theturbulence dissipation length-scale defined as S t = k / t ε t , (21)and ε t is the specific total turbulence dissipation. Theturbulence quantities k t and S t are obtained from thefollowing evolution equations, ∂k t ∂t + ˜ V j ∂k t ∂x j = P b t + P s t − k / t S t + 1 ρ ∂∂x j (cid:18) ρν t σ k ∂k t ∂x j (cid:19) , (22) ∂S t ∂t + ˜ V j ∂S t ∂x j = S t k t ( c P b t + c P s t ) − c (cid:112) k t + 1 ρ ∂∂x j (cid:18) ρν t σ S ∂S t ∂x j (cid:19) , (23) ∂a i t ∂t + ˜ V j ∂a i t ∂x j = b t ρ ∂P∂x j + R ( V i , V j ) ρ ∂ρ∂x j + ∂ ( a i t a j t ) ∂x j − c a a i t √ k t S t − a j t ∂V i ∂x j + ∂∂x j (cid:18) ν t σ a ∂a i t ∂x j (cid:19) , (24) ∂b t ∂t + ˜ V j ∂b t ∂x j = 2 a j t ∂b t ∂x j − a j t ( b t + 1) ρ ∂ρ∂x j − c b b t √ k t S t + ρ ∂∂x j (cid:18) ν t ρσ b ∂b t ∂x j (cid:19) , (25)where P b t and P s t are the specific total production ofturbulence kinetic energy by buoyancy and shear mech-anisms, P b t = a j t ρ ∂P∂x j , (26) P s t = − τ ( V i , V j ) ∂ ˜ V i ∂x j , (27)and c , c , c , c a , c b , σ k , σ S , σ a , and σ b are coefficientsof the original RANS BHR model, whose values are givenin table I.Equations 22 to 25 have been designed to operate ex-clusively with RANS variables: Reynolds averaged, Φ,density-weighted averaged, ˜Φ, and total turbulence, Φ t ,quantities. We now derive their PANS counterpart byextending the framework proposed by Girimaji [27] tovariable-density flow. To this end, the parameters defin-ing the ratios of modeled-to-total specific turbulence ki-netic energy, f k , dissipation length-scale, f S , mass fluxvelocity, f a i , and density-specific volume correlation, f b , f k ≡ k u k t , f S ≡ S u S t , f a i ≡ a i u a i t , f b ≡ b u b t , (28)need to be included in equations 22 to 25. These enablethe closure to operate at any degree of physical resolu-tion, i.e., from RANS to DNS. f S can also be calculatedas a function of f k and f ε , f S ≡ S u S t = (cid:32) k / u ε u (cid:33) (cid:32) ε t k / t (cid:33) = f / k f ε , (29)where f ε is the ratio of modeled-to-total turbulence dissi-pation. Since f ε is physically more intuitive than f S , theevolution equations for k u , S u , a i u , and b u are derived interms of f k , f ε , f a i , and f b . II.1.1. k u evolution equation It has been demonstrated by Girimaji [27] and Sumanand Girimaji [34] that the scale-invariant form of k u equa-tion can be written as, ∂k u ∂t + { V j } ∂k u ∂x j = P b u + P s u − ε u + T u , (30)TABLE I: Coefficients of BHR-LEVM closure [29]. c c c c a c b c µ σ a σ b σ c σ k σ S where P b u and P s u are the production of k u by buoyancyand shear mechanisms, P b u = a j u (cid:104) ρ (cid:105) ∂ (cid:104) P (cid:105) ∂x j , (31) P s u = − τ ( V i , V j ) ∂ { V i } ∂x j , (32)and ε u and T u represent the dissipation and transportof modeled turbulence kinetic energy. For constant f k ,differentiation commutes in time and space and so it pos-sible to establish a relationship between the equations for k u (PANS) and k t (RANS), ∂k u ∂t + ˜ V j ∂k u ∂x j = f k (cid:20) ∂k t ∂t + ˜ V j ∂k t ∂x j (cid:21) . (33)Since PANS calculates filtered or partially-averaged de-pendent variables, the former relationship can be rewrit-ten as follows, ∂k u ∂t + { V j } ∂k u ∂x j = f k (cid:20) ∂k t ∂t + ˜ V j ∂k t ∂x j (cid:21) + (cid:16) { V j } − ˜ V j (cid:17) ∂k u ∂x j . (34)Next, replacing the term between brackets by the right-hand side of equation 22 leads to ∂k u ∂t + { V j } ∂k u ∂x j = f k [ P b t + P s t − ε t + T t ]+ (cid:16) { V j } − ˜ V j (cid:17) ∂k u ∂x j , (35)and applying self-similarity considerations to the left-hand side, the following relation is obtained, P b u + P s u − ε u + T u = f k [ P b t + P s t − ε t + T t ]+ (cid:16) { V j } − ˜ V j (cid:17) ∂k u ∂x j . (36)This equation shows the formal similarity between PANS(left-hand side) and RANS (right-hand side) production,dissipation and transport terms. Hence, it is possible torelate the source and sink (local processes), and transportterms as follows, P b u + P s u − ε u = f k [ P b t + P s t − ε t ] , (37a) T u = f k T t + (cid:16) { V j } − ˜ V j (cid:17) ∂k u ∂x j . (37b) Now, we define the weighting functions ω s and ω b , ω s = P s u P b u + P s u , (38) ω b = P b u P b u + P s u , (39)which define the relative weight of the shear, ω s , andbuoyancy, ω b , mechanisms to the total production of spe-cific turbulence kinetic energy. Thus, their sum is equalto unity, ω s + ω u = 1. Using f k , f ε and these weightingfunctions, we can rewrite equation 37a, P b u + P s u − ε u ω s − ε u ω b = f k (cid:20) P b t + P s t − ω s ε u f ε − ω b ε u f ε (cid:21) , (40)and obtain the relationships, P b t = P b u f k − ω b ε u (cid:18) f k − f ε (cid:19) , (41) P s t = P s u f k − ω s ε u (cid:18) f k − f ε (cid:19) . (42)The former relations are used to derive the evolutionequation for S u . On the other hand, the transport termsof the k u and k t equations can be related as follows, T u = f k [ T t ] + (cid:16) { V j } − ˜ V j (cid:17) ∂k u ∂x j = f k (cid:20) ρ ∂∂x j (cid:18) ρν t σ k ∂k t ∂x j (cid:19)(cid:21) + (cid:16) { V j } − ˜ V j (cid:17) ∂k u ∂x j = 1 (cid:104) ρ (cid:105) ∂∂x j (cid:18) (cid:104) ρ (cid:105) ν u σ k f (cid:15) f k ∂k u ∂x j (cid:19) + (cid:16) { V j } − ˜ V j (cid:17) ∂k u ∂x j , (43)where ν u = c µ S u √ k u . Using scaling arguments, Girimaji[27] showed that (cid:16) { V j } − ˜ V j (cid:17) ∂k u ∂x j ≈ , (44)leading to the so-called zero-transport model (ZTM). Theaccuracy of this model has been confirmed in the recentstudy of Tazraei and Girimaji [46]. Also, it is importantto highlight that the velocity difference term tends tozero in the limit of f k = 0 .
00 and 1 .
00 since ∂k u ∂x j = 0 at f k = 0 . , { V j } − ˜ V j = 0 at f k = 1 . . (45)The derivation of the evolution equation for k u concludesby combining equations 35, 43, and 44, this leading to itsfinal form, ∂k u ∂t + { V j } ∂k u ∂x j = P b u + P s u − ε u + 1 (cid:104) ρ (cid:105) ∂∂x j (cid:18) (cid:104) ρ (cid:105) ν u σ k f ε f k ∂k u ∂x j (cid:19) , (46)where, ε u = k / u S u . (47)We recall that the derivation of equation 46 assumes that f k and f ε are constant. If this property does not hold, themodel’s derivation needs to consider additional terms [33]and the modeled-to-total ratio of density, f ρ . Despite be-ing commonly neglected, this requirement holds for anybridging and hybrid formulation. II.1.2. S u evolution equation The derivation of the evolution equation for S u is sim-ilar to that for k u . From f S , it is possible to establishthe following relationship between the evolution equa-tions for S u (PANS) and S t (RANS), ∂S u ∂t + ˜ V j ∂S u ∂x j = f S (cid:20) ∂S t ∂t + ˜ V j ∂S t ∂x j (cid:21) , (48)which can be rewritten as ∂S u ∂t + { V j } ∂S u ∂x j ≈ f S (cid:20) ∂S t ∂t + ˜ V j ∂S t ∂x j (cid:21) , (49)using the zero transport model [27]. Now, we replacethe material derivative of S t by the right-hand side ofequation 23, ∂S t ∂t + ˜ V j ∂S t ∂x j = f S S t k t ( c P b t + c P s t ) − f S c (cid:112) k t + f S ρ ∂∂x j (cid:18) ρν t σ S ∂S t ∂x j (cid:19) . (50)Using the parameters f k and f ε , the definition of f S , andrelationships 41 and 42, we get, ∂S u ∂t + { V j } ∂S u ∂x j = 1 (cid:104) ρ (cid:105) ∂∂x j (cid:18) (cid:104) ρ (cid:105) ν u σ S f ε f k ∂S u ∂x j (cid:19) − c f k f ε (cid:112) k u + S u k u c f k (cid:32) P b u f k − ω b k / u S u (cid:20) f k − f ε (cid:21)(cid:33) + S u k u c f k (cid:32) P s u f k − ω s k / u S u (cid:20) f k − f ε (cid:21)(cid:33) . (51) This equation can be rearranged by introducing the co-efficient c ∗ , ∂S u ∂t + { V j } ∂S u ∂x j = 1 (cid:104) ρ (cid:105) ∂∂x j (cid:18) (cid:104) ρ (cid:105) ν u σ S f ε f k ∂S u ∂x j (cid:19) − c ∗ (cid:112) k u + S u k u ( c P b u + c P s u ) , (52) c ∗ = c f k f ε + ( c ω b + c ω s ) (cid:18) − f k f ε (cid:19) . (53) II.1.3. a u i evolution equation The production terms of k u and S u in PANS BHR-LEVM closure require the calculation of the velocity massflux, a i , which is obtained from an additional evolutionequation. The derivation of the PANS equation for a i starts by establishing the following relationship, ∂a i u ∂t + ˜ V j ∂a i u ∂x j = f a i (cid:20) ∂a i t ∂t + ˜ V j ∂a i t ∂x j (cid:21) . (54)Following the approach used for k u and S u equations, theleft-hand side of equation 54 can be approximated as, ∂a i u ∂t + { V j } ∂a i u ∂x j ≈ f a i (cid:20) ∂a i t ∂t + ˜ V j ∂a i t ∂x j (cid:21) , (55)using the zero transport model [27]. Next, we replace theright-hand side of this equation by that of equation 24, ∂a i u ∂t + { V j } ∂a i u ∂x j = f a i b t ρ ∂P∂x j + f a i R ( V i , V j ) ρ ∂ρ∂x j − f a i a j t ∂V i ∂x j + f a i ∂ ( a i t a j t ) ∂x j − f a i c a a i t √ k t S t + f a i ∂∂x j (cid:18) ν t σ a ∂a i t ∂x j (cid:19) . (56)The final step to derive the a u i equation is to express theright-hand side of equation 56 in terms of filtered andunresolved quantities. This can be accomplished throughthe parameters f k , f ε , f a i , f b , and key closure simplifi-cations, f a i b t ρ ∂P∂x j ≈ f a i b u f b (cid:104) ρ (cid:105) ∂ (cid:104) P (cid:105) ∂x j , (57) f a i R ( V i , V j ) ρ ∂ρ∂x j ≈ f a i f k τ ( V i , V j ) (cid:104) ρ (cid:105) ∂ (cid:104) ρ (cid:105) ∂x j , (58) f a i a j t ∂V i ∂x j = f a i f a j a j u ∂V i ∂x j = f a i f a j a j u ∂ (cid:104) V i (cid:105) ∂x j + f a i f a j a u,j (cid:18) ∂V i ∂x j − ∂ (cid:104) V i (cid:105) ∂x j (cid:19) = f a i f a j a j u ∂ (cid:104) V i (cid:105) ∂x j (ZTM) , (59) f a i ∂ ( a i t a j t ) ∂x j = 1 f a j ∂ ( a i u a j u ) ∂x j , (60) f a i c a a i t √ k t S t = c a a i u √ k u S u f k f ε , (61) f a i ∂∂x j (cid:18) ν t σ a ∂a i t ∂x j (cid:19) = ∂∂x j (cid:18) ν u σ a f ε f k ∂a i u ∂x j (cid:19) . (62)These six terms allow us to rearrange equation 56 andobtain its final form, ∂a i u ∂t + { V j } ∂a i u ∂x j = f a i b u f b (cid:104) ρ (cid:105) ∂ (cid:104) P (cid:105) ∂x j + f a i f k τ ( V i , V j ) (cid:104) ρ (cid:105) ∂ (cid:104) ρ (cid:105) ∂x j − f a i f a j a j u ∂ (cid:104) V i (cid:105) ∂x j + 1 f a j ∂ ( a i u a j u ) ∂x j − c a a i u √ k u S u f k f ε + ∂∂x j (cid:18) ν u σ a f ε f k ∂a i u ∂x j (cid:19) . (63) II.1.4. b u evolution equation The derivation of PANS BHR-LEVM closure concludeswith the evolution equation for the unresolved density-specific volume correlation, b u . Once again, we start byestablishing the following relationship using the parame-ter f b ∂b u ∂t + ˜ V j ∂b u ∂x j = f b (cid:20) ∂b t ∂t + ˜ V j ∂b t ∂x j (cid:21) , (64)which, using the ZTM model [27], can be simplified asfollows, ∂b u ∂t + { V j } ∂b u ∂x j ≈ f b (cid:20) ∂b t ∂t + ˜ V j ∂b t ∂x j (cid:21) . (65)Replacing the term between brackets by the right-handside of equation 25, ∂b u ∂t + { V j } ∂b u ∂x j = 2 f b a j t ∂b t ∂x j − f b a j t ( b t + 1) ρ ∂ρ∂x j − f b c b b t √ k t S T + f b ρ ∂∂x j (cid:18) ν t ρσ b ∂b t ∂x j (cid:19) , (66)and converting total to partial quantities using f φ , weget the final form of the b u equation, ∂b u ∂t + { V j } ∂b u ∂x j = 2 a j u f a j ∂b u ∂x j − a j u f a j ( b u + f b ) (cid:104) ρ (cid:105) ∂ (cid:104) ρ (cid:105) ∂x j − c b b u f k f ε √ k u S u + (cid:104) ρ (cid:105) ∂∂x j (cid:18) ν u (cid:104) ρ (cid:105) σ b f ε f k ∂b u ∂x j (cid:19) . (67) Thus, the PANS BHR-LEVM closure is composed byequations 46, 52, 63, and 67. III. FILTER CONTROL PARAMETER
The efficiency of bridging and hybrid formulations isdetermined by the degree of physical resolution, i.e., therange of resolved scales or filter’s cut-off. As the modelresolves a wider range of flow scales, both the cost and ac-curacy of the simulations are expected to grow. Whereasexcessive physical resolution reduces the computationalefficiency by increasing the cost without commensurateimprovement in accuracy, insufficient resolution can com-promise accuracy by precluding the model from resolv-ing the scales not amenable to modeling [22, 42, 47, 48].Hence, the success of such SRS methods is dictated bythe parameters controlling their physical resolution.As discussed in [42], there are three main factors toconsider when determining the physical resolution neededfor a given model, flow configuration, and quantities ofinterest: i ) the length- and time-scales that need to be re-solved; ii ) the smallest flow scales that the selected spatio-temporal grid resolution and numerical setup canaccurately resolve; iii ) the effect of the physical resolution on the depen-dent variables of the turbulence closure.The authors have recently investigated the first pointthrough the analysis of flows around cylinders and theTaylor-Green vortex [22, 42, 47, 48]. These studies haveshown that the accurate prediction of these complexproblems is determined by the mathematical model’sability to resolve the instabilities and coherent structuresgoverning the flow physics. This flow physics is dom-inated by non-local effects, which most one-point clo-sures cannot represent accurately. Thus, accurate com-putations of such problems require resolving the Kelvin-Helmholtz rollers observed in flows past cylinders in thesub-critical regime [49, 50], and the vortex-reconnectionprocess of the TGV. The remaining flow scales can beaccurately modeled through an adequate turbulence clo-sure model. These studies also illustrate the importanceof understanding the flow physics and its main featuresto select the physical resolution and obtain high-fidelitysolutions.The second aspect is ideally addressed through verifi-cation exercises [38, 39, 51]. We have illustrated the cru-cial role of verification of RANS and SRS predictions in[42, 44, 48, 52, 53]. It is possible to obtain a reasonable a-priori estimate of the maximum physical resolution thata given grid can support and, conversely, the dependencyof the grid requirements from the physical resolution. Ir-respective of the numerical scheme, using Kolmogorovarguments, and assuming high-Re flow ( f ε = 1 . f k and( f k ) ref = 0 . r ∆ [42].value of f k that a spatial grid resolution can accuratelyresolve [42], f k ≥ (cid:18) c µ (cid:19) / (cid:18) ∆ k . t /ε t (cid:19) / . (68)This expression can be rearranged to provide the ratiobetween the smallest grid size for two values of f k , r ∆ = ∆ f k ∆ ( f k ) ref = (cid:18) f k ( f k ) ref (cid:19) / , (69)where ∆ is the grid resolution or size, and the subscript“ref” denotes a reference f k . This expression enables theevaluation of the relative evolution of ∆ with the physicalresolution. Figure 1 depicts the evolution of ∆ with f k relative to the case ( f k ) ref = 0 .
10. The results show theclose dependence between the grid and physical ( f k ) res-olutions. As f k increases, the minimum grid resolutioncoarsens as (10 f k ) . . Considering the cases of f k = 0 . .
40, this represents reducing the grid resolution re-quirements by a factor of four and eight when comparedto the case at f k = 0 .
10. Since SRS computations areinherently three dimensional and unsteady, figure 1 con-firms the potential of bridging methods to predict com-plex flows efficiently. It also emphasizes the importanceof selecting an adequate physical resolution for a givenproblem and quantities of interest. These estimates arebased on canonical turbulence field. The grid require-ments are expected to be less severe for under developedturbulence.The third factor is caused by the fact that an SRSmodel’s physical resolution does not affect all turbulencedependent quantities equally, i.e., a given range of re-solved scales does not lead to equal ratios modeled-to-total for all dependent variables of the closure model.For instance, it is not expected that the turbulence ki-netic energy and dissipation possess similar spectral sig-natures in a fully-developed turbulent flow [54, 55]. This has been recently addressed in Pereira et al. [42]. Weemphasize that despite most SRS models neglecting thisaspect, PANS can consider the spectral signature of eachdependent quantity of the closure through f φ . This bene-fit comes at the expense of having to determine the othercontrol parameters in such a manner that consistencybetween the various physical processes can be preserved.PANS BHR-LEVM model relies on f k , f ε , f a i , and f b to set the physical resolution. The parameter f k can beeither set constant [27, 52, 56, 57] or dynamically [56, 58–60] in space and time. Although the second approachmay enhance the simulation’s efficiency, we choose us-ing constant values of f k to prevent commutation errors[32, 33] and enable robust verification and validation ex-ercises where one can evaluate numerical and modelingerrors separately to avoid possible error canceling. Re-garding f ε , this parameter is commonly defined constantand equal to one. This modeling assumption stems fromthe fact that most turbulence dissipation in high-Re flowsoccurs at the smallest scales [54, 55]. For this reason, f ε = 1 .
00 is often used in practical PANS simulations.The validity of this option has been recently confirmedby the authors [42], and it is discussed in Section III.1.1.The remaining parameters, f a i and f b , have never beenused before and so their definition needs to be investi-gated, Section III.1.2.The remainder of this section addresses the specifica-tion of the parameters f φ in PANS BHR-LEVM, in amanner that is consistent with the implicit-filter corre-sponding to f k . Unfortunately, canonical turbulence the-ories cannot be used for this purpose. Instead, this ob-jective is accomplished through a-priori testing, in whichthe parameters f φ are calculated at successively smallerphysical resolutions - from DNS ( f k = 0 .
00) to RANS( f k = 1 . λ = 140 and 300,and the buoyancy driven homogeneous variable-densityturbulence (HVDT [62–65]) of Aslangil et al. [66, 67] atAt=0 .
75. These flow problems have been simulated bymeans of DNS, and their streamwise velocity (FHIT) anddensity (HVDT) fields are depicted in figure 2. The com-prehensive description of these data sets is given in Silvaet al. [61], Livescu and Ristorcelli [64, 65], and Aslangilet al. [66, 67].The ratios modeled-to-total, f φ , of the dependentquantities of PANS BHR-LEVM at distinct physical res-olutions are computed as follows. The DNS flow fieldsare filtered using the operator [68, 69], (cid:104) Φ (cid:105) ( x ) = (cid:90) +∆ / − ∆ / (cid:90) +∆ / − ∆ / (cid:90) +∆ / − ∆ / Φ( x ) G ∆ ( x − x (cid:48) ) d x (cid:48) , (70)where bold symbols denote vectors, ∆ is the filter’swidth, and G ∆ is the kernel of the filtering operator.Here, we use a box filter so that, G ∆ ( x − x (cid:48) ) = (cid:26) ∆ − , | x − x (cid:48) | < . , otherwise . (71) (a) FHIT - V ( x )(b) HVDT - ρ ( x , t o ) FIG. 2: FHIT and HVDT DNS velocity and density( t = 0) flow fields.Figure 3 illustrates the effect of varying this operator fil-ter’s width on the velocity field of the FHIT flow. As thefilter’s width increases, the magnitude of the filtered tur-bulent velocity field asymptotes to zero, and its gradientsget smoother. This is the reasoning for the cost reduc-tion observed from DNS to RANS [44]. It is importantto note that the shape of the filter implied by the PANSdecomposition is not known, however the utilization of abox filter is not expected to alter the conclusions of thisclass of studies [70–72] (see Section III.1.2 for HVDT).This is confirmed by comparing our HVDT results withthose of Saenz et al. [73] obtained with a Gaussian fil-ter. In contrast, cut-off filters are not suitable for thisexercise since practical SRS computations do not rely onsuch operators [70, 74].For any given filter { · } or (cid:104) · (cid:105) chosen, the unresolveddependent variables of the BHR-LEVM closure are cal-culated from relations 14 and 15 [19, 34, 75], k u = 12 (cid:88) i =1 ( { V i V i } − { V i }{ V i } ) , (72) ε u = ν (cid:18)(cid:26) ∂V i ∂x j ∂V i ∂x j (cid:27) − (cid:26) ∂V i ∂x j (cid:27) (cid:26) ∂V i ∂x j (cid:27)(cid:19) , (73) a i u = (cid:104) ρ (cid:48) v (cid:48) i (cid:105) − (cid:104) ρ (cid:48) (cid:105)(cid:104) v (cid:48) i (cid:105)(cid:104) ρ (cid:105) , (74) b u = (cid:104) ρ (cid:48) (cid:105) (cid:28) ρ (cid:48) (cid:29) − . (75)It is crucial to emphasize that ρ (cid:48) and v (cid:48) i in equations 74and 75 consider the fluctuating component of the den-sity and velocity fields, i.e., these quantities can compriseboth the coherent and stochastic fields [17, 18]. Yet, thestochastic field is expected to be the main contributorto a u i and b u at late times when the flow exhibits fully-developed and high-intensity turbulence features. Thequantities given by relations 72 to 75 are calculated with n = ∆ / ∆ η up to 349, being ∆ η the grid size used in theDNS simulations. This exercise is performed with thecode used in [44, 53].The outcome of the a-priori exercises is now discussedin Section III.1. However, before presenting the results,note that the FHIT problem is an archetypal problemwidely utilized to investigate the dynamics and modelingof fully-developed incompressible turbulence. For thisreason, we use this flow to analyze the dependence of f k and f ε on the range of resolved scales [42]. On the otherhand, the HVDT flow is a canonical problem used tostudy the fundamental physics and modeling of variable-density flow. Hence, we use the HVDT case to investigatethe evolution of f k , f a i , and f b with the physical resolu-tion. III.1.
A-priori testing results
III.1.1. Forced homogeneous isotropic turbulence
Figure 4 presents the variation of f k and f ε with therelative filter size, n = ∆ / ∆ η . n indicates how large is thefilter size when compared to DNS resolution (so n = 1 isDNS). As expected, the f k ( n ) results indicate that mostturbulence kinetic energy is contained at the largest flowscales. This behavior gets more pronounced with increas-ing Re λ . It is observed that to filter only 20% of the totalturbulence kinetic energy we need n = 29 (Re λ = 140)and 49 (Re λ = 300). This clearly illustrates the poten-tial of bridging models to efficiently compute complexflow problems. Also, note that bridging models are usu-ally not used at f k < .
20 (LES range [54]). Regarding f ε , the results of figure 4 confirm that most turbulencedissipation occurs at the smallest scales and, as such, f ε grows significantly more rapidly than f k with n . Thedata show f ε = 0 .
22 (Re λ = 140) and 0 .
34 (Re λ = 300)for n = 5, and f ε > .
99 for n = 199 and both Reynoldsnumbers.Next, figure 5 presents f ε as a function of f k . The re-sults indicate that f ε is only weakly dependent on f k at0 (a) n = 1 (b) n = 17 (c) n = 33(d) n = 66 (e) n = 199 (f) n = 349 FIG. 3: Evolution of x -component of the FHIT velocity field with the relative filter length size, n = ∆ / ∆ η .coarser f k values. At f k = 0 . f ε = 0 .
87 for Re λ = 140and f ε = 0 .
93 for Re λ = 300. Considering that practi-cal simulations of turbulence are expected to operate at f k ≥ .
20 due to the inherent cost and availability of LESformulations, figure 5 confirms that f ε = 1 .
00 is a goodassumption for practical PANS computations.
III.1.2. Homogeneous variable-density turbulence
The HVDT is a transient flow and, as such, the apriori tests are conducted at the three distinct and representa-tive time instants shown in figure 6: at t = 1 .
8, the flowis in the so-called explosive growth regime [67] and thekinetic energy of the system is rapidly increasing throughthe conversion of potential into kinetic energy. As shownin figure 6a, the flow does not exhibit small scale tur-bulence, and the two fluids ( ρ = 7 . ρ = 1 .
0) aremostly unmixed. At t = 2 .
8, the flow kinetic energygrows, reaching close to its peak. This leads to flow re-gions characterized by small scale turbulence, where thetwo fluids mix. Finally, t = 4 . f k , f a , and f b with the relative filter size, and the energy spectra of k t , a t , and b t for the unfiltered fields. Due to the HVDTflow properties, the quantities a and a are equal tozero, so we consider f a = f a . The results for f k ( n )indicate that before the peak of k t ( t = t ), most of thekinetic energy is contained in the largest coherent flowscales (blobs of laminar fluid). For this reason, f k doesnot exceed 0 .
38 when n = 99. As the flow and turbulencefield develop, t = t , k t increases to a value close to itsmaximum, altering the spectral properties of the kineticenergy field. In addition to the conversion of potentialinto kinetic energy, the energy of the largest scales istransferred to the smallest ones, widening the spectraso that larger fractions of k t are modeled for the same n . At later times, t = t , when the flow exhibits high-intensity turbulence and mixing, f k ( n ) exhibits a slightreduction. This stems from the dissipation of k t at thesmallest scales.Turning our attention to f a ( n ), the results indicatethat the evolution of this quantity with the filter widthat t = t and t is nearly identical. This shows thatthe morphological flow changes occurring at these earlyinstants significantly affect the turbulence kinetic energybut not the velocity fluctuations uncorrelated with thedensity field (figures 8b and 8c). The data of figure 7balso indicate that approximately 60% of the energy of a atthese instants is contained in the smallest flow scales, n =99. At t = t , the values of f a are reduced approximatelyfifty percent. Note that this is when the flow exhibits1 (a) f k ( n )(b) f ε ( n ) FIG. 4: Variation of f k and f ε with the relative filtersize, n , at Re λ = 140 and 300.FIG. 5: f ε as a function of f k at Re λ = 140 and 300.high-intensity turbulence, a more homogeneous mixture,and diminishing influence of the coherent field (figure 7).The variation of f b with n shows that the magnitude ofthis quantity increases from t = t to t , reaching values (a) t = t (b) t = t (c) t = t FIG. 6: Density field of HVDT flow at distinct timeinstants: before ( t ), at ( t ), and after ( t ) the peak ofkinetic energy.of 0 .
63 at t = t and 0 .
77 at t . Such result indicatethat the density-specific volume correlation is dominatedby the smallest wavelengths at these early stages. The2 (a) f k ( n )(b) f a ( n )(c) f b ( n ) FIG. 7: Variation of f k , f a , and f b with the relativefilter size, n , at distinct time instants.observed high-intensity turbulence and enhanced mixingat t = t leads to a significant reduction of f b . For n ≤ f b does not exceed 0 . (a) E ( k )(b) E ( a )(c) E ( b ) FIG. 8: Energy spectra of k t , a t , and b t at distinct timeinstants.efficiently. Considering f k , the data indicate that simu-lations at f k = 0 .
50 and 0 .
25 ( t = t ) can run on gridresolutions 89 and 33 times coarser (in each direction)than those required by DNS. This constitutes a signifi-cant cost reduction.3 (a) f a ( f k )(b) f b ( f k ) FIG. 9: Variation of f a and f b with f k at distinct timeinstants.Finally, figure 9 depicts the variation of f a and f b as afunction of f k . The results show that the ratio f a /f k getssmaller in time, and f a ≈ f k / t = t (the instant whenthe flow is characterized by high-intensity turbulence).In contrast, f b has a small temporal dependence until t ,and f b > f k . Such a behavior is not observed at t = t ,where f b ≈ f k . We attribute this result to the breakdowninto turbulence and dissipation of the coherent field. III.1.3. Parameter selection
The above a-priori tests have been conducted to de-termine the parameters f φ of PANS BHR-LEVM, andpropose guidelines toward their efficient selection. Nev-ertheless, we reiterate that the present paper’s primaryobjectives are to extend PANS methodology to variable-density flow and provide the resulting PANS BHR-LEVMmodel’s initial validation space. Closures using differentdependent variables may require similar studies to de-fine f φ (only k and ε tend to be used in most closures[14, 76, 77]). The FHIT results have shown that prescribing f ε =1 .
00 is a good strategy because most dissipation occursat the smallest flow scales. These are usually mod-eled in practical PANS computations ( f k ≥ . f k and f ε can get closer in transient and/or tran-sitional flows, the results available in the extensive PANSliterature have shown that this approach is still good[27, 42, 46, 47, 52, 56, 57, 59, 78–80]. In these cases,the resulting modeling shortcomings need to be compen-sated by slighter finer values of f k . Referring to f a and f b , selecting these parameters is more complex and hasnever been done before. Our a-priori tests suggest set-ting f a = f k / f b = f k at late times when turbulenceis closer to fully-developed, and the coherent field has adiminishing impact on the flow dynamics. However, notethat these quantities are inherently time-dependent andcontain a meaningful coherent component at early times.These aspects are particularly important for b since thisquantity is highest at t = 0. Once again, these calibra-tion issues can be overcome through a proper selectionof f φ , and slightly lower values of f k . Although oftenneglected, it is crucial to emphasize that these issues arecommon to most SRS models.Considering the previous points and the results of the a-priori tests, the present simulations utilize one of thefollowing strategies to define f k , f ε , f a ( f a = f a i ), and f b : i ) Upon inspection of the governing equations of thePANS BHR-LEVM closure, it is possible to inferthat f k and f ε are the major influence on the pro-duction of modeled turbulence by shear and buoy-ancy effects ( k u and S u equations), and, conse-quently, on the modeled turbulent stresses. Thus,we prescribe f k , and set the remaining parametersequal to one. We expect this approach to be morerobust (numerically) and general so it can be ap-plied to most flows. Yet, it might require slightlysmaller values of f k to compensate for possible cal-ibration deficits. ii ) Assume that the density fluctuations are domi-nated by the coherent field, and define f a = √ f k , f ε = 1 .
00, and f b = 1 . i ) Based on the outcome of the a-priori exercises ofthe HVDT, set f a ≈ f k / f ε = 1 .
00, and f b = f k .This strategy is optimized for HVDT type of flows,and best suited for instants characterized by fully-developed turbulence.These approaches are summarized in table II and testedin Section V.2. IV. FLOWS AND SIMULATIONS DETAILSIV.1. Taylor-Green vortex
The Taylor-Green vortex flow [36] is a canonical testcase used to investigate the modeling and simulation of4TABLE II: Modeled-to-total ratios, f φ , used in the RTPANS computations at 0 . < f k < . f k .
25 0 .
35 0 . f ε .
00 1 .
00 1 . S f a .
00 1 .
00 1 . f b .
00 1 .
00 1 . f ε .
00 1 .
00 1 . S f a .
50 0 .
59 0 . f b .
00 1 .
00 1 . f ε .
00 1 .
00 1 . S f a .
15 0 .
20 0 . f b .
25 0 .
35 0 . onset, development, and decay of turbulence [42, 81–95].The flow is initially characterized by multiple laminar,well-characterized, and single-mode vortices. These areillustrated in figure 10, and defined by [36, 82], V ( x , t o ) = V o sin( x ) cos( x ) cos( x ) , (76) V ( x , t o ) = − V o cos( x ) sin( x ) cos( x ) , (77) V ( x , t o ) = 0 , (78)where V o is the initial velocity magnitude. The corre-sponding pressure field is obtained from solving the Pois-son equation, P ( x , t o ) = P o + ρ o V o
16 [2 + cos (2 x )] [cos (2 x ) + cos (2 x )] , (79)where P o and ρ o are the pressure and density at t = 0.The vortical structures of figure 10 interact and evolvein time, and vortex stretching processes generate vortex-sheets that gradually get closer. Afterward, these vortex-sheets roll-up and reconnect [81, 96], leading to the onsetof turbulence and subsequent intensification of vorticity.The coherent structures breakdown and high-intensityturbulence appears. Finally, the turbulence kinetic en-ergy dissipates rapidly by the action of viscous effects.The analyzed TGV flow is characterized by a Reynoldsnumber Re ≡ ρL o V o /µ = 3000 [82, 84], and an initialMach number Ma o = 0 .
28. Such a Ma o leads to max-imum instantaneous and averaged ( L norm) variationsof ρ smaller than 11 .
0% and 1 .
4% of ρ o for f k = 0 .
00, re-spectively. The computational domain of this problem isa cube with length equal to L = 2 πL o . Periodic bound-ary conditions are applied on all boundaries. The ini-tial thermodynamic and flow properties are the following: V o = 10 cm/s, L o = 1 . ρ o = 1 . × − g/cm , P o = 10 Pa, µ = 3 . × − g/(cm.s), heat capacity ra-tio γ = 1 . k o = 10 − cm /s , and S o = 6 . × − cm. FIG. 10: Vortical structures present in theTaylor-Green vortex flow at t = 0. IV.2. Raleigh-Taylor flow
The RT flow [3, 4] is a benchmark problem of variable-density turbulent mixing, which has been intensely stud-ied through numerous numerical experiments [97–106].Its importance to the variable-density flow commu-nity motivated diverse validation initiatives such as theAlpha-Group collaboration [98].The flow is initially characterized by a perturbed in-terface separating two fluids of different densities, figure11. These materials are at rest, and the dense fluid, ρ h ,is on top of the light medium, ρ l . The Atwood numberof the flow is defined as At ≡ ( ρ h − ρ l ) / ( ρ h + ρ l ). Af-ter this instant, the heavy fluid starts accelerating down-wards by the action of gravity, whereas the light ma-terial moves upwards. The interface perturbations cre-ate a misalignment between the density gradient and thepressure, which induces the RT instability. The result-ing upward moving structures, named bubbles, are thepenetration of heavy fluid into the light medium and,conversely, downward moving spikes are the penetrationof light fluid into the heavy medium. The shearing mo-tion on the edges of these coherent structures triggersa Kelvin-Helmholtz instability leading to the onset anddevelopment of turbulence. As a result, the mixing rateof the two materials and the mixing-layer width increase.The temporal evolution of the RT flow comprises a linear(laminar flow) and a non-linear (laminar and turbulentflow) regimes. A comprehensive description of the flowis given by Sharp [107], Zhou [76, 77], and Boffeta andMazzino [108].The flow configuration analyzed here is based on theDNS of Livescu et al. [109] at At = 0 .
5. The com-putational domain is a rectangular prism defined in aCartesian coordinate system ( x , x , x ), figure 11. Itscross-section is L = 2 π cm wide, and the height is 3 L toensure a negligible influence of the vertical boundaries5FIG. 11: Density field of the Rayleigh-Taylor flow at t = 0.on the simulations during the simulated time T = 25time-units (time normalized by t ∗ = (cid:112) ( L/ (32 g At) [109]).The bulk Reynolds number defined as Re ≡ h ˙ h/ν canreach Re ≈
500 for the current settings ( h and ˙ h arethe mixing-layer height and its temporal derivative). Pe-riodic boundary conditions are applied on the lateralwalls and reflective conditions on the vertical boundaries, x = ± . L . The domain height and simulation timeguarantee that the simulations are not disturbed by thelatter boundary condition.The location of the interface between the two fluids isperturbed by, h p ( x , x ) = (cid:88) n,m cos (cid:104) π (cid:16) n x L + r (cid:17)(cid:105) cos (cid:104) π (cid:16) m x L + r (cid:17)(cid:105) . (80)These perturbations possess wavelengths ranging frommodes ( n and m ) 30 to 34, and amplitudes with standard-deviation not exceeding 0 . L [94, 101], figure 12. Notethat m and n are selected to include the most unstablemode of the linearized problem [109, 110]. In equation80, r and r are random numbers between 0 and 1. Thenumerical experiments rely on the ideal gas equation ofstate, and the initial temperature is set to maintain theflow Ma < .
10 and guarantee incompressible flow. Theinitial thermodynamic and flow properties are definedas follows: µ l = 0 . µ h = 0 . ρ l = 1 . , ρ h = 3 . , γ l = γ h = 1 . g = − , k o = 10 − cm / s , S o = 10 − cm, andSchmidt and Prandtl numbers are set equal to one. (a) Wave-number space(b) Physical space FIG. 12: Initial perturbations at the interface ( x = 0)in wave-number (modes) and physical space. IV.3. Numerical settings
All calculations are conducted with the flow solverxRAGE [111]. This code utilizes a finite volume ap-proach to solve the compressible and multi-material con-servation equations for mass, momentum, energy, andspecies concentration. The resulting system of governingequations is resolved [112] through the Harten-Lax-vanLeer-Contact [113] Riemann solver using a directionallyunsplit strategy, direct remap, parabolic reconstruction[114], and the low Mach number correction proposed byThornber et al. [115]. The equations are discretized withsecond-order accurate methods: the spatial discretizationis based on a Godunov scheme, whereas the temporaldiscretization relies on an explicit Runge-Kutta schemeknown as Heun’s method. The time-step, ∆ t , is definedby prescribing the maximum instantaneous CFL number,∆ t = ∆ x. CFL3( | V | + c ) , (81)where c is the speed of sound, and ∆ x is the grid cellsize. The CFL is set equal to 0 .
45 for the TGV and 0 . (1024 for f k = 0 .
00) elements for theTGV, and 128 ×
384 cells for the RT.xRage models miscible material interfaces and highconvection-driven flows with a van-Leer limiter [117],without artificial viscosity, and no material interfacetreatments [118, 119]. The solver uses the assumptionthat cells containing more than one material are in pres-sure and temperature equilibrium as a mixed cell closure.The effective kinematic viscosity in multi-material prob-lems [95] is defined as ν = n t (cid:88) n =1 ν n f n , (82)where n is the material index, n t is the number of ma-terials, and f n is the volume fraction of material n . Forthe RT flow, the diffusivity D and thermal conductiv-ity κ are defined by imposing Schmidt (Sc ≡ ν/ D ) andPrandtl (Pr ≡ c p µ/κ ) numbers equal to one.The RT computations test all three strategies for set-ting the parameters f φ given is Section III.1.3. For theTGV simulations, where there are no a i u and b u equa-tions, the three strategies are equivalent. V. RESULTS AND DISCUSSION
This section summarizes preliminary results for theTGV and RT results to illustrate the accuracy and poten-tial of the proposed variable-density PANS formulation.Additional details of the TGV results are given in Pereiraet al. [42]. More detailed analysis and finer resolutionresults for the RT will be the subject of a subsequentmanuscript.
V.1. Taylor-Green vortex
As previously mentioned, the TGV initially featuresthe laminar, single-mode, and well-defined vortical struc-tures depicted in figure 10. Immediately after t = 0,these coherent structures start interacting and deform-ing, leading to vortex-stretching processes that generatethe pairs of long sheet-like vortices observed in figure 13at t = 4 .
0. Between t = 4 . .
0, these structures getcloser and undergo a complex vortex-reconnection mech-anism [42, 82, 96] between pairs of counter-rotating vor-tices, figure 13b. This triggers the onset of the turbulenceat t ≈ .
0. Figure 13c shows the bursts of small turbu-lence scales at this instant. Afterward, turbulence fur-ther develops and eventually decays. This is illustratedin figures 13d to 13f. Considering the flow evolution, thephysics of the first nine time-units is expected to posethe greatest challenges to modeling and simulation of theTGV flow. Figure 14 presents the temporal evolution of the totalkinetic energy, k t , predicted by PANS BHR-LEVM atdifferent degrees of physical resolution, f k . Note that k t comprises a resolved, k r , and unresolved, k u , component, k t = k r + k u , (83)which are obtained from the resolved velocity field andthe turbulence closure, respectively. The results indicatethat k t is initially nearly constant, and independent of f k until t = t c ≈ .
0. At this instant, in which the flowundergoes vortex-reconnection processes, the simulationsbecome strongly dependent on f k . This result allows usto categorize the simulations into high- (HPR, f k < . f k ≥ .
50) physical resolution. The datashow that LPR simulations lead to a pronounced non-physical decay of k t . As discussed later, this is caused bya rapid increase and overprediction of the modeled tur-bulent stresses. In contrast, HPR computations exhibitsmaller energy decay rates and, as such, larger values of k t at late times. Yet, the most significant result of fig-ure 14 is that the solutions convergence upon physicalresolution refinement ( f k → t = 10.Next, figure 15 depicts the temporal evolution of thedissipation of total kinetic energy, ε t , ε t = − ∂k t ∂t , (84)and compares the results against the DNS of Brachet etal. [82] at Ma = 0 (DNS ) and Drikakis et al. [84] atMa = 0 .
28 (DNS ). The results exhibit similar tenden-cies to those of k t . Until t = t c , all PANS predictionsare independent of the physical resolution and in excel-lent agreement with the reference DNS solutions. Afterthis instant, the simulations become closely dependenton f k , but their solutions converge toward the referenceDNS data upon physical resolution refinement, f k → f k , and a good agreementwith the reference numerical experiments [82, 84]. Con-sidering the cases at f k ≤ .
25, the maximum values of ε t at 8 . ≤ t ≤ . .
145 to 0 . .
143 and 0 . k t , the differencesgrow as the physical resolution coarsens, f k → .
0. Forinstance, the magnitude of the dissipation peak can reach0 .
218 ( f k = 1 . t ≈ t = 9and 9 . (a) t = 3 . t = 5 . t = 7 . t = 9 . t = 12 . t = 20 . FIG. 13: Temporal evolution of the coherent and turbulent structures of the TGV predicted with f k = 0 .
00 [42].FIG. 14: Temporal evolution of the total kinetic energy, k t , for predictions at different f k .The results of figures 14 and 15 suggest that LPR sim-ulations prematurely predict the onset of turbulence, andoverpredict the turbulence. This would explain the rapiddecay of k t observed in figure 14. These ideas are sup-ported by the ratio of unresolved-to-total kinetic energy, k u /k t , depicted in figure 16. The data indicate that k u /k t FIG. 15: Temporal evolution of the total kinetic energydissipation, ε t , for predictions at different f k .is negligible and independent of f k until t = t c . Afterthis time, k u /k t grows considerably, and its magnitudebecomes closely dependent on f k . Also, it is visible thatthe growth of k u /k t starts earlier and it is more rapidfor simulations at f k = 1 .
00 than at f k = 0 .
25. For ex-ample, k u /k t predicted at t = 20 and f k = 1 .
00 is six8FIG. 16: Temporal evolution of the ratiomodeled-to-total kinetic energy, k u /k t , for predictionsat different f k .times larger than that at f k = 0 .
25. Considering theresults for ε t , this shows that simulations at f k = 1 . k u /k t only starts growing rapidlyat t ≥ . t = 6 . f k . The plots show that the LPRsimulation ( f k = 1 .
00) dissipates the coherent structuresinvolved in the vortex-reconnection processes. This iscaused by the overprediction of turbulence since ν t (seeequation 17) can exceed its laminar counterpart by a fac-tor of 30. This does not occur at f k = 0 .
35 (highest HPR f k ) nor at f k = 0 .
00. A comprehensive assessment of thisflow is given in [42].In summary, the results indicate that the proposedPANS BHR-LEVM model can accurately predict theshear driven TGV flow using values of f k smaller than0 .
50. Next, we evaluate the performance of the modelpredicting the buoyancy driven RT flow.
V.2. Rayleigh-Taylor
The present RT flow is initialized with the perturbedinterface shown in figures 11 and 12. Immediately afterthis instant, the two fluids accelerate, and the interfaceperturbations create a misalignment between the densitygradient and the pressure. This leads to the generationof coherent structures called spikes and bubbles [76, 77]with the mushroom-like shape illustrated in figures 18a (a) f k = 0 . f k = 0 . f k = 1 . FIG. 17: Temporal evolution of the coherent andturbulent structures of the TGV at t = 6 . f k [42].and 19. During this period, the flow is laminar, and itis in the so-called linear regime [76, 77]. In the follow-ing instants, figures 18b to 18d, the mixing-layer contin-ues growing, and the initially linear structure and theKelvin-Helmholtz secondary instability will eventuallytrigger the onset and development of turbulence. Thisphenomenon increases the mixing rate and enhances themixture homogeneity, occurring in the non-linear regime[76, 77]. It is particularly pronounced at t = 20 . (a) t = 2 . t = 5 . t = 10 . t = 20 . FIG. 18: Temporal evolution of the RT density field predicted at f k = 0 . t = 2 . f k = 0 . h , and density field, χ , predicted atdifferent f k . Also, the simulations test the three strate-gies proposed in Section III.1.3 to define the relationshipbetween the parameters f φ ( f k , f ε , f a , and f b ). Here, thequantity χ used to analyze the density field is defined as χ = ρ − ρ l ρ h − ρ l . (85)The mixing layer width h is defined as the distance be-tween the locations ¯ χ = 0 .
05 and 0 .
95, where here thebar indicates a planar average normal to gravity. Both h and x are normalized by L . Since the physics and statis-tics of RT flow are highly dependent on initial conditions and settings, we use the solutions obtained at f k = 0 . f k , and converge upon this pa-rameter’s refinement, f k → .
00. Comparing the solutionat f k = 0 .
00 against that at f k = 1 .
00 shows that the firstcase leads to a significantly thicker mixing-layer, charac-terized by a longer linear region ( t ≈ f k = 1 .
00 (RANS)prematurely predicts the onset of turbulence by overpre-dicting the total turbulent stresses. The refinement of f k improves the simulations by reducing the discrepanciesagainst the reference solution ( f k = 0 . f k ≤ .
25 are in excel-lent agreement.Regarding the strategies to define the relationship be-tween the parameters f φ ( f k , f ε , f a , f b ), the data indicatethat the approach S permits using larger values of f k forthe same degree of accuracy. For instance, it is observedthat simulations using S lead to reduced discrepanciesbetween solutions at f k ≤ .
50, whereas the remainingstrategies require f k ≤ .
25. This indicates that strategy S can run on lower physical resolutions (larger f k ) than S and S . Note that figure 20c shows some fluctuationsin h using S at late times. We attribute this result tonumerical uncertainty [52].Figure 21 shows how the density field evolves in timeand space with f k . For conciseness, we only show threerepresentative cases: f k = 0 . f k = 0 .
25 ( S ), and f k = 1 .
00. The remaining cases are in line with theresults of figure 20. As for the quantity h , the resultsshow distributions of χ quite similar between solutionsobtained at f k ≤ .
25. In contrast, those obtained fromsimulations at f k = 1 .
00 exhibit steeper and smoother0 (a) S .(b) S .(c) S . FIG. 20: Temporal evolution of the mixing-layer height, h , predicted at different f k . (a) f k = 0 . f k = 0 .
25 using S .(c) f k = 1 . FIG. 21: Temporal evolution of the density field, χ ,predicted at different f k .1profiles. This behavior stems from the caveats of RANS( f k = 1 .
00) closures predicting transient turbulence andthe fact that this modeling strategy does not resolve tur-bulence (lower numerical requirements).Next, figure 22 presents the evolution of the maximumplanar ( x − x ) averaged value of ν u obtained at differentphysical resolutions. Note that ν u is here utilized to eval-uate how the unresolved turbulent stresses evolve with f k (see equation 17). Also, it important to emphasize thatratios ν t /ν exceeding O (1) are usually attributed to tur-bulent flow. As expected, the results show that ν u /ν decreases with f k . Considering t = 25, ν u /ν reducesfrom 6128 . f k = 1 .
00 to 479 . S ), 122 . S ), and63 . S ) at f k = 0 .
25. This also shows that the strategyused to define the relationship between the different f φ has an important impact on the magnitude of ν u and,consequently, τ ( V i , V j ).However, the most significant result in figure 22 is thefact that ν u /ν does not exceed 3 . f k = 0 .
25 and t ≤ . f k = 1 .
00. Such a result indicates thatRANS ( f k = 1 .
00) predicts turbulent flow at this early in-stant, this being the consequence of the premature onsetof turbulence. On the other hand, PANS at f k = 0 .
25 cancapture the coherent structures that later are involved onthe onset and development of turbulence. This outcomeexplains the results of figures 20 and 21 and can be seenin figure 23, in which the flow coherent structures aredepicted for f k = 0 .
25 ( S ) and 1 .
00 at t = 2 .
5. Com-pared to figure 19, the results show that simulations at f k = 1 .
00 dissipate all the coherent structures. This iscaused by the overprediction of turbulence, i.e., the mag-nitude of ν t or τ ( V i , V j ). In clear contrast, simulationsat f k = 0 .
25 can accurately predict these coherent struc-tures, this being the reason for the excellent agreementbetween simulations at f k ≤ .
25 (for all S i ). Hence,the RT computations reinforce the importance of resolv-ing the instabilities and coherent structures governing theflow dynamics to obtain efficient high-fidelity simulations[22].Overall, the results have shown the ability of the PANSBHR-LEVM model to predict the current RT flow ac-curately. Regarding the strategies to prescribe f φ , S leads to the lowest values of ν u , allowing the utilizationof larger f k than S and S . A comprehensive analysisof this problem will be given in a subsequent paper. VI. CONCLUSIONS
We extended the framework of the PANS model tovariable-density flow, i.e., multi-material and/or com-pressible mixing problems including density fluctuationsand production of turbulence kinetic energy by shear andbuoyancy mechanisms. The framework was then utilizedto derive the PANS version of the k − S − a i − b equa-tion BHR-LEVM closure. The parameters defining thephysical resolution of the model ( f k , f ε , f a , and f b ) havebeen studied through a-priori testing. Three strategies (a) S .(b) S .(c) S . FIG. 22: Temporal evolution of the ratio ν u /ν predictedat different f k .2 (a) f k = 0 .
25 and S .(b) f k = 1 . FIG. 23: RT structures at t = 2 . f k = 0 .
25 and 1 .
00. are proposed to set these parameters as a function of f k : i) define f ε = f a i = f b = 1 . ii) prescribe f a i = 0 . f k and f ε = f b = 1 .
0; and iii) set f ε = 1 . f a = 0 . f k ,and f b = f k . All strategies lead to high-fidelity simula-tions. However, the first two may reveal more robust atearly times because the third approach considers devel-oped turbulence. Future studies will further investigatethis aspect. The initial validation space of the PANSBHR-LEVM model comprises the TGV at Re = 3000and the RT at At = 0 . max ≈
500 flows. Theresults confirm the ability of the model to calculate theserepresentative flows accurately. Hence, this initial valida-tion space and the theoretical justification demonstratethe new methodology’s potential to predict complex flowsof variable-density flow. Nevertheless, subsequent stud-ies will be needed to study and extend the validationspace of the model further. Finally, all simulations in-dicate the importance of resolving the instabilities andcoherent structures governing the flow dynamics to ob-tain high-fidelity simulations.
ACKNOWLEDGMENTS
We would like to thank C .B. da Silva, D. Aslangil, andD. Livescu for sharing their DNS data sets. Los AlamosNational Laboratory (LANL) is operated by TRIAD Na-tional Security, LLC for the US DOE NNSA. This re-search was funded by LANL Mix and Burn project underthe DOE ASC, Physics and Engineering Models program. [1] H. von Helmholtz, On Discontinuous Movements of Flu-ids, The London, Edinburgh, and Dublin PhilosophicalMagazine and Journal of Science , 337 (1868).[2] W. Thomson, Hydrokinetic Solutions and Observa-tion, The London, Edinburgh, and Dublin PhilosophicalMagazine and Journal of Science , 362 (1871).[3] L. Rayleight, Investigation of the Character of the Equi-librium of an Incompressible Heavy Fluid of VariableDensity, Proceedings of the London Mathematical Soci-ety , 170 (1882).[4] G. Taylor, The Instability of Liquid Surfaces when Ac-celerated in a Direction Perpendicular to their Planes,Proceeding of the Royal Society A , 192 (1950).[5] R. Richtmyer, Taylor Instability in Shock Accelerationof Compressible Fluids, Communications on Pure andApplied Mathematics , 297 (1960).[6] E. Meshkov, Instability of the Interface of Two gasesAccelerated by a Shock Wave, Soviet Fluid Dynamics , 151 (1969).[7] J. Smagorinsky, General Circulation Experiments withthe Primitive Equations I. The Basic Experiment,Monthly Weather Review , 99 (1963).[8] O. Reynolds, On the Dynamical Theory of Incompress-ible Viscous Fluids and the Determination of the Crite-rion, Philosophical Transactions of the Royal Society ofLondon , 123 (1985). [9] A. Favre, ´Equations Statistiques des Gaz Turbulens,Comptes Rendus de l’Acad´emie des Sciences Paris (1958).[10] A. Favre, ´Equations Statistiques des Gaz Turbu-lens Compressibles. I. Formes G´en´erales, Journal deM´ecanique , 361 (1965).[11] A. Favre, ´Equations Statistiques des Gaz TurbulensCompressibles. II. M´ethode des Vitesses Moyennes;M´ethod des Vitesses Macroscopiques Pond´er´ees par laMasse Volumique, Journal de M´ecanique , 391 (1965).[12] A. Favre, ´Equations Statistiques aux Fluctuationsd´Entropie, de Concentration, de Rotationnel dansles ´Ecoulements Compressibles, Comptes Rendus del’Acad´emie des Sciences Paris , 1289 (1971).[13] D. Besnard, F. H. Harlow, R. M. Rauenzahn,and C. Zemach, Turbulence Transport Equations forVariable-Density Turbulence and Their Relationship toTwo-Field Models , Tech. Rep. LA-12303-MS, DE92017292 (Los Alamos National Laboratory, 1992).[14] P. Chassaing, The Modeling of Variable Density Turbu-lent Flow, Flow, Turbulence and Combustion , 293(2001).[15] P. Chassaing, R. Antonia, F. Anselmet, L. Joly, andS. Sarkar, Variable Density Fluid Turbulence , editedby R. Moreau, Fluid Mechanics and Its Applications,Vol. 69 (Springer-Science+Business Media, B.V., 2002). [16] D. Wilcox, Turbulence Modeling for CFD , 3rd ed.(DCW Industries, La Ca˜nada, United States of Amer-ica, 2006).[17] A. K. M. F. Hussain and W. C. Reynolds, The Me-chanics of an Organized Wave in Turbulent Shear Flow,Journal of Fluid Mechanics , 241 (1970).[18] R. Schiestel, Multiple-Time-Scale Modeling of Turbu-lent Flows in One-Point Closures, Physics of Fluids ,722 (1987).[19] M. Germano, Turbulence: the Filtering Approach, Jour-nal of Fluid Mechanics , 325 (1992).[20] M. Germano, From RANS to DNS: Towards a BridgingModel, in Direct and Large-Eddy Simulation III, ER-COFTAC series (Springer Netherlands, 1999).[21] C. Speziale, Computing Non-Equilibrium TurbulentFlows With Time-Dependent RANS and VLES, in th International Conference on Numerical Methods inFluid Dynamics , Lecture Notes in Physics, Vol. 490,edited by P. Kutler and J. F. J. Chattot (Springer,Berlin, Germany, 1997).[22] F. S. Pereira, L. E¸ca, G. Vaz, and S. S. Girimaji,Challenges in Scale-Resolving Simulations of TurbulentWake Flows with Coherent Structures, Journal of Com-putational Physics , 98 (2018).[23] P. Batten, U. Goldberg, and S. Chakravarthy, Sub-GridTurbulence Modeling for Unsteady Flow With Acoun-stic Resonance, in th American Institute of Aeronau-tics and Astronautics Aersopaces Science Meeting andExhibit , AIAA-00-0473 (Reno, United States of Amer-ica, 2000).[24] H. Fasel, J. Seidel, and S. Wernz, A Methodology forSimulations of Complex Turbulent Flows, Journal ofFluids Engineering , 993 (2002).[25] R. Schiestel and A. Dejoan, Towards a New PartiallyIntegrated Transport Model for Coarse Grid and Un-steady Turbulent Flow Simulations, Theoretical andComputational Fluid Dynamics , 443 (2005).[26] B. Chaouat and R. Schiestel, A New Partially In-tegrated Transport Model for Subgrid-Scale Stressesand Dissipation Rate for Turbulent Developing Flows,Physics of Fluids (2005).[27] S. S. Girimaji, Partially-Averaged Navier-Stokes Modelfor Turbulence: A Reynolds-Averaged Navier-Stokes toDirect Numerical Simulation Bridging Method, Journalof Applied Mechanics , 413–421 (2005).[28] F. S. Pereira, L. E¸ca, and G. Vaz, Investigating the Ef-fect of the Closure in Partially-Averaged Navier-StokesEquations, Journal of Fluids Engineering , 121402,FE (2019).[29] K. Stalsberg-Zarling and R. Gore, The BHR2 Turbu-lence Model: Incompressible Isotropic Decay, Rayleigh-Taylor, Kelvin-Helmholtyz and Homogenous VariableDensity Turbulence , Tech. Rep. LA-UR-11-04773 (LosAlamos National Laboratory, 2011).[30] J. D. Schwarzkopf, D. Livescu, J. R. Baltzer, R. A. Gore,and J. R. Ristorcelli, A Two-Length Scale TurbulenceModel for Single-Phase Multi-Fluid Mixing, Flow, Tur-bulence and Combustion , 1 (2016).[31] S. S. Girimaji, E. Jeong, and R. Srinivasan, PartiallyAveraged Navier-Stokes Method for Turbulence: FixedPoint Analysis and Comparison with Unsteady PartiallyAveraged Navier-Stokes, Journal of Applied Mechanics , 422 (2006).[32] F. Hamba, Analysis of Filtered Navier-Stokes Equationfor Hybrid RANS/LES Simulation, Physics of Fluids (2011).[33] S. S. Girimaji and S. Wallin, Closure Modeling in Bridg-ing Regions of Variable-Resolution (VR) TurbulenceComputations, Journal of Turbulence , 72 (2013).[34] S. Suman and S. S. Girimaji, On the Invariance of Com-pressible Navier-Stokes and Energy Equations Subjectto Density Weighted Filtering, Flow, Turbulence andCombustion , 383 (2010).[35] A. Banerjee, R. A. Gore, and M. J. Andrews, Devel-opment and Validation of a Turbulent-Mix Model forVariable-Density and Compressible Flows, Physical Re-view E , 046309 (2010).[36] G. Taylor and A. Green, ”mechanism of the productionof small eddies from large ones”, in Proceedings of theRoyal Society A , Vol. 158 (1937) pp. 499–521.[37] L. Rayleigh, Investigation of the Character of the Equi-librium of an Incompressible Heavy Fluid of VariableDensity, Proceedings of the London Mathematical Soci-ety , 170 (1982).[38] T. Trucano, M. Pilch, and W. Oberkampf, General Con-cepts for Experimental Validation of ASCI Code Ap-plications , Technical Report SAND2002-0341 (SandiaNational Laboratories, Albuquerque, United States ofAmerica, 2002).[39] W. L. Oberkampf and C. J. Roy,
Verification and Vali-dation in Scientific Computing , 1st ed. (Cambridge Uni-versity Press, Cambridge, United Kingdom, 2010).[40] F. A. Williams,
Combustion Theory: The Fundamen-tal Theory of Chemically Reacting Flow Systems , 1sted. (Addison-Wesley Publishing Company, Inc, Read-ing, Massachussetts, USA, 1965).[41] A. W. Cook, Enthalpy Diffusion in MulticomponentFlows, Physics of Fluids , 055109 (2009).[42] F. S. Pereira, F. F. Grinstein, D. M. Israel,R. Rauenzahn, and S. S. Girimaji, Modelingand Simulations of Transitional Taylor-GreenVortex Flow with PANS Method, (submitted)(https://arxiv.org/abs/2012.01304).[43] J. Boussinesq, Essai sur la Th´eorie des Eaux Courantes,M´emoires Pr´esent´es par Divers Savants `a L’Acad ´mie desSciences , 335 (1877).[44] F. S. Pereira, L. E¸ca, G. Vaz, and S. S. Girimaji, To-ward Predictive RANS and SRS Computations of Tur-bulentExternal Flows of Practical Interest, Archives ofComputational Mechanics and Engineering (2021).[45] J. D. Schwarzkopf, D. Livescu, R. A. Gore, and R. M.R. J. R. Ristorcelli, Application of a Second-MomentClosure Model to Mixing Processes Involving Multi-component Miscible Fluids, Journal of Turbulence ,1 (2011).[46] P. Tazraei and S. S. Girimaji, Scale-Resolving Simula-tions of Turbulence: Equilibrium Boundary Layer Anal-ysis Leading to Near-Wall Closure Modeling, PhysicalReview Fluids , 104607 (2019).[47] F. S. Pereira, L. E¸ca, G. Vaz, and S. S. Girimaji, Onthe Simulation of the Flow Around a Circular Cylinderat Re=140,000, International Journal of Heat and FluidFlow , 40 (2019).[48] F. S. Pereira, G. Vaz, and L. E¸ca, Evaluation of RANSand SRS Methods for Simulation of the Flow Arounda Circular Cylinder in the Sub-Critical Regime, OceanEngineering (2019).[49] C. H. K. Williamson, Vortex Dynamics in the Cylin-der Wake, Annual Review of Fluid Mechanics , 477(1996). [50] M. M. Zdravkovich, Flow Around Circular Cylinders.Volume 1: Fundamentals , 1st ed. (Oxford Science Pub-lications, Oxford, United Kingdom, 1997).[51] P. J. Roache,
Verification and Validation in Computa-tional Science and Engineering , 1st ed. (Hermosa Pub-lishers, Albuquerque, United States of America, 1998).[52] F. S. Pereira, G. Vaz, L. E¸ca, and S. S. Girimaji, Sim-ulation of the Flow Around a Circular Cylinder at Re=3900 with Partially-Averaged Navier-Stokes Equations,International Journal of Heat and Fluid Flow , 234(2018).[53] F. S. Pereira, Towards Predictive Scale-Resolving Simu-lations of Turbulent External Flows , phdthesis, InstitutoSuperior T´ecnico, Lisbon, Portugal (2018).[54] S. B. Pope,
Turbulent Flows , 1st ed. (Cambridge Uni-versity Press, Cambridge, United Kingdom, 2000).[55] P. A. Davidson,
Turbulence: An Introduction for Scien-tists and Engineers (Oxford University Press, 2006).[56] S. S. Girimaji and K. S. Abdol-Hamid, Partially-Averaged Navier Stokes Model for Turbulence: Imple-mentation and Validation, in
Proceedings of the 43 rd American Institute of Aeronautics and Astronautics(AIAA) Aerospace Sciences Meeting and Exhibit (AIAA2005-502, Reno, United States of America, 2005).[57] S. Lakshmipathy and S. S. Girimaji, Partially AveragedNavier–Stokes (PANS) Method for Turbulence Simula-tions: Flow Past a Circular Cylinder, Journal of FluidsEngineering (2010), 121202.[58] A. Elmiligui, K. S. A.-H. S. J. Massey, and S. P. Pao,Numerical Study of Flow Past a Circular Cylinder Us-ing RANS, Hybrid RANS/LES and PANS Formula-tions, in nd Applied Aerodynamics Conference and Ex-hibit , AIAA 2004-4959 (Providence, Rhode Island, USA,2004).[59] B. Basara and Z. P. S. Girimaji, A New Approach forthe Calculation of the Cut-Off Resolution Parameter inBridging Methods for Turbulent Flow Simulation, Inter-national Journal of Heat and Fluid Flow , 76 (2018).[60] L. Davidson and C. Friess, A New Formulation of f k forthe PANS model, Journal of Turbulence , 322 (2019).[61] T. S. Silva, M. Zecchetto, and C. B. da Silva, The Scal-ing of the Turbulent/Non-Turbulent Interface at HighReynolds Numbers, Journal of Fluid Mechanics (2018).[62] G. Batchelor, V. Canuto, and J. Chasnov, HomogeneousBuoyancy-Generated Turbulence, Journal of Fluid Me-chanics , 349 (1992).[63] D. L. Sandoval, T. T. Clark, and J. J. Riley, Buoyancy-Generated Variable-Density Turbulence, in IUTAMSymposium on Variable Density Low-Speed TurbulentFlows , Fluid Mechanics and Its Applications, Vol. 41,edited by L. Fulachier, J. L. Lumley, and F. Ansel-met (Springer Netherlands, Dordrecht, The Nether-lands, 1997) pp. 173–180.[64] D. Livescu and J. R. Ristorcelli, Buyancy-DrivenVariable-Density Turbulence, Journal of Fluid Mechan-ics , 43 (2007).[65] D. Livescu and J. R. Ristorcelli, Variable-Density Mix-ing in Buyancy-Driven Turbulence, Journal of Fluid Me-chanics , 145 (2008).[66] D. Aslangil, D. Livescu, and A. Banerjee, Variable-Density Buoyancy-Driven Turbulence with AsymmetricInitial Density Distribution, Physica D: Nonlinear Phe-nomena , 132444 (2020). [67] D. Aslangil, D. Livescu, and A. Banerjee, Effects ofAtwood and Reynolds Numbers on the Evolution ofBuoyancy-Driven Homogeneous Variable-Density Tur-bulence, Journal of Fluid Mechanics (2020).[68] C. da Silva,
The Role of Coherent Structures in theControl and Insterscale Interactions of Round, Planeand Coaxial Jets , Phd thesis, Institut National Poly-technique de Grenoble, Grenoble, France (2001).[69] C. B. da Silva, S. Rego, and C. F. Pereira, Analysis ofthe Viscous/Molecular Subgrid-Scale Dissipation Termsin LES Based on Transport Equations:
A Priori
Tests,Journal of Turbulence , N25 (2008).[70] C. B. da Silva and J. C. F. Pereira, Analysis of theGradient-Diffusion Hypothesis in Large-Eddy Simula-tions Based on Transport Equations, Physics of Fluids (2007).[71] V. Borue and S. A. Orszag, Local Energy Flux andSubgrid-Scale Statistics in Three-Dimensional Turbu-lence, Journal of Fluid Mechanics , 1 (1998).[72] S. Liu, C. Meneveau, and J. Katz, On the Properties ofSimilarity Subgrid-Scale Models as Deduced from Mea-surements in a Turbulent Jet, Journal of Fluid Mechan-ics , 83 (1994).[73] J. A. Saenz, D. Aslangil, and D. Livescu, Filtering, Av-eraging and Scale Dependency in Homogeneous VariableDensity Turbulence, Physics of Fluids (2021).[74] B. Vreman, B. Geurts, and H. Kuerten, RealizabilityConditions for the Turbulent Stress Tensor in Large-Eddy Simulation, Journal of Fluid Mechanics , 351(1994).[75] P. Moin, K. Squires, W. Cabot, and S. Lee, A dynamiocSubgrid-Scale Model for Compressible Turbulence andScalar Transport, Physics of Fluids A: Fluid Dynamics , 2746 (1991).[76] T. Zhou, Rayleigh-Taylor and Richtmyer-Meshkov In-stability Induced Flow, Turbulence, and Mixing. I,Physics Reports , 1 (2017).[77] T. Zhou, Rayleigh-Taylor and Richtmyer-Meshkov In-stability Induced Flow, Turbulence, and Mixing. II,Physics Reports , 1 (2017).[78] F. S. Pereira, G. Vaz, and L. E¸ca, An Assessment ofScale-Resolving Simulation Models for the Flow Arounda Circular Cylinder, in Proceedings of the 8 th Inter-national Symposium on Turbulence, Heat and MassTransfer (THMT15) (Sarajevo, Bosnia and Herzegov-ina, 2015).[79] C. Kamble and S. Girimaji, Characterization of Coher-ent Structures in Turbulent Wake of a Sphere UsingPartially Averaged Navier–Stokes (PANS) Simulations,Physics of Fluids (2020).[80] T. S. Fowler, F. D. Witherden, and S. S. Girimaji,Partially-Averaged Navier–Stokes Simulations of Tur-bulent Flow Past a Square Cylinder: Comparative As-sessment of Statistics and Coherent Structures at Dif-ferent Resolutions, Physics of Fluids , 125106 (2020).[81] Y. Yang and D. I. Pullin, Evolution of Vortex-SurfaceFields in Viscous Taylor-Green and Kida-Pelz Flows,Journal of Fluid Mechanics , 146 (2011).[82] N. Brachet, D. Meiron, S. Orszag, B. Nickel, R. Morf,and U. Frish, Small Scale Structure of the Taylor-GreenVortex, Journal of Fluid Mechanics , 411 (1983).[83] C.-W. Shu, W.-S. Don, D. Gottlieb, O. Schilling, andL. Jameson, Numerical Convergence Study of Nearly In-compressible, Inviscid Taylor–Green Vortex Flow, Jour-nal of Scientific Computing , 1 (2005). [84] D. Drikakis, C. Fureby, F. F. Grinstein, and D. Youngs,Simulation of Transition and Turbulence Decay in theTaylor-Green Vortex, Journal of Turbulence (2007).[85] Y. Yang and D. I. Pullin, On Lagrangian and Vortex-Surface Fields for Flows with Taylor-Green and Koda-Pelz Initial Conditions, Journal of Fluid Mechanics ,446 (2010).[86] J.-B. Chapelier, M. Plata, and F. Renac, Inviscidand Viscous Simulations of the Taylor-Green VortexUsing a Model Discontinuous Galerkin Approach, in nd American institute of Aeronautics and Astronautics(AIAA) Fluid Dynamics Conference and Exhibit , AIAA2012-3073 (New Orleans, USA, 2012).[87] J. R. DeBonis,
Solutions of the Taylor-Green VortexProblem Using High-Resolution Explicit Finite Differ-ence Methods , techreport NASA/TM-2013-217850 (Na-tional Aeronautics and Space Administration (NASA),Glenn Research Center, Cleveland, USA, 2013).[88] I. A. Shirokov and T. G. Elizarova, Simulation ofLaminar-Turbulent Transition in Compressible Taylor-Green Flow Basing on Quasi-Gas Dynamic Equation,Journal of Turbulence , 707 (2014).[89] J. R. Bull and A. Jameson, Simulation of the Taylor-Green Vortex using high-Order Flux ReconstructionSchemes, American Institute of Aeronautics and Astro-nautics (AIAA) Journal , 2750 (2015).[90] T. Dairay, E. Lamballais, S. Laizet, and J. C. Vassili-cos, Numerical Dissipation vs. Subgrid-Scale Modellingfor Large Eddy Simulation, Journal of ComputationalPhysics , 252 (2017).[91] R. E. Moura, G. Mengaldo, and J. P. amd S. J. Sher-win, On the Eddy-Resolving Capability of High-Order Discontinuous Galerkin Approaches to ImplicitLES/Under-Resolved DNS of Euler Turbulence, Jour-nal of Computational Physics , 615 (2017).[92] N. Peng and Y. Yang, Effects of the Mach Number onthe Evolution of Vortex-Surface Fields in CompressibleTaylor-Green Flows, Physical Review Fluids (2018).[93] N. Sharma and T. K. Sengupta, Vorticity Dynamics ofthe Three-Dimensional Taylor-Green Vortex Problem,Physics of Fluids (2019).[94] F. F. Grinstein, J. A. Saenz, J. C. Dolence, T. O.Masser, R. M. Rauenzahn, and M. M. Francois, Ef-fects of Operator Splitting and Low Mach-Number Cor-rection in Turbulent Mixing Transition Simulations,Computers and Mathematics with Applications , 437(2019).[95] F. S. Pereira, F. F. Grinstein, D. M. Israel, andR. Rauenzahn, Molecular Viscosity and Diffusivity Ef-fects in Transitional and Shock-Driven Mixing Flows,Physical Review E , 013106 (2021).[96] S. Kida and M. Takaoka, Vortex Reconnection, AnnualReviews of Fluid Mechanics , 169 (1994).[97] A. W. Cook and P. E. Dimotakis, Transition Stagesof Rayleigh–Taylor Instability Between Miscible Fluids,Journal of Fluid Mechanics , 69 (2001).[98] G. Dimonte, D. L. Youngs, A. Dimits, S. Weber,M. Marinak, S. Wunsch, C. Garasi., A. Robinson, M. J.Andrews, P. Ramaprabhu, A. C. Calder, B. Fryxell,J. Biello, L. Dursi, P. MacNeice, K. Olson, P. Ricker,R. Rosner, F. Timmes, H. Tufo, Y.-N. Young, andM. Zingale, A Comparative Study of the Turbu-lent Rayleigh–Taylor Instability using High-ResolutionThree-Dimensional Numerical Simulations: The Alpha-Group Collaboration, Physics of Fluids , 1668 (2004). [99] J. R. Ristorcelli and T. T. Clark, Rayleigh–TaylorTurbulence: Self-Similar Analysis and Direct Numeri-cal Simulations, Journal of Fluid Mechanics , 213(2004).[100] W. Cabot and A. Cook, Reynolds Number Effects onRayleigh–Taylor Instability with Possible Implicationsfor Type Ia Supernovae, Nature Physics , 562 (2006).[101] A. Banerjee and M. J. Andrews, 3D Simulations to In-vestigate Initial Condition Effects on the Growth ofRayleigh–Taylor Mixing, International Journal of Heratand Mass Transfer , 3906 (2009).[102] D. Livescu, J. R. Ristorcelli, R. A. Gore, S. H. Dean,W. H. Cabot, and A. W. Cook, High-Reynolds NumberRayleigh–Taylor Turbulence, Journal of Turbulence ,N13 (2009).[103] N. Vladimirova and M. Chertkov, Self-Similarity andUniversality in Rayleigh–Taylor, Boussinesq Turbu-lence, Physics of Fluids , 015102 (2009).[104] D. Livescu, Numerical Simulations of Two-Fluid Turbu-lent Mixing at Large Density Ratios and Applicationsto the Rayleigh-Taylor Instability, Philosophical Trans-actions of the Royal Society A: Mathematical, Physicaland Engineering Sciences , 20120185 (2013).[105] D. L. Youngs, Rayleigh–Taylor Mixing: Direct Numer-ical Simulation and Implicit Large Eddy Simulation,Physica Scripta , 074006 (2017).[106] I. W. Konnikakis, D. Drikakis, and D. L. Youngs, Model-ing of Rayleigh-Taylor Mixing Using Single-Fluid Mod-els, Physical Review E , 013104 (2019).[107] D. H. Sharp, An Overview of Rayleigh-Taylor Instabil-ity, Physica D: Nonlinear Phenomena , 3 (1984).[108] G. Boffeta and A. Mazzino, IncompressibleRayleigh–Taylor Turbulence, Annual Review ofFluid Mechanics , 119 (2017).[109] D. Livescu, T. Wei, and P. T. Brady, Rayleigh–TaylorInstability with Gravity Reversal, Physica D: NonlinearPhenomena (2021).[110] R. E. Duff, F. H. Harlow, and C. W. Hirt, Effects of Dif-fusion on Interface Instability Between Gases, Physicsof Fluids , 417 (1962).[111] M. Gittings, R. Weaver, M. Clover, T. Betlach,N. Byrne, R. Coker, E. Dendy, R. Hueckstaedt, K. New,W. R. Oakes, D. Ranta, and R. Stefan, The RAGERadiation-Hydrodynamic Code, Computational Science& Discovery (2008).[112] F. S. Pereira, F. F. Grinstein, and D. Israel, Effect of theNumerical Discretization Scheme in Shock-Driven Tur-bulent Mixing Simulations, Computers & Fluids ,104487 (2020).[113] E. F. Toro, M. Spruce, and W. Speares, Restoration ofthe Contact Surface in the HLL-Riemann Solver, ShockWaves , 25 (1994).[114] P. Colella and P. R. Woodward, The PiecewiseParabolic Method (PPM) for Gas-Dynamical Simula-tions, Journal of Computational Physics , 174 (1984).[115] B. Thornber, A. Mosedale, D. Drikakis, D. Youngs,and R. J. R. Williams, An Improved ReconstructionMethod for Compressible Flows with Low Mach Num-ber Features, Journal of Computational Physics ,4873 (2008).[116] F. S. Pereira, Verification of ReFRESCO with theMethod of Manufactured Solutions , mathesis, InstitutoSuperior T´ecnico, Lisbon, Portugal (2012).[117] B. van Leer, Towards the Ultimate Conservative Differ-ence Scheme, Journal of Computational Physics ,
229 (1997).[118] F. F. Grinstein, A. A. Gowardhan, and A. J. Wach-tor, Simulations of Richtmyer-Meshkov Instabilities inPlanar Shock-Tube Experiments, Physics of Fluids (2011). [119] B. M. Haines, F. F. Grinstein, and J. R. Fincke, Three-Dimensional Simulation Strategy to Determine the Ef-fects of Turbulent Mixing on Inertial-Confinement-Fusion Capsule Performance, Physical Review E (2014).[120] D. Virk, F. Hussain, and R. M. Kerr, Compressible Vor-tex Reconnection, Journal of Fluid Mechanics304