Large-eddy simulation of Rayleigh-Taylor turbulence with compressible miscible fluids

Large-eddy simulation of Rayleigh-Taylor turbulence with compressible miscible fluids J. P. Mellado and S. Sarkara兲 Department of Mechanical and Aerospace Engineering, University of California San Diego, La Jolla, California 92093-0411 Y. Zhou Lawrence Livermore National Laboratory, University of California, Livermore, California 94551 共Received 25 March 2004; accepted 17 May 2005; published online 13 July 2005兲 Turbulence developed from Rayleigh-Taylor instability between two compressible miscible fluids in an unbounded domain is addressed in this paper. It is demonstrated that the turbulent Mach number in the turbulent core has an upper bound, independent of the density ratio under a broad range of initial mean configurations. The initial thermodynamic state of the system determines the amount of potential energy per unit mass involved in the turbulent mixing stage, and thus the characteristic level of turbulent fluctuations that is achievable is linked to the characteristic speed of sound such that the turbulent Mach number is limited. For the particular case of an ideal gas, this bound on the turbulent Mach number is found to be between 0.25 and 0.6, depending on the particular initial thermodynamic state. Hence, intrinsic compressibility effects 共those owing to large Mach number兲 are likely to be limited in the turbulent stage of a pure Rayleigh-Taylor problem. This result is confirmed by large-eddy simulations 共LES兲of systems with density jumps at the interface of 3:1, a density ratio for which there is extensive data available in the literature. The LES predictions of the mixing depth growth and overall mixing agree with results previously obtained in incompressible configurations with a negligibly small Mach number, and the data fully describing the Reynolds stresses and the budget of the 共resolved兲turbulent kinetic energy equation are provided. © 2005 American Institute of Physics.关DOI: 10.1063/1.1965130兴 I. INTRODUCTION Rayleigh-Taylor instability occurs when the interface between two fluids with different densities is subjected to a pressure force 共−⵱p兲pointing toward the heavy fluid, and this instability can eventually lead to a turbulent flow. The standard configuration is the hydrostatic equilibrium under the presence of a volumetric force g, directed toward the light fluid, ⵱p= ␳ g.共1兲 The problem was described by Rayleigh1and the first linear stability analysis was carried out by Taylor,2an analysis extended later to include diffusion effects3and, more recently, compressibility effects.4Interest in the topic reappeared strongly in 1980s, due to the inertial confinement fusion 共ICF兲programs.5–8 Literature covers theoretical,3,4,7,9 experimental,5,8,10,11 and numerical6,8,12–15 approaches to the problem. A detailed comparison between experiments and simulations has also been documented.16,17 Qualitatively, the flow evolves as follows.6,7 If the perturbation of the interface is small enough, then the linear analysis is valid and describes the exponential growth in this linear stage. Then follows a nonlinear stage, where asymmetric structures in the form of rising bubbles and falling spikes start forming owing to the baroclinic production of vorticity, organized in rings around the protruding fingers. This stage is shown in Fig. 1 共top row兲, where isosurfaces of the density field and the horizontal vorticity field are represented. The latter is defined by ␻ h 2= ␻ x 2+ ␻ y 2if the OZ axis is chosen parallel to the volumetric force. The vortex rings are slightly distorted in the vertical direction due to the shear between adjacent fingers. This organized distribution of coherent structures is eventually lost due to the nonlinear interaction among them, and more intertwined distributions of both fields are observed in the bottom row of Fig. 1. Larger structures appear, either by the amalgamation of smaller ones or by their presence in the initial condition. At the same time, KelvinHelmholtz instabilities appear at the sides of the fingers. By the end of this stage the memory of the initial conditions 共at least “small scales”兲can be potentially lost, with the turbulent stage taking over. Larger scales continue to be formed, viscosity having little effect on them and nonlinearity driving their energy toward the smallest scales. If the previous description holds and the statistics of the flow become truly independent of the initial conditions and viscosity, then the width of the mixing region between the two layers, denoted by h共t兲, in the low Mach number case usually considered in the literature, depends only on ␳ L, ␳ H, g, and t. Thus, dimensional analysis yields h/gt2=f共A兲 which, for the case AⰆ1, can be approximated by the linear relationship a兲Author to whom correspondence should be addressed. Telephone: 共858兲 534-8243. Fax: 共858兲534-7599. Electronic mail: [email protected] PHYSICS OF FLUIDS 17, 076101 共2005兲 1070-6631/2005/17共7兲/076101/20/$22.50 © 2005 American Institute of Physics17, 076101-1 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp h= ␣ Agt2,共2兲 the Atwood number being defined in terms of the ratio between the densities of the heavy fluid ␳ Hand the light fluid ␳ Las A= ␳ H/ ␳ L−1 ␳ H/ ␳ L+1.共3兲 This scaling leads to a temporal variation of the Reynolds number as Reh=hh ˙/ ␯ ⬃t3, ␯ being the kinematic viscosity, which implies that the Kolmogorov scale ␩ =共 ␯ 3/␧兲1/4, where ␧is the mean dissipation rate of turbulent kinetic energy per unit mass, continually decreases with time13 as ␩ ⬃t−1/4. Therefore, direct numerical simulation 共DNS兲is restricted in the temporal interval that can be used. As an alternative to DNS, the problem can be closed using numerics alone, for example, Youngs12 and Linden et al.16 solved the Euler equations with a monotone scheme. Cook et al.15 performed large-eddy simulation 共LES兲employing a subgrid eddy diffusivity that is strongly nonlinear, being dependent on the eighth-order derivative of the velocity. Here we adopt the LES strategy as well, but using a dynamic mixed model, already demonstrated to provide the temporal evolution of the large-scale three-dimensional fields with reasonable accuracy in different flows.18 The reader is referred to any of the several review articles available in the literature18–21 for a general description of the LES approach. A large extent of the research on the Rayleigh-Taylor problem has focused on the growth of the mixing zone. The total thickness his usually split into the bubble penetration distance hband the spike penetration height hs. According to the literature, ␣ bseems to be independent of the Atwood number, whereas ␣ sexhibits a slight increase with it, causing asymmetry in the flow for sufficiently large density ratios. A recent thorough review of the available data22 shows ␣ b varying between 0.03 and 0.08, the numerical simulations generally indicating smaller growth rates than the experimental results. In this regard, the strong influence that largescale initial conditions 共large meaning comparable to the final mixing thickness兲may have in the subsequent flow development has been reported by Linden et al.16 and Cook et al.,13 among others. However, an analysis of the role of compressibility in the turbulent stage of the Rayleigh-Taylor problem has not been reported. This paper considers this turbulent stage with miscible fluids in an unbounded domain. As mentioned before, work has been done on the effect of compressibility on the initial linear stage, an example for immiscible fluids being provided by Livescu,4but our paper is concerned with the stage of fully developed turbulence.A priori, these compressibility effects could appear if the speed of sound of the fluid is sufficiently reduced 共e.g., reducing the reference pressure of the system兲so that a characteristic Mach number is sufficiently large, and this is the motivation for the present study. It should be noted that the Richtmyer-Meshkov problem is fundamentally different in this respect, since the velocity scale, imposed there by the impulsive deposition of vorticity at the front, can be set independently of the thermodynamic state. A general overview containing the latest results on compressible turbulence can be found in the work by Chassaing et al.23 In the case of free shear flows, there is a strong intrinsic compressibility effect: the growth rate of the shear layer thickness and the turbulence intensity are significantly reduced as the Mach number increases. The cause has been the subject of a considerable amount of study during the past years,24–29 leading to the following picture: the production term in the turbulent kinetic energy equation reduces as a consequence of the decrease in the pressure-strain correlation, which diminishes the transfer of energy from the streamwise fluctuations to the cross-stream fluctuations. It is reasonable, then, to formulate the same questions for the Rayleigh-Taylor problem, where the input of energy is essentially different. Is the growth rate ␣ of the mixing depth ha function of the Mach number? And if so, is the physical mechanism the same as in the case of a shear flow? The focus of this paper is on intrinsic compressibility, density variations due to pressure variations, that is measured by the flow Mach number. There is no mean flow in the problem, and the compressibility associated with the turbulent fluctuations is measured by the turbulent Mach number, Mt=q 具c典,共4兲 where q=冑2K,Kbeing the turbulent kinetic energy, and 具c典 the average speed of sound. The average of any variable ␾ is written as 具 ␾ 典and it is computed as a plane average at a fixed inhomogeneous location z. They denote Reynolds averages for quantities ␾ per unit volume and Favre averages for quantities ␾ per unit mass, unless otherwise stated. The quantity Mt共z兲varies across the mixing zone. The value of Mtin the core is of interest where the speed of sound is that of the mixed interfacial fluid, c0. The compressibility of the fluid also modifies the stability inside of the fluid layers with respect to the incompressible situation. This stability is measured by the buoyancy-frequency30,31 FIG. 1. Isosurfaces of density field 共left兲, ␳ =2, and of horizontal vorticity 共right兲, ␻ h=0.2 of maximum. Top—t=2.5 ␶ 共nonlinear stage兲; bottom—t =10 ␶ 共turbulent stage兲. The characteristic time ␶ =冑 ␦ ␳ ,0/共Ag兲is defined using the initial thickness of the mixing layer. Gravity is acting downward. Only 1/4 of the whole domain is shown for clarity. 076101-2 Mellado, Sarkar, and Zhou Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp N2=g ␪ d ␪ dz ,共5兲 with ␪ being the potential temperature. For an ideal gas, it follows that N2=g ␳ p1/ ␥ d dz 冉 p1/ ␥ ␳ 冊 ,共6兲 where pis the pressure, ␳ is the density, and ␥ is the ratio of specific heats. The coordinate zincreases in the opposite direction of the volumetric force g. The buoyancy-frequency has been written in terms of the pressure and the density because it is more suitable for the analysis to be presented below. The sign of N2determines the static stability of the configuration. Therefore, in addition to the Rayleigh-Taylor instability at the interface between the two layers originating from the density jump, each layer can individually be buoyancy-stable if N2⬎0 or buoyancy-unstable if N2⬍0. Although the topic of the investigation reported here is intrinsic compressibility in a turbulent flow, the buoyancycompressibility coupling may affect the evolution of Rayleigh-Taylor turbulence and is included in the analysis. Note that one could also study the evolution of small 共nonturbulent兲fluctuations in cases with positive, negative, and zero values of N2, but this is not what is being done here. The paper is organized as follows. In Sec. II, the problem is defined and a theoretical analysis of the role of compressibility, giving definite bounds on the turbulent Mach number, is presented. The dynamic mixed model used to close the subgrid terms in the filtered equations is presented in Sec. III along with the general formulation of the problem. Sections IV and V contain the results from the LES, confirming the theoretical predictions obtained in Sec. II. Several statistics, such as mean profiles, Reynolds stresses, and the budget of the resolved turbulent kinetic energy, are also discussed. II. COMPRESSIBILITY OF THE TURBULENCE We consider the hydrostatic equilibrium of a layer of heavy fluid on top of a layer of lighter fluid. Choosing the OZ axis parallel to and in the opposite direction of the constant volumetric force g, the equilibrium is determined by dp dz =− ␳ g.共7兲 The function ␳ = ␳ 共p兲has to be provided in order to solve the problem, and different configurations will be explored below. The equation of state for a mixture of two ideal gases reads p ␳ =R0T W=R0T 冉 YL WL +YH WH 冊 ,共8兲 where Yidenotes the species mass fraction of the fluid iand R0is the universal gas constant. The molecular weight of the mixture, W, has been written in terms of the molecular weight of the heavy fluid and the light fluid, WHand WL, respectively. The intrinsic compressibility is determined by the turbulent Mach number, defined by Eq. 共4兲. The local speed of sound of an ideal gas is given by c=冑 ␥ R0T W=冑 ␥ p ␳ ,共9兲 where ␥ , the ratio of specific heats, lies in the range 1⬍ ␥ ⬍5/3. With respect to the turbulent kinetic energy, turbulent theory predicts and experiments confirm that, in the incompressible limit, a self-similar state is achieved after a sufficiently long time. A thorough dimensional analysis is presented in the work by Cook et al.13 关see Eq. 共6兲therein兴, which yields q0= ␤ 冑Agh 共10兲 as a characteristic scale of the turbulent velocity fluctuations at each time. In this expression, Ag represents the constant external force per unit mass, with Abeing the Atwood number and h共t兲the thickness of the increasing mixing depth. With miscible fluids, an 共smaller兲effective Atwood number instead of the nominal Atwood number Awould be more appropriate, as discussed by Cook et al.15 The coefficient ␤ , of order 1, has to be provided by the experimental data. A priori, it is not known how Eq. 共10兲is modified by compressibility. Hence an analysis based on the available potential energy is preferred to study the problem and Eq. 共10兲is not used. However, the fact that the fluid velocity is determined by the potential energy per unit mass 共set by the initial thermodynamic variables兲links the maximum of q0 over the time to c0, the characteristic speed of sound, such that Mtis bounded from above. The consequence is that Mt, small at early times because of an initial state that is steady, might not become large enough for intrinsic compressibility effects to be strong, so that Eq. 共10兲might be reasonably well satisfied even in the compressible case of RayleighTaylor turbulence. This section is devoted to obtain analytically an upper bound on Mtfor three different initial configurations ␳ 共p兲 leading to the three different types of static stability: 共1兲constant ␳ /pin each layer, in Sec. II A, which is buoyancystable 共N2⬎0兲;共2兲constant ␳ /p1/ ␥ in each layer, in Sec. II B, which is buoyancy-neutral 共N2=0兲,apure Rayleigh-Taylor problem; 共3兲constant ␳ in each layer, in Sec. II C, which is buoyancy-unstable 共N2⬍0兲. A. Two buoyancy-stable layers This section considers two layers which are buoyancystable, i.e., a relation ␳ 共p兲in each of them such that the buoyancy-frequency is positive. The particular configuration, shown in Fig. 2, is formed by two layers, each composed of a pure homogeneous fluid. The ratio between the molecular weight and the temperature of the mixture varies between two well-defined levels, WL/TLat the bottom and WH/TH 共larger兲at the top. This setup implies a relation ␳ ⬀p, and Eq. 共6兲provides N2⬎0. This case can be easily set up, for instance, by depositing a layer of heavy fluid on top of a second lighter fluid at the same temperature. 076101-3 Large-eddy simulation of Rayleigh-Taylor Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp The pressure is obtained by substituting the expression for the density from the equation of state, Eq. 共8兲, into the balance equation, Eq. 共7兲, which yields p共z兲=p0exp 冉 −gWH R0TH 冕 0 zW共 ␨ 兲/WH T共 ␨ 兲/TH d ␨ 冊 ,共11兲 p0being the value of the pressure at the middle plane, z=0. This equation shows that the pressure decays exponentially with increasing height with a characteristic length scale LH=R0TH gWH ,共12兲 the scale-height,32 which is directly related to the speed of sound, in particular, proportional to cH 2/g. There is a similar length scale associated with the light fluid LL=R0TL gWL .共13兲 The density profile is then computed from the equation of state, which yields ␳ 共z兲= ␳ 0 +W共z兲/WH T共z兲/TH exp 冉 −1 LH 冕 0 zW共 ␨ 兲/WH T共 ␨ 兲/TH d ␨ 冊 ,共14兲 ␳ 0 +being related to p0by ␳ 0 +=p0WH R0TH =p0 gLH .共15兲 Either ␳ 0 +or p0has to be provided. It is interesting to note that they are also linked to the amount of mass deposited on the top, mH, a relation obtained by integrating the density profile along the OZ axis. Sketches of the different profiles are presented in Fig. 2. The density has a jump at the interface z=0 from ␳ 0 −at z=0−to ␳ 0 +at z=0+given by ␳ 0 + ␳ 0 −=WHTL WLTH =LL LH ,共16兲 where the last equality follows from the definition of the two scale-heights. This equation shows that, if ␳ 0 +/ ␳ 0 −⬎1, a condition required to have Rayleigh-Taylor instability, then TH/WH⬍TL/WLand, consequently, LH⬍LL. 1. Energy analysis In order to estimate the characteristic intensity of turbulent fluctuations q0, the available potential energy of the system is calculated, since this constitutes an upper limit to the kinetic energy released to the flow. The length over which the profile jumps from WL/TLto WH/THis assumed to be small compared to the scale-heights, the latter being defined by Eqs. 共12兲and 共13兲. The density profile is ␳ 共z兲= 再 ␳ 0 +e−z/LH,z⬎0 ␳ 0 −e−z/LL,z⬍0 冎 共17兲 according to Eq. 共14兲. A mass of light fluid mLis deposited over a domain depth DLand, on top of it, a mass of heavy fluid mHover a depth DH. These four lengths and two masses completely define the system, since mH= 冕 0 DH ␳ 共z兲dz = ␳ 0 +LH共1−e−DH/LH兲, 共18兲 mL= 冕 −DL 0 ␳ 共z兲dz = ␳ 0 −LL共eDL/LL−1兲 can be used to obtain ␳ 0 +and ␳ 0 −. The equation of state provides the pressure, if desired. This system is not in stable equilibrium, and, either slowly by diffusion or more rapidly by Rayleigh-Taylor instability, it evolves toward a totally mixed configuration. A final buoyancy-stable configuration with TF/WFconstant is assumed. The final density distribution is ␳ f共z兲= ␳ Fe−共z+DL兲/LF,共19兲 where the density at the bottom, ␳ F, is the final density at z =−DLand the scale-height of the final state is defined by LF=R0TF gWF .共20兲 The constants ␳ Fand LFneed to be determined. One equation is the conservation of total mass mH+mL= ␳ FLF共1−e−共DL+DH兲/LF兲.共21兲 A second equation is derived from the conservation of total mass of each component separately. Working with the heavy fluid, its mass fraction is given, in general, by YH共z兲=1−WL/W共z兲 1−WL/WH ,共22兲 and the mass of heavy fluid is the integral of the profile ␳ 共z兲YH共z兲between −DLand DH, which yields the relation mH mH+mL =1− ␺ 1−WL/WH ,共23兲 with ␺ =WL mH+mL 冕 −DL DH ␳ Wdz.共24兲 Conservation of the mass of heavy fluid, mH, implies then that the function ␺ is equal between the initial and the final states. For the initial configuration, which is buoyancy-stable and isothermal in each of the two layers, it is FIG. 2. Profiles of the thermodynamic variables in a system formed by two layers of homogeneous fluids, which is Rayleigh-Taylor unstable. 076101-4 Mellado, Sarkar, and Zhou Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp ␺ 1,i= 共eDL/LL−1兲+共1−e−DH/LH兲TLLH THLL eDL/LL−e−DH/LH,共25兲 and for the final distribution ␺ 1,f=WL WF .共26兲 Thus, the conservation of mass of the heavy species provides the equation WL WF = ␺ 1,i,共27兲 and from the definitions of the scales heights, Eqs. 共13兲and 共20兲, we have LF LL =WLTF WFTL =TF TL ␺ 1,i.共28兲 Equations 共21兲and 共28兲provide ␳ F/ ␳ 0 +, and LF/LLas a funcion of DH/LL,DL/LL, and LH/LL, once TF/TLand TH/TL are known. The potential energy of the fluid particle located at zis given by ␳ 共z兲gz, the reference level chosen at z=0, and the total amount of potential energy Vof a certain region in space is the integral of ␳ 共z兲gz over that region. Hence, for the heavy fluid layer of the initial configuration we have VH= ␳ 0 +gLH 2关1−e−DH/LH共1+DH/LH兲兴 =mHgLH− ␳ 0 +gLHDHe−DH/LH.共29兲 The light fluid region yields VL= ␳ 0 −gLL 2关e+DL/LL共1−DL/LL兲−1兴 =mLgLL− ␳ 0 −gLLDLeDL/LL.共30兲 Lastly, the final configuration provides VF= ␳ FgLF 2e−DL/LF关e+DL/LF共1−DL/LF兲−e−DH/LF共1+DH/LF兲兴 =共mH+mL兲gLF− ␳ FgLF共DL+DHe−共DH+DL兲/LF兲.共31兲 The difference between the total initial potential energy and the final one is the available potential energy, Ep=VH+VL−VF=Ep,t+Ep,i,共32兲 which has been decomposed as the sum of two parts Ep,t=共mHgLH+mLgLL兲−共mH+mL兲gLF, 共33兲 Ep,i= ␳ FgLF共DL+DHe−共DH+DL兲/LF兲− ␳ 0 +gLH共DLe+DL/LL −DHe−DH/LH兲. It is the specific 共per unit mass兲available potential energy Ep 共mH+mL兲共34兲 which determines the velocity attained by the fluid. Equations 共25兲and 共28兲, under the assumption TH=TL, allow to write Ep,t 共mH+mL兲gLL = ␺ i,1 −LF LL = 冉 1−TF TL 冊 ␺ i,1,共35兲 where gLLis used to normalize the specific potential energy. This is a thermal contribution to the change in potential energy, nonzero only in the case of change in temperature between the initial and the final states. The second part, which is nonzero even for an isothermal case, can be written as Ep,i 共mH+mL兲gLL = ␾ 1 冉 DL LL ,DH LL ,LH LL 冊 ,共36兲 where the function ␾ 1is defined by ␾ 1=共DL/LL兲eDL/LF+共DH/LL兲e−DH/LF eDL/LF−e−DH/LF −共DL/LL兲eDL/LL+共DH/LL兲e−DH/LH eDL/LL−e−DH/LH,共37兲 and LF/LLis given by Eq. 共28兲in terms of the length ratios of the problem shown as the arguments of ␾ 1. The characteristic intensity of the turbulent fluctuations is given by q0/冑gLL=冑2 ␾ 1. The characteristic speed of sound c0is required at this stage in order to construct a characteristic turbulent Mach number. The speed of sound in the final mixed configuration 共miscible fluids are being considered兲, ␥ R0TF/WF, is used for this purpose, which normalized by gLLprovides c0 2 gLL = ␥ LF LL 共38兲 according to the definition of the final scale-height, Eq. 共20兲. Hence, a characteristic turbulent Mach number is given by Mt,0 =冑2 ␾ 1 ␥ 共LF/LL兲,共39兲 which is a function of DL/LL,DH/LL, and LH/LL. The mixing region 共−DL,DH兲must now be estimated. A first estimate in the turbulent case for the upper limit of the final mixing region is given by the distance from the initial density jump until the point in the upper layer where the density becomes equal to ␳ 0 −. This reasoning yields DH=LHln LL LH .共40兲 Similarly, the turbulent motion can develop toward the lower layer at most until the downward position z=−DLat which the initial density profile equals the value ␳ 0 +. This distance is given by DL=LLln LL LH .共41兲 However, bubbles expand as they rise, and this could imply larger values of DHand DLif turbulent mixing is not fast enough to eliminate the density gradients through molecular diffusion. We explore this other limit now. Consider pure fluid bubbles that rise through the upper layer without turbulent mixing. Close to the interface there are bubbles 076101-5 Large-eddy simulation of Rayleigh-Taylor Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp with density ␳ b= ␳ 0 −surrounded by a fluid with density ␳ 0 +. These bubbles expand as they rise because the surrounding pressure is lower, and their density diminishes. If this process is adiabatic, then ␳ b⬀p1/ ␥ , a smaller rate of decrease with height than the surrounding density, which varies proportionally to the pressure. Therefore, at a certain height, the density of the rising bubble is no longer smaller than the surrounding density and the buoyancy force is zero. This height determines the upper limit of the mixing zone and is given by the equation ␳ 0 − 冉 p共DH兲 p0 冊 1/ ␥ H = ␳ 0 + 冉 ␳ 共DH兲 ␳ 0 + 冊 ,共42兲 where the density and the pressure follow the known exponential decay given by Eq. 共17兲. The solution is DH= ␥ H ␥ H−1LHln LL LH .共43兲 Similar reasoning provides an estimate for DL, DL= ␥ L ␥ L−1LLln LL LH .共44兲 These two lengths are equal to those defined by Eqs. 共40兲 and 共41兲multiplied by the same constant ␥ /共 ␥ −1兲, greater than one, when a mean value of ␥ is considered. Therefore, we can consider the mixing zone given by 共− ␣ DL, ␣ DH兲, with DLand DHgiven by Eqs. 共43兲and 共44兲, and this approach covers also the former estimate of the thickness of the mixing layer for an intermediate value of ␣ . The parameter ␣ varies between 0 and 1, increasing in time corresponding to the evolution of h共t兲. It has to be noted that the function ␾ 1 gives the maximum specific available potential energy that would be released if the instantaneous mixing thickness h共t兲 would be held constant sufficiently long to allow complete mixing. The actual potential energy being released is smaller because the characteristic time for mixing inside the region is the same as the characteristic time to engulf new mass into the mixing zone, namely, the turbulent time scale determined by the turbulent kinetic energy and the rate of turbulent dissipation. To conclude the analysis, the maximum of Mt,0, Eq. 共39兲, along the curves DL/LL= ␣ 关 ␥ L/共 ␥ L−1兲兴ln共LL/LH兲and DH/LL= ␣ 关 ␥ H/共 ␥ H−1兲兴共LH/LL兲ln共LL/LH兲for all values of ␣ in the range 0⬍ ␣ ⬍1 has to be calculated. When a mean value of ␥ between 1 and 5/3 is substituted in that equation, the characteristic turbulent Mach number is found to be bounded by 0.62. Note that this value is a conservative estimate, because a considerable amount of the available potential energy 共of the order of 50% according to the data regarding the incompressible case兲is dissipated and does not contribute to the turbulent kinetic energy. B. Two buoyancy-neutral layers This section considers two layers which are buoyancyneutral, i.e., a relation ␳ 共p兲in each of them such that the buoyancy frequency is zero. Equation 共6兲with N=0 implies ␳ ⬀p1/ ␥ . Hence, the thermodynamic state is defined by the group p/ ␳ ␥ varying between two uniform values 共p/ ␳ ␥ 兲Land 共p/ ␳ ␥ 兲H, the profile W/Tbeing as required by the equation of state in order to obtain such a distribution. When the layers are pure fluids, then the entropy is constant inside each layer. This initial condition has been studied in the past by Chen et al.,33 although they analyzed the two-dimensional problem with immiscible fluids, fundamentally different from the fully turbulent configuration with miscible fluids considered here. For notational convenience, the profile p1/ ␥ / ␳ is denoted by s共z兲, which then varies from sLat the bottom to sHat the top. Substituting ␳ =p1/ ␥ /sinto the balance equation, Eq. 共7兲, the pressure distribution is p共z兲=p0 冉 1− ␥ −1 ␥ g p0 共 ␥ −1兲/ ␥ sH 冕 0 zd ␨ s共 ␨ 兲/sH 冊 ␥ /共 ␥ −1兲 ,共45兲 where p0is again the pressure at the center plane. The ratio of specific heats ␥ has been taken as constant. The characteristic length scale is LH=p0 共 ␥ −1兲/ ␥ sH g=p0 g ␳ 0 +=R0T0 + W0 +g,共46兲 having defined ␳ 0 +=p0 1/ ␥ /sH. Similarly, the scale-height of the lower layer is LL=p0 共 ␥ −1兲/ ␥ sL g=p0 g ␳ 0 −=R0T0 − W0 −g.共47兲 The density profile is ␳ 共z兲= ␳ 0 +1 s共z兲/sH 冉 1− ␥ −1 ␥ 1 LH 冕 0 zd ␨ s共 ␨ 兲/sH 冊 1/共 ␥ −1兲 .共48兲 The density has a jump at the interface z=0 from ␳ 0 −at z =0−to ␳ 0 +at z=0+given by ␳ 0 + ␳ 0 −=sL sH =LL LH .共49兲 It can be observed that the derivation is parallel to that of Sec. II A. In fact, that section considers ␳ ⬀p, which is the limit ␥ →1 of the case considered here, ␳ ⬀p1/ ␥ . Since ␥ is always of order one, this convergence suggests that there is no qualitative change in the bounding of Mtwith respect to the buoyancy-stable configuration presented in Sec. II A. It is also observed that the density and the pressure fall down to zero along the upper layer on a distance of order O关 ␥ /共 ␥ −1兲LH兴. 1. Energy analysis The calculation of specific potential energy is now presented, neglecting once more the thickness of the initial interface with respect to the scale heights LHand LL, Eqs. 共46兲 and 共47兲, respectively. The initial mean density profile is 076101-6 Mellado, Sarkar, and Zhou Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp ␳ 共z兲= 冦 ␳ 0 + 冉 1− ␥ H−1 ␥ H z LH 冊 1/共 ␥ H−1兲 ,z⬎0 ␳ 0 − 冉 1− ␥ L−1 ␥ L z LL 冊 1/共 ␥ L−1兲 ,z⬍0. 冧 共50兲 It can be verified that Eq. 共50兲tends toward Eq. 共17兲in the limit ␥ →1. In a fashion similar to that of the previous sections, a mass of light mixed fluid mLis deposited over a domain depth DL, mL= 冕 −DL 0 ␳ 共z兲dz = ␳ 0 −LL 冋 冉 1+ ␥ L−1 ␥ L DL LL 冊 ␥ L/共 ␥ L−1兲 −1 册 , 共51兲 and, on top of it, a mass of heavy mixed fluid mHover a depth DH, mH= 冕 0 DH ␳ 共z兲dz = ␳ 0 +LH 冋 1− 冉 1− ␥ H−1 ␥ H DH LH 冊 ␥ H/共 ␥ H−1兲 册 . 共52兲 The temperature is supposed to be uniform inside each layer, as done by Chen et al.33 and as assumed in Sec. II A as well. This condition implies that the two layers are not composed of pure fluid, but that the molecular weight varies, according to the equation of state, as 1 W共z兲= 冦 1 W0 + 冉 1− ␥ H−1 ␥ H z LH 冊 ,z⬎0 1 W0 − 冉 1− ␥ L−1 ␥ L z LL 冊 ,z⬍0. 冧 共53兲 The same reasoning as in the preceding section is now applied. The final state is considered to be a complete homogeneous mixture, buoyancy-stable, which provides the minimum final potential energy and therefore a conservative estimate for the potential energy released into the flow. The density is given by Eq. 共19兲and the conservation of total mass is expressed again by Eq. 共21兲. The corresponding function ␺ 2,i, defined by Eq. 共24兲, is given by W0 − WL ␺ 2,i= ␥ L 2 ␥ L−1 冋 冉 1+ ␥ L−1 ␥ L DL LL 冊 共2 ␥ L−1兲/共 ␥ L−1兲 −1 册 + ␥ H 2 ␥ H−1 冋 1− 冉 1− ␥ H−1 ␥ H DH LH 冊 共2 ␥ H−1兲/共 ␥ H−1兲 册 LHTL LLTH 冉 1+ ␥ L−1 ␥ L DL LL 冊 ␥ L/共 ␥ L−1兲 − 冉 1− ␥ H−1 ␥ H DH LH 冊 ␥ H/共 ␥ H−1兲.共54兲 It can be again verified that Eq. 共54兲recovers Eq. 共25兲in the limit ␥ →1. Then the conservation of the mass of heavy fluid implies, using Eq. 共26兲, WL WF = ␺ 2,i.共55兲 The final scale-height is given by Eq. 共20兲, which can be written explicitly as LF LL =TF TL W0 − WF =TF TL W0 − WL ␺ 2,i共56兲 with Eq. 共47兲. Calculating the potential energies as explained in Sec. II A 共TH=TLis again assumed兲, the normalized specific available potential energy due to temperature changes becomes Ep,t 共mH+mL兲gLL =W0 − WL ␺ 2,i−LF LL = 冉 1−TF TL 冊 W0 − WL ␺ 2,i.共57兲 The remaining part is written as Ep,i 共mH+mL兲gLL = ␾ 2 冉 DL LL ,DH LL ,LH LL 冊 ,共58兲 where ␾ 2is defined by ␾ 2=共DL/LL兲eDL/LF+共DH/LL兲e−DH/LF eDL/LF−e−DH/LF− DL LL 冉 1+ ␥ L−1 ␥ L DL LL 冊 ␥ L/共 ␥ L−1兲 +DH LL 冉 1− ␥ H−1 ␥ H DH LH 冊 ␥ H/共 ␥ H−1兲 冉 1+ ␥ L−1 ␥ L DL LL 冊 ␥ L/共 ␥ L−1兲 − 冉 1− ␥ H−1 ␥ H DH LH 冊 ␥ H/共 ␥ H−1兲.共59兲 076101-7 Large-eddy simulation of Rayleigh-Taylor Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp The characteristic intensity of the turbulent fluctuations is given by q0/冑gLL=冑2 ␾ 2. The characteristic speed of sound c0is now needed. The squared speed of sound in the final mixed configuration is ␥ R0TF/WF, which normalized by gLLprovides c0 2 gLL = ␥ LF LL 共60兲 according to the definition of the final scale-height, Eq. 共20兲. Finally, a characteristic turbulent Mach number is given by Mt,0 =冑2 ␾ 2 ␥ 共LF/LL兲,共61兲 similar to Eq. 共39兲but with the corresponding function ␾ 2. The final mixing region 共−DL,DH兲has to be estimated. A first estimate for the upper limit of the final mixing region is given by the distance from the initial density jump until the point in the upper layer where the density becomes equal to ␳ 0 −. This reasoning yields DH= ␥ H ␥ H−1LH 冋 1− 冉 LH LL 冊 ␥ H−1 册 .共62兲 Similarly, the turbulent motion can develop toward the lower layer until the downward position z=−DLat which the initial density profile equals the value ␳ 0 +. This distance is given by DL= ␥ L ␥ L−1LL 冋 冉 LL LH 冊 ␥ L−1 −1 册 .共63兲 Now, the role of the neutral stability of each initial layer has to be considered, as it was done in Sec. II A. The same reasoning here implies that pure bubbles continually rise through all the available upper layer, whose limit is given by the point where the density becomes zero. From Eq. 共50兲, this singularity occurs at a distance DH=LH ␥ H ␥ H−1.共64兲 Therefore, we consider the mixing zone given by 共− ␣ DL, ␣ DH兲, with ␣ varying between 0 and 1, DHgiven by Eq. 共64兲, and DL=DH 冉 LL LH 冊 ␥ .共65兲 Such a choice for the mixing zone covers that described by Eqs. 共62兲and 共63兲for some intermediate value of ␣ , and includes the whole upper layer when ␣ =1. Note that Eq. 共65兲 recovers the ratio DL/DHof the buoyancy-stable initial configuration, Sec. II A, in the limit ␥ →1. The maximum of Mt,0, Eq. 共61兲, along the previous curves DL/LL= ␣ 关 ␥ /共 ␥ −1兲兴共LL/LH兲共 ␥ −1兲and DH/LL = ␣ 关 ␥ /共 ␥ −1兲兴共LH/LL兲for 0⬍ ␣ ⬍1, is now to be calculated. When a mean value of ␥ between 1 and 5/3 is substituted in this equation, the characteristic turbulent Mach number is bounded by 0.57. The reader is again reminded that all the potential energy is assumed to be converted to turbulent kinetic energy without any loss to molecular dissipation; therefore, this upper bound of Mt,0 is a conservative value. C. Two buoyancy-unstable layers This section considers two layers which are buoyancyunstable, i.e., a relation ␳ 共p兲in each layer such that the buoyancy-frequency is negative. The particular configuration, shown in Fig. 3, is formed by two layers of constant density with a jump from ␳ Lto ␳ H共larger兲at the interface z=0, the profile T/Wvarying as required by the equation of state in order to obtain such a distribution. Equation 共6兲with ␳ constant implies N2⬍0 inside each layer. This case might not be set up as easily as the one considered in the previous sections 共it is unstable兲, and it will be not considered with such generality, but it is of interest for several reasons. First, because it represents a hydrostatic profile that is unstable. Second, the Atwood number based on the mean density profile is constant instead of decreasing with time 共although the effective Atwood number increases, as it is later shown in the simulations兲. Third, there are data available in the literature for the incompressible limit of this configuration that can be utilized to validate the LES results. Given ␳ 共z兲, the integration of Eq. 共7兲yields p共z兲=p0 冉 1− ␳ Hg p0 冕 0 z ␳ 共 ␨ 兲 ␳ H d ␨ 冊 .共66兲 The scales that characterize the thermodynamic variation are LH=p0 ␳ Hg,共67兲 similar to Eq. 共15兲with ␳ 0 +in that equation corresponding to ␳ Hin Eq. 共67兲, and LL=p0 ␳ Lg.共68兲 The required ratio T共z兲/W共z兲is calculated from the equation of state, T共z兲 W共z兲= 冉 T W 冊 0 +1− 冕 0 z ␳ 共 ␨ 兲/ ␳ Hd ␨ /LH ␳ 共z兲/ ␳ H ,共69兲 where 冉 T W 冊 0 + =p0 R0 ␳ H .共70兲 The profiles are sketched in Fig. 3. The constancy of density over each 共heavy or light兲fluid column can be achieved by either decreasing the temperature or increasing the molecular weight as zincreases. Needless to say, the distribution W共z兲 FIG. 3. Profiles of the thermodynamic variables in a system formed by two layers of constant density, which is Rayleigh-Taylor unstable. 076101-8 Mellado, Sarkar, and Zhou Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp determines the two species to be used in the problem, since the relation WL⬍W共z兲⬍WH共71兲 holds. In this configuration, the buoyancy frequency, Eq. 共6兲,is negative, and this instability implies that the whole domain will be involved in the mixing process. However, it is immediately observed from the linear decrease in the pressure that the maximum domain length cannot be larger than O共LH兲, given by Eq. 共67兲, because the pressure reaches zero at that limit 共we would need either zero absolute temperature or a species of infinite molecular weight兲. The problem is studied up to this instant of time. 1. Energy analysis An analysis based on the specific available potential energy similar to that of the preceding section is performed again. The thickness of the density jump is neglected and the density profile is ␳ 共z兲= 再 ␳ H,z⬎0 ␳ L,z⬍0. 冎 共72兲 We consider a mixing region between −DLand DH. Based on the available incompressible data,13,15 the density profiles at each time are estimated approximately by a linear variation from ␳ Lat z=−DLto ␳ Hat z=+DH. This assumption is later verified in the simulation for the compressible case. Conservation of mass inside this mixing region then implies DH=DL. Experimental results show, however, that spikes evolve significantly faster than bubbles when the density ratio is large enough,8therefore, the analysis presented here for this configuration is applicable only to relatively small Atwood numbers, Aⱗ0.5. Utilizing the same formulation as in the preceding section, the potential energies are VH= ␳ HgDH 2/2, VL=− ␳ LgDH 2/2, 共73兲 VF=共 ␳ H− ␳ L兲gDH 2/3, the available potential energy is Ep=VH+VL−VF, and the normalized specific available potential energy is Ep 共mH+mL兲gLL = ␾ 3 冉 DH LL ,LH LL 冊 ,共74兲 where ␾ 3=1 6ADH LL =1 6 LL/LH−1 LL/LH+1 DH LL .共75兲 A characteristic turbulent intensity is then given by q0/冑gLL=冑2 ␾ 3. On the other hand, an estimate of the characteristic speed of sound at the center plane, c0, can be obtained now using a mean value 共p0/ ␳ H+p0/ ␳ L兲/2 in Eq. 共9兲and the definitions Eqs. 共67兲and 共68兲, which yields c0 2 gLL = ␥ 2 冉 1+LH LL 冊 .共76兲 The minimum of the speed of sound corresponds to the upper limit of the domain, where the pressure is minimum, and it can be made as small as desired by reducing p0. However, the turbulent motion develops around the middle plane and it is the value of cthere that matters. Hence, a characteristic turbulent Mach number is given by Mt,0 =冑4 ␾ 3 ␥ 共1+LH/LL兲,共77兲 which is a function of DH/LLand LH/LL, the function ␾ 3 defined by Eq. 共75兲. The final size of the mixing region DH/LLmust now be expressed in terms of the given configuration LH/LL = ␳ L/ ␳ H. As already mentioned before, the upper domain height DHcannot exceed LHso that the maximum value that the mixing zone can reach is precisely DH LL =LH LL .共78兲 The maximum of the function Mt,0, defined in Eq. 共77兲, along the curve DH/LL=f共LH/LL兲gives an upper bound for the turbulent Mach number of 0.25, having used again a mean value of ␥ between 1 and 5/3. Note that ␾ 3, Eq. 共75兲, increases monotonically with the mixing region DH/LL, and the final value of DH/LLused above maximizes the intermediate states ␣ DH/LL,0⬍ ␣ ⬍1. D. Summary of the analytical results The major result of Sec. II is that the turbulent Mach number has an upper bound, independent of the density ratio under a wide range of initial thermodynamic configurations, which may not be large enough for intrinsic compressibility effects to be important in Rayleigh-Taylor turbulence. This result holds for three different choices of static stability inside each layer: stable, neutral, and unstable. An assumption underlying the analysis is that the flow is fully turbulent. In this respect, if a large scale perturbation, O共LH兲, is imposed initially at the interface then there could be compressibility effects as a blob of pure fluid rises/falls into the opposite pure fluid layer; this is not the case studied here. Also, miscible fluids subject to a turbulent flow have been considered, which allows the characteristic speed of sound to be estimated by that of the mixed fluid. This aspect of the problem could be different in the case of immiscible fluids. Another assumption is that of an ideal gas. The fundamental cause of the limitation of Mtis independent of this latter assumption, however, the particular upper bound of the turbulent Mach number found in the analysis depends on the details of each equation of state. It is also worth noticing that Secs. II A and II B consider 076101-9 Large-eddy simulation of Rayleigh-Taylor Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp come, ␦ ␳ starts to accelerate again, seemingly toward the quadratic scaling. The linear fit to the data beyond t=10 ␶ provides ␦ ˙ ␳ /冑Ag ␦ ␳ ,0 = 0.25共t/ ␶ − 3.96兲.共101兲 Again, intrinsic compressibility does not influence the usual scaling law of the growth of ␦ ␳ that neglects any Mach number influence. At the nondimensional time of 16, the thickness of the mixing layer is ␦ ␳ *⯝0.61共h ␳ *⯝0.72兲. The simulation is then stopped because the vertical extent is comparable to the homogeneous dimensions 共recall that the domain size is LH⫻LH⫻2LH兲and the number of structures in the domain does not allow a reasonable average. Figure 11 plots the density profile across the flow for two different times, t=10 ␶ and t=15 ␶ . The vertical coordinate is scaled with the mixing depth ␦ ␳ , and the collapse of the two curves suggests an approach to similarity. There is a small deviation, however, at the two fronts; it is not clear if the deviation is due to the LES model or truly due to the lack of self-similarity. This deviation will be shown to be stronger in the Reynolds stresses profiles. The last density-related quantity to be analyzed is the density fluctuation variance, cast in terms of an effective Atwood number, Eq. 共92兲. Its temporal evolution at the center plane is shown in Fig. 12. The peak occurring at t ⯝2.5 ␶ represents an increase in segregation taking place in the nonlinear stage 共minimum of the mixing parameters in Fig. 9兲. After that, the sustained increase of density fluctuations beyond transition is due to the unstable stratification 共negative buoyancy-frequency兲of the upper and the lower layers, consistently the opposite tendency from the earlier case, Fig. 5, where the stratification was stable. This result is also consistent with the sustained decrease in the mixing parameters shown in Fig. 9. The information given by Fig. 12 is relevant for the topic investigated in this paper. First, miscible fluids are characterized by an effective Atwood number significantly smaller than immiscible fluids. For immiscible fluids, i.e., completely segregated states, the effective Atwood number is approximately15 equal to A, or possibly larger because of the unstable stratification in each layer. Miscible fluids are characterized by a smaller value in the turbulent stage, in particular, we obtain Ae=0.5A, which is in agreement with the incompressible case considered by Cook et al.15 Figure 13, the visual counterpart of the information conveyed by the curve in Fig. 12, shows that the late-time mixing 共bottom panel兲is turbulent; there is a central region of intermediate fluctuating values of density, and the probability of bubbles or spikes enclosing pure fluid and cleanly penetrating in the opposite layer is small. Turbulent mixing also justifies the choice of the characteristic speed of sound in the theoretical analysis to be a mean value of the speed of sound in the turbulent core. Second, the effect of compressibility can be directly observed comparing Fig. 12 with Fig. 19 in the work of Cook et al.15 These authors considered the incompressible case and they observed a constant behavior of Aefor late times due to stable stratification. Compressibility imposes a negative buoyancy-frequency and thus unstable stratification, which FIG. 10. Mixing width ␦ ␳ 共dash-dotted line兲and its time derivative ␦ ˙ ␳ 共dashed line兲as a function of time. The linear fit to ␦ ˙ ␳ ,Eq.共101兲, is also shown as a solid line. FIG. 11. Mean density profile using self-similar variables: 共-兲t =10冑 ␦ ␳ ,0/共Ag兲;共---兲t=15冑 ␦ ␳ ,0/共Ag兲. FIG. 12. Temporal evolution of the normalized intensity of the density fluctuations or effective Atwood number, Eq. 共92兲, at the center plane. FIG. 13. Density fields. Top—t=2.5 ␶ 共nonlinear stage兲; bottom—t=10 ␶ 共turbulent stage兲. Gravity is acting downward. 076101-16 Mellado, Sarkar, and Zhou Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp results in the increase of the density fluctuations. The increase in the effective Atwood number due to this buoyancycompressibility coupling is of the order of 20%. E. Reynolds stresses The Reynolds stress tensor is split into resolved, and subgrid parts, Rij =Rij r+Rij sg,共102兲 where Rij r=具 ␳ ¯ u ¯ i ⬘u ¯ j ⬘典/具 ␳ ¯ 典, 共103兲 Rij sg =具qu,ij典/具 ␳ ¯ 典, the subgrid contribution being estimated from the same dynamic mixed model used in the LES. The previous relation between the total, resolved, and subgrid stresses holds if the filter size ⌬fis small enough to allow the approximation 具 ␾ 典⯝具 ␾ 典. Figure 14 shows the temporal evolution at the center plane of the total Reynolds stresses, as defined above, nondimensionalized with the velocity scale 冑Ag ␦ ␳ . The horizontal Reynolds stress is defined by Rh=Rxx+Ryy. The temporal evolution shows the transient seen in the mixing depth, with a strong peak in the vertical Reynolds stress Rzz. It is also observed that the horizontal stress has certain delay with respect to the vertical one, the peak occurring at later times, because the input of energy is through Rzz followed by a transfer to Rh. The tendency toward a self-similar stage is also approximately seen beyond the nondimensional time 8.0, where the plateau in the curves confirms the expected scaling with the velocity 冑Ag ␦ ␳ . However, there is still a slight increase in anisotropy. Figure 15 plots the profiles at the times t=10 ␶ and t =15 ␶ using the similarity variable z/ ␦ ␳ . This figure shows a clear asymmetry in the flow, the peak of the Reynolds stresses being slightly displaced toward the light fluid. It is also observed that the profiles do not collapse very well, especially in the spikes front 共negative z兲, similar to the density profiles. The level of anisotropy of the Reynolds stresses measured by the ratio between the vertical and the horizontal fluctuations is Rzz/Rh⯝1.8, as can be observed in Fig. 15. This result is in agreement with the value 1.7 reported by Cook et al.15 and Dimonte et al.,22 and 1.5 reported by Linden et al.16 This anisotropy is larger than that found in free jets or mixing layers,50 where the ratio of the streamwise intensity to the cross-sectional fluctuation is closer to unity, suggesting that a directional volumetric force induces more anisotropy than a directional shear as driving mechanism. The subgrid Reynolds stresses, as obtained from the dynamic mixed model, Eq. 共85兲, are shown in Fig. 14. The contribution of the Smagorinsky part to the Reynolds stresses is negligible, ⬇1% of the scale similarity part. Their temporal evolution is similar to that of the resolved counterparts. The subgrid kinetic energy at the end of the simulation is about 7% of the total amount, decreasing with time. It is interesting to see how the latter ratio Ksg/Kvaries with L␧/⌬f, where the integral dissipation length is defined by L␧=K3/2/␧sg,␧sg being the mean subgrid-scale dissipation. This dissipation scale achieves an approximate constant value of L␧=0.3 ␦ ␳ beyond t=10 ␶ . Figure 16 shows the expected decay, as the energy is continuously displaced toward the larger resolved sales. Bagget et al.51 reports smaller subgrid contribution for smaller ratios of L␧/⌬f⯝10 in a channel flow, but this is consistent with the fact that their spectral cut-off filter yields more energy in the resolved part of the spectrum. Besides, some dependence on the specifics of the flow might be expected. FIG. 14. Temporal evolution at the center plane of the total Reynolds stresses: 共—兲vertical, Rzz;共---兲horizontal, Rh=Rxx+Ryy. The bottom curves denote, with the same line pattern, the corresponding subgrid parts. FIG. 15. Total Reynolds stresses profiles using self-similar variables: 共—兲 t=10冑 ␦ ␳ ,0/共Ag兲;共---兲t=15冑 ␦ ␳ ,0/共Ag兲. FIG. 16. Variation of the percentage of subgrid energy as the filter size ⌬f decreases in comparison with the dissipation length scale, L␧=K3/2/␧sg. 076101-17 Large-eddy simulation of Rayleigh-Taylor Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp F. Turbulent kinetic energy budget The transport equation for the resolved turbulent kinetic energy per unit mass, Kr, reads 具 ␳ ¯ 典 ⳵ Kr ⳵ t+具 ␳ ¯ 典具w ¯ 典 ⳵ Kr ⳵ z+ ⳵ Tz r ⳵ z=⌸r+具 ␳ ¯ 典共Pr−␧sg兲+Br +⌺r,共104兲 where Tz r=具 ␳ ¯ w ¯ ⬘u ¯ i ⬘2典/2 + 具p ¯ ⬘w ¯ ⬘典+具w ¯ ⬘qu,zz ⬘典, ⌸r=具p ¯ ⬘u ¯ i,i ⬘典, Pr=−Rzz具w ¯ 典,z, 共105兲 ␧sg =−具qu,ij ⬘u ¯ i,j ⬘典/具 ␳ ¯ 典, Br=具 ␳ ¯ ⬘共w ¯ 兲R ⬘典具p ¯ 典,z/具 ␳ ¯ 典, ⌺r=具 ␳ ¯ ⬘共w ¯ 兲R ⬘典具qu,zz典,z/具 ␳ ¯ 典. In these expressions, 共w ¯ 兲R ⬘denotes the Reynolds fluctuation of the vertical velocity and Bis the buoyancy term. This equation is obtained from the first two transport equations in the system, Eq. 共81兲. In comparison with the transport equation for the total turbulent kinetic energy, it does not contain any of the terms due to molecular diffusion, and it includes instead the subgrid-scale part in the transport Tz rand the mean flux ⌺r, in addition to the mean subgrid-scale dissipation ␧sg. Figure 17 shows the temporal evolution at the center plane of the most significant terms in Eq. 共104兲. The subgrid part of the transport term is approximately constant beyond t=10 ␶ while the resolved part increases so that the ratio follows a trend similar to Ksg/K, Fig. 16, being about 7% at the end of the simulation. It is clear that the main input of energy in this flow is the buoyancy term Br. Note that 具p ¯ 典,z ⯝−具 ␳ ¯ 典gis negative, and so is the mass flux 具 ␳ ¯ ⬘共w ¯ 兲R ⬘典because heavy fluid pockets in lighter ambient fluid drop while light pockets in heavy fluid rise. Therefore, the buoyancy term is positive and represents turbulent potential energy being transferred to turbulent kinetic energy 共in other cases, e.g., stable stratification, the transfer could be in the opposite direction兲. Production and convection play a small role because the mean velocity and its gradient are small. The pressure-dilatation term is negligible, which further indicates the quasi-incompressible character of the flow. The energy is distributed in space by the transport term Tz,z r, moving energy from the center of the mixing zone toward the upper and lower fronts, as shown by Fig. 18, where the profiles across the flow are plotted. A similar result was obtained from the spectral analysis performed by Cook et al.14 The energy is finally transferred toward the subfilter scales through the subgrid-scale dissipation ␧sg. The dissipation caused by the numerical dealiasing filter described in Sec. III is 5.5% of ␧sg at t=10 ␶ and 5.0% at t=15 ␶ . The input of energy by the buoyancy term goes directly into the Reynolds stress in the inhomogeneous direction, Rzz, and it is transferred by the pressure-strain correlation, ⌸ij r=具p ¯ ⬘共u ¯ i,j ⬘ +u ¯ j,i ⬘兲典, toward the turbulent fluctuations in the homogeneous plane. It is worth noticing as well that the subgrid-scale model shows backscatter 共owing to the scale similar part兲during the initial time, indicated by the positive values of −␧sg for t ⬍3 ␶ . The curves are nondimensionalized with the quantity 冑共Ag兲3 ␦ ␳ , formed by the length scale ␦ ␳ and the velocity scale 冑Ag ␦ ␳ . It is observed in Figs. 17 and 18 that the temporal evolution of the the buoyancy flux and the transport terms do not reach a plateau after the transient, though the subgrid-scale dissipation clearly does. This is not surprising since different statistical quantities achieve self-similarity at different times. Longer simulations would be interesting to ascertain this behavior. Lastly, the integral of Eq. 共104兲along the vertical direction implies that the input of energy through the buoyancy flux is split between the gain of turbulent kinetic energy and the dissipation. Figure 19 shows the ratio dissipation buoyancy production = 冕 具 ␳ ¯ 典␧sgdz 冕 Brdz .共106兲 Instead of the integral of Br, the loss of potential energy based on the mean density profile is used in the literature. FIG. 17. Temporal evolution at the center plane of the budget of the resolved turbulent kinetic energy equation per unit mass, nondimensionalized by 冑共Ag兲3 ␦ ␳ :共—兲buoyancy production, Br/具 ␳ ¯ 典;共---兲subgrid dissipation, −␧sg;共¯¯¯兲transport, −Tj,j r/具 ␳ ¯ 典;共-·-兲pressure-dilatation, ⌸r/具 ␳ ¯ 典. FIG. 18. Budget of the resolved turbulent kinetic energy equation, nondimensionalized by 冑共Ag3兲 ␦ ␳ :共—兲t=10冑 ␦ ␳ ,0/共Ag兲;共---兲t=15冑 ␦ ␳ ,0/共Ag兲. 076101-18 Mellado, Sarkar, and Zhou Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp Both are equivalent when the domain of integration is closed, as obtained from mass conservation, but the present case allows outflow at the top and the bottom and the amount of potential energy released that goes into the turbulent motion is given by B, as indicated by Eq. 共104兲. Figure 19 shows a tendency toward an asymptotic value which, however, is not achieved at the end of the simulation. The observed value, 0.42, is slightly smaller than the value of 0.5 found in experiments11,16 and a recent numerical simulation.15 VI. CONCLUSIONS Compressibility effects in Rayleigh-Taylor turbulence with miscible fluids in an unbounded domain have been studied using analysis and LES. Three configurations are considered in the theoretical analysis. The first one is a two-layer system formed by a step-like distribution of the ratio between molecular weight and temperature. The density decreases exponentially with increasing height in each layer and each layer is buoyancy-stable. The second configuration has a density jump between two layers each being buoyancyneutral. The third configuration is defined by a step-like profile of the density itself, while the pressure decreases linearly with height in each layer. Each layer is buoyancy-unstable. It has been shown analytically that the turbulent Mach number is bounded from above, in the first two cases independently of the density jump at the interface and, for the last case, for moderate Atwood numbers 共A艋0.5兲. The reason is that the initial thermodynamic state of the system determines the amount of potential energy per unit mass involved in the turbulent mixing stage, and thus the level of turbulent fluctuations that is achievable is linked to the characteristic speed of sound such that the turbulent Mach number is limited. In the particular case considered here of an ideal gas, this bound on the turbulent Mach number is derived to be Mt,max ⯝0.25 in the constant density configuration. However, this initial configuration in a compressible case is buoyancyunstable and may be difficult to set up. The bound is larger in the buoyancy-stable and buoyancy-neutral systems, for which Mt,max ⯝0.6, but the result is independent of the density ratio. It has to be noted that these values are conservative because the amount of potential energy that is dissipated 共which is up to 50% in the incompressible case兲is retained in the estimate of the turbulent kinetic energy. In all situations, Mtis small enough so that compressibility effects may be relatively small. LES performed with a particular density jump at the interface of 3:1 for the stable and unstable configurations considered in the analysis indeed confirm that Mt does not exceed the analytical bounds. The compressibility effects have been studied in the LES decomposing the density fluctuation into the entropic part, due to variation of composition, and the acoustic part, due to intrinsic compressibility. This latter is found to be less than 10% of the total density fluctuation, indicating that the intrinsic compressibility effects are indeed small. Consequently, key features such as the quadratic time evolution of the mixing depth, the anisotropy of the Reynolds stresses, and the value of the mixing parameters compare well with those observed in the incompressible cases reported in the literature. The LES is performed using a dynamic mixed model. The general evolution of the flow has been studied in the configuration with constant density. Mean profiles and mixing parameters are in good agreement with the available incompressible data to the extent that they can be compared. The structural changes in the flow as it evolves from the initial ordered finger-like structure to the disordered turbulent stage are manifest in practically all the discussed quantities. Peaks in density and velocity fluctuations, along with a minimum in mixing that corresponds to a more segregated state, are observed at times t⯝2.5冑 ␦ ␳ ,0/共Ag兲, where ␦ ␳ ,0 is the initial thickness of the mixing region. Subsequently, the vorticity field, in the form of rings around the density fingers, starts to break the ordered density structures, mixing increases, and the subgrid model activates to provide the required dissipation. The Reynolds stresses and the budget of the turbulent kinetic energy have been fully described. The level of anisotropy in the Reynolds stresses, measured as the ratio between the vertical to the horizontal fluctuations, is ⬇1.8. This is higher than that for shear driven flows, where the ratio of streamwise to cross-sectional fluctuations is closer to or smaller than unity. The vertical fluctuations gain energy from the available potential energy, and then there is a transfer to the horizontal fluctuations by the pressure-strain terms. This phenomenon peaks at the center 共approximately兲of the mixing width. The energy is then transported spatially toward the upper and the lower fronts of the turbulent core. Finally, the subgrid-scale dissipation transfers the energy toward the subfilter scales. Although some statistics, notably the thickness of the mixing region, show signs of self-similarity, others continue to slowly evolve in time; the size of the problem, the initial thickness being ⬇4% of the final value, is still too small for complete self-similarity. ACKNOWLEDGMENTS Partial support for J.P.M. was provided by Lawrence Livermore National Laboratory through the Student Employee Graduate Research Fellowship Program. This work FIG. 19. Temporal evolution of the ratio of depth-integrated values of subgrid-scale dissipation and buoyancy production. 076101-19 Large-eddy simulation of Rayleigh-Taylor Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp was supported in part by a grant of HPC time from the Naval Oceanographic Office Department of Defense Major Shared Resource Center and by the National Partnership for Advanced Computational Infrastructure. We thank an anonymous reviewer for constructive comments. 1Lord Rayleigh, “Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density,” Proc. London Math. Soc. 14, 170 共1883兲. 2G. I. Taylor, “The instability of liquid surfaces when accelerated in in a direction perpendicular to their planes,” Proc. R. Soc. London, Ser. A 201, 192 共1950兲. 3R. E. Duff, F. H. Harlow, and C. W. Hirt, “Effects of diffusion on interface instability between gases,” Phys. Fluids 5, 417 共1962兲. 4D. Livescu, “Compressibility effects on the Rayleigh-Taylor instability growth between inmiscible fluid,” Phys. Fluids 16,118共2004兲. 5K. I. Read, “Experimental investigation of turbulent mixing by RayleighTaylor instability,” Physica D 12,45共1984兲. 6D. L. Youngs, “Numerical simulation of turbulent mixing by RayleighTaylor instability,” Physica D 12,32共1984兲. 7D. H. Sharp, “An overview of Rayleigh-Taylor instability,” Physica D 12, 3共1984兲. 8D. L. Youngs, “Modelling turbulent mixing by Rayleigh-Taylor instability,” Physica D 37,270共1989兲. 9Y. Zhou, “A scaling analysis of turbulent flows driven by Rayleigh-Taylor and Richtmyer-Meshkov instabilities,” Phys. Fluids 13, 538 共2001兲. 10P. F. Linden and J. M. Redondo, “Molecular mixing in Rayleigh-Taylor instability. Part I: Global mixing,” Phys. Fluids A 31269 共1991兲. 11P. Ramaprabhu and M. J. Andrews, “Experimental investigation of Rayleigh-Taylor mixing at small Atwood numbers,” J. Fluid Mech. 502, 233 共2004兲. 12D. L. Youngs, “Three-dimensional numerical simulation of turbulent mixing by Rayleigh-Taylor instability,” Phys. Fluids A 3, 1312 共1991兲. 13A. W. Cook and P. E. Dimotakis, “Transition stages of Rayleigh-Taylor instability between miscible fluids,” J. Fluid Mech. 443,69共2001兲. 14A. W. Cook and Y. Zhou, “Energy transfer in Rayleigh-Taylor instability,” Phys. Rev. E 66, 026312 共2002兲. 15A. W. Cook, W. Cabot, and P. L. Miller, “The mixing transition in rayleigh-taylor instability,” J. Fluid Mech. 511, 333 共2004兲. 16P. F. Linden, J. M. Redondo, and D. L. Youngs, “Molecular mixing in Rayleigh-Taylor instability,” J. Fluid Mech. 265,97共1994兲. 17S. B. Dalziel, P. F. Linden, and D. L. Youngs, “Self-similarity and internal structure of turbulence induced by Rayleigh-Taylor instability,” J. Fluid Mech. 399,1共1999兲. 18C. Meneveau and J. Katz, “Scale-invariance and turbulence models for large-eddy simulations,” Annu. Rev. Fluid Mech. 32,1共2000兲. 19R. S. Rogallo and P. Moin, “Numerical simulation of turbulent flows,” Annu. Rev. Fluid Mech. 16,99共1984兲. 20U. Piomelli and J. R. Chasnov, in Large-Eddy Simulations: Theory and Applications, Turbulence and Transition Modelling, edited by M. Hallbäck, D. S. Henningson, A. V. Johansson and P. H. Alfredsson 共Kluwer, Dordrecht, 1996兲, pp. 269–336. 21M. Lesieur and O. Métais, “New trends if large-eddy simulation of turbulence,” Annu. Rev. Fluid Mech. 28,45共1996兲. 22G. 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. Young, and M. Zingale, “A comparative study of the turbulent Rayleigh-Taylor instability using high-resolution threedimensional numerical simulations: The Alpha-group collaboration,” Phys. Fluids 16, 1668 共2004兲. 23P. Chassaing, R. A. Antonia, F. Anselmet, L. Joly, and S. Sarkar, Variable Density Fluid Turbulence 共Kluwer Academic, Dordrecht, 2002兲. 24S. Sarkar, G. Erlebacher, M. Y. Hussaini, and H. O. Kreiss, “The analysis and modeling of dilatational terms in compressible turbulence,” J. Fluid Mech. 227,473共1991兲. 25S. Sarkar, “The stabilizing effect of compressibility in turbulent shear flow,” J. Fluid Mech. 282,163共1995兲. 26A. W. Vreman, N. D. Sandham, and K. H. Luo, “Compressible mixing layer growth rate and turbulence characteristics,” J. Fluid Mech. 320,235 共1996兲. 27S. Sarkar, “On density and pressure fluctuations in uniformly sheared compressible flow,” Proceedings, IUTAM Symposium on Variable Density Low-Speed Flows, Marseille, edited by L. Fulachier, J. L. Lumley, and F. Anselmet 共Kluwer Academic, Dordrecht, 1996兲. 28J. B. Freund, S. K. Lele, and P. Moin, “Compressibility effects in a turbulent annular mixing layer. 1. Turbulence and growth rate,” J. Fluid Mech. 421, 229 共2000兲. 29C. Pantano and S. Sarkar, “A study of compressibility effects in the highspeed, turbulent shear layer using direct simulation,” J. Fluid Mech., 451, 329 共2002兲. 30A. E. Gill, Atmosphere-Ocean Dynamics 共Academic, New York, 1982兲. 31J. Pedloski, Geophysical Fluid Dynamics 共Springer, Berlin, 1987兲. 32G. K. Batchelor, An Introduction to Fluid Dynamics 共Cambridge University Press, Cambridge, 1967兲. 33Y. Chen, Y. Deng, J. Glimm, G. Li, and Q. Zang, “A renormalization group scaling analysis for compressible two-phase flow,” Phys. Fluids A 5, 2929 共1993兲. 34B. Vreman, B. Geurts, and H. Kuerten, “Subgrid-modelling in LES of compressible flow,” Appl. Sci. Res. 54,191共1995兲. 35B. Vreman, B. Geurts, and H. Kuerten, “Large-eddy simulation of the turbulent mixing layer,” J. Fluid Mech. 339, 357 共1997兲. 36A. Yoshizawa, “Statistical theory for compressible turbulent shear flows, with application to subgrid modeling,” Phys. Fluids 29, 2152 共1986兲. 37C. G. Speziale, G. Erlebacher, T. A. Zang, and M. Y. Hussaini, “The subgrid-scale modeling of compressible turbulence,” Phys. Fluids 31,940 共1988兲. 38P. Moin, K. Squires, W. Cabot, and S. Lee, “A dynamic subgrid-scale model for compressible turbulence and scalar transport,” Phys. Fluids A 3, 2746 共1991兲. 39G. Erlebacher, M. Y. Hussaini, C. G. Speziale, and T. A. Zang, “Toward the large-eddy simulation of compressible flows,” J. Fluid Mech. 238,155 共1992兲. 40T. A. Zang, R. B. Dahlburg, and J. P. Dahlburg, “Direct and large-eddy simulation of three-dimensional compressible Navier-Stokes turbulence,” Phys. Fluids A 4, 127 共1992兲. 41S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” J. Comput. Phys. 103,16共1992兲. 42F. K. Chow and P. Moin, “A further study of numerical errors in largeeddy simulations,” J. Comput. Phys. 18, 366 共2003兲. 43J. H. Williamson, “Low-storage Runge-Kutta schemes,” J. Comput. Phys. 35,48共1980兲. 44K. W. Thompson, “Time-dependent boundary conditions for hyperbolic systems, II,” J. Comput. Phys. 89, 439 共1990兲. 45A. G. Kravchenko and P. Moin, “On the effect of numerical errors in large-eddy simulations of turbulent flows,” J. Comput. Phys. 131,310 共1997兲. 46Y. Morinishi, T. S. Lund, O. V. Vasilyev, and P. Moin, “Fully conservative higher order finite difference schemes for incompressible flows,” J. Comput. Phys. 143,90共1998兲. 47G. A. Blaisdell, E. T. Spyropoulos, and J. H. Qin, “The effect of the formulation of nonlinear terms on aliasing errors in spectral methods,” Appl. Numer. Math. 21, 207 共1996兲. 48C. Le Ribault, S. Sarkar, and S. A. Stanley, “Large eddy simulation of a plane jet,” Phys. Fluids 11, 3069 共1999兲. 49C. Le Ribault, S. Sarkar, and S. A. Stanley, “Large eddy simulation of evolution of a passive scalar in plane jet,” AIAA J. 39 1509 共2001兲. 50S. B. Pope, Turbulent Flows 共Cambridge University Press, Cambridge, 2000兲. 51J. S. Bagget, J. Jiménez, and A. G. Kravchenko, “Resolution requirements in large-eddy simulations of shear flows,” CTR Annual Research Briefs 共Stanford University, Stanford, 1997兲. 076101-20 Mellado, Sarkar, and Zhou Phys. Fluids 17, 076101 共2005兲 Downloaded 27 Jul 2005 to 129.187.68.124. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp