Large eddy simulation of high atwood number rayleigh-taylor mixing

. Large eddy simulation of Rayleigh-Taylor instability at high Atwood numbers is performed using recently developed, kinetic energy-conserving, non-dissipative, fully-implicit, ﬁnite volume algorithm. The algorithm does not rely on the Boussinesq assumption. It also allows density and viscosity to vary. No interface capturing mechanism is requried. Because of its advanced features, unlike the pure incompressible ones, it does not su ↵ er from the loss of physical accuracy at high Atwood numbers. Many diagnostics including local mole fractions, bubble and spike growth rates, mixing e ﬃ ciencies, Taylor micro-scales, Reynolds stresses and their anisotropies are computed to analyze the high Atwood number e ↵ ects. The density ratio dependence for the ratio of spike to bubble heights is also studied. Results show that higher Atwood numbers are characterized by increasing ratio of spike to bubble growth rates, higher speeds of bubble and especially spike fronts, faster development in instability, similarity in late time mixing values, and mixing asymmetry.


Introduction
Rayleigh-Taylor Instability (RTI) occurs when a heavy fluid of density ⇢ h on top is supported against the gravity, g, by a light fluid of density ⇢ l on bottom [1][2][3]. When it is subject to perturbations, the fluids which are initially in hydrostatic equilibrium start to interpenetrate each other due to the baroclinic vorticity generated by the opposite density and pressure gradients. It can be observed in flows such as type Ia supernovae, Inertial Confinement Fusion (ICF), cavitation bubbles, oceanic and atmospheric currents and many other natural and engineering flows. In terms of flow modeling, RTI serves as very challenging test case for algorithms, as it includes many aspects such as body force treatment, turbulent mixing, di↵usion and interface capturing. RTI is characterized by the Atwood number, A, which is given as A = ⇢ h −⇢ l ⇢ h +⇢ l . Penetration of the light fluid into the heavy one as bubbles and penetration of the heavy fluid into the light one as spikes can be modeled in terms of the penetration lengths as h b = ↵ b Agt 2 and h s = −↵ s Agt 2 . ↵ b and ↵ s are the growth rates for bubble and spike, respectively. The mixing zone and its rate of change in time are then found as, h = h b − h s = ↵Agt 2 andḣ = 2↵Agt, where ↵ = ↵ b + ↵ s is the total growth rate factor.
In some flows in engineering and astrophysics, such as ICF and supernovae collapse, the Atwood number can reach very high values. Although RTI at low to intermediate Atwood numbers have been investigated extensively, experimental and numerical studies on high and very high Atwood number three-dimensional multi-mode RTI are rarely available in the literature. It was first reported by Youngs [4]. Suitability of the numerical method used in ⇤ e-mail: ilyas.yilmaz@bilgi.edu.tr his paper was discussed by Burton [5]. Additionally, only the overall characteristics of growth and mixing were studied. Dimonte and Schneider [6] studied RTI up to A = 0.96 using Linear Electric Motor. They provided scaling laws between spike to bubble penetration ratio and heavy to light density ratio, and observed greater asymmetry at high Atwood numbers. They also found that the bubble sizes increase while the spikes become narrower. Large Eddy Simulation (LES) of RTI up to A = 0.96 were also presented by Burton [5]. It was one of the most comprehensive works on this subject and provided many diagnostics and measures. A series of high-resolution simulations with an aim to provide benchmark solutions for use in engineering models were also performed by Youngs [7]. Recently, another numerical study was reported by Shmony et. al. [8] to examine the density ratio e↵ects on the mixing process and the e↵ective Atwood number. However, the literature on this area is limited and not complete, and it is worth to enhance it by adding more results obtained from advanced simulation methods.
The aim of this study is twofold. Firstly, to show that the present LES algorithm is able to capture the evolution and characteristics of the RTI at moderate Atwood numbers in a physically more correct way and then to investigate the e↵ects of increasing Atwood numbers (i.e., higher density ratios) on the development of the RTI and contribute to the literature.
second-order, pressure-correction type algorithm based on an iterative predictor-corrector approach to update the flow variables [10]. A pressure relaxation procedure was also added to remove the possible oscillations and speed-up the solution [11]. The well-known Wall Adapting Local Eddy (WALE) viscosity model was employed to model the subgrid-scale viscosity [12].

Solver
An in-house, fully-implicit, fully-parallel, structured, single-block, finite-volume LES solver, written using PETSc [13] with Fortran syntax, is used. It was previously validated for various transitional and turbulent flows. See for example [9,11,14,15]. Incomplete LU preconditioned-GMRES with no fill-in, ILU(0), solves the linear systems stemming from the implicit discretization. Domain decomposition and load balancing are determined at the beginning by PETSc. See [9] for other details.

Problem setup
A series of runs are performed in order to properly analyze the Atwood number (i.e., density ratio) e↵ects. In each run, the flow is initialized with a specified Atwood number, ranging from 0.5 to 0.9, which corresponds to the density ratios from 3 to 19 respectively. There are five Atwood numbers studied in total.
The fluids are initially at rest with constant downward gravitational acceleration and in hydrostatic equilibrium p = p 0 − ⇢gy. p 0 is the reference pressure. The reference density ⇢ 0 is the light fluid density. The height of the domain is taken as the reference length l 0 . The free-fall velocity is used as the reference velocity u 0 = p gl 0 . This leads to the reference time scale ⌧ = q l 0 Ag . The Eckert number can be approximated as Ec = (γ−1)M 2 r . The reference Mach number, M r , is also defined as, and is set to 0.1. p 0 is computed using M r . The reference dynamic viscosity µ 0 here can be obtained using the following grid-dependent relation, where ! is the coefficient taken as 0.21 for A = 0.5 [18], and is increased up to 3.6 for A = 0.9. These choices lead to the scaling(working) Reynolds number, Re = ⇢ 0 u 0 l 0 µ 0 , ranging from 1.7 ⇥ 10 4 to 8 ⇥ 10 3 . The Froude number, Fr, which is the ratio of the reference velocity to the free-fall velocity is equal to unity. The initial non-dimensional temperature field, T , is found from the non-dimensional equation of state. Then, the non-dimensional dynamic viscosity, µ, can be computed using µ = T n . The value of n is set to 0.71 in the simulations. It is also possible to obtain the nondimensional fields for T and µ at initial by dividing their reference values. However it is not preferred. Since the approach followed here allows not only the density but also the temperature and the dynamic viscosity to vary during the simulation, and it is expected to give physically more accurate results, especially at high Atwood numbers.
The domain size is H ⇥ 3H ⇥ H in the stream-wise, the normal-wise and the span-wise directions respectively, with a resolution of 64 ⇥ 192 ⇥ 64. H is set to 2⇡ . The interface is located at y=0.
The stream-wise and the span-wise boundary conditions are periodic. The no-slip wall boundary condition is applied to the velocity at top and bottom. The homogeneous Neumann boundary condition for the thermodynamic variables are used at walls.
Usually, the interface is perturbed to trigger the instability. However, this may require an additional grid smoothing, if the amplitude of the perturbation is smaller than the grid resolution. In addition, a small initial perturbations for velocity may be necessary for consistency at the inter-facial region [16].
Jun et. al. [17] showed that directly perturbing the velocity field is superior to perturbing the interface position, since no smoothing is required and resolution-independent perturbations are obtained which yield almost identical results. Following the latter way, perturbation is added to the vertical velocity component as v 0 (y) = A 0 Rand (1 + cos(2⇡y/L y )) with the amplitude A 0 is 0.005 and the random number, Rand, between -1 and +1. The perturbation spectrum is a white noise almost down to the grid scale, meaning that at early times high-k modes are dominant. The perturbations are also small amplitude, zoneto-zone perturbations and decrease toward to boundaries. The random number generator implemented in PETSc, which generates purely real random number and creates normally (Gaussian) distributed values within the given range, is used to generate random values.

Results and discussion
As the primary flow diagnostic, the local mole fraction (LMF) field based on heavy fluid is computed as Then, it is averaged over the homogeneous directions to give the mixing zone growth as Using the traditional measures based on LMF which are hχ(x, y, z, t)i xz  1 − ✏ and hχ(x, y, z, t)i xz ≥ ✏ where ✏ is traditionally set to 0.01, the bubble and spike penetration lengths are obtained, respectively. Now, it is easy to compute the total extent of mixing zone, its rate of change in time, the corresponding growth rates and the total growth rate from the equations given in Introduction. The mixing efficieny ✓ which quantifies the amount of mixing can also be derived from LMF via which produces a value between 0 (immiscible) and 1 (complete mixing). Fig.1 shows the evolution of RTI in terms of LMF isosurfaces. Three of the five Atwood numbers are selected to represent. The development of RTI is clearly observed. Initially small structures merge and evolve into large ones. Higher Atwood numbers characterized by higher bubble and spike velocities. As a result, they reach boundaries earlier. Additionally, the morphologies of the structures change, as Atwood number increases. While the bubbles are getting large in diameter, the spikes are becoming narrower. This also agrees well with Dimonte and Schneider [6] and Burton [5].
How the bubble and spike penetration lengths vary in time is shown in Fig.2. It suggests that there two phases of the development of the RTI: an initial linear phase which The corresponding growth rates are plotted in Fig.3. For A = 0.5, ↵ b is found as 0.027 at late times. The reported values are between 0.02 and 0.03 [5,16,18]. The late time total growth rate (not shown here) is also obtained as 0.067 which falls into the range of 0.01-0.07 [5,16]. Unlike the many other numerical methods suffering from excessive dissipation, the present algorithm takes the advantage of non-dissipativeness and produces a higher value close to the experimental data (i.e., the upper limit). No significant e↵ect of the higher Atwood numbers on the bubble growth rate is observed. However, the spike growth rates increase slightly. A late time steadystate value is reached.
In Fig.4, it is seen that mixing decreases rapidly during the di↵usive phase. Then, it starts to rise, beginning from the non-linear phase, and reaches a nearly steady value at very late times. For A = 0.5 case, ✓ is around 0.85 consistently with the values in the literature [5,16,19]. This result indicates that molecular mixing is well-captured by the algorithm even on coarser grids compared to highly resolved simulations. Fig.4 also suggests that ✓ min moves back to earlier times as A increases, and is getting smaller  and smaller. In the non-linear phase, it recovers itself and reaches a plateau around 0.8. Compared to the intermediate Atwood numbers, the recovery period takes longer (i.e., larger time-span).
In the second group of flow diagnostics, the Reynolds stress tensor R i j , the anisotropy tensor b i j , the potential and the kinetic energies are computed using the following equations respectively. k sgs is the turbulent kinetic energy. give an insight about the direction of the turbulent momentum transport. In Fig.5, vertical distribution of the root mean square of the vertical component of the Reynolds stress tensor R 22 , which is a↵ected mostly, is plotted at the end of the each simulation. As the Atwood number increases, R 22 value increases and the approximate symmetry observed at lower Atwood numbers is broken. The peak locations are also getting closer to the bottom wall and move into the spike-dominated region. Considering the larger spike velocities and the penetration lengths at higher Atwood numbers, this result is in consistency with the previous observations presented here. Intensity of the spike interactions with their surroundings is higher than those of the bubbles and this can be considered as the primary mechanism that leads the instability to transition to turbulence and a state of turbulence eventually.
Consistently, the largest variations are observed in b 22 which is given in Fig.6. As the Atwood number increases, the initial fluctuations in di↵usive phase are getting stronger and approach to the upper limit given as 2/3 in the literature. However, at later times (i.e., in the nonlinear phase), b 22 values reach a steady-state value around 0.3 [21] How the ratio of the kinetic to the potential energy evolves in time is shown in Fig.7. After an initial settling period, a rapid increase is observed around t/⌧ = 1.5. Then, at later times, as observed in many other flow diagnostics, an almost steady-state value around 0.4 is reached. This is in consistent with the range found in previous Direct Numerical Simulation (DNS) and LES studies which is 0.4 − 0.5 [5,18,20].

Conclusions
RTI is simulated at high Atwood numbers using an advanced, fully-featured LES algorithm. In addition to the unique features of the algorithm such as nondissipativeness, kinetic energy conserving, fully-implicit discretization, no dependence on Boussinesq assumption low-Mach number handling capability and variable property, to the best knowledge of the author, this is the first time the WALE subgrid-scale model is applied to multimode RTI at high Atwood numbers. It should also be considered that the grid resolution is coarser than many other numerical studies previously reported. The Atwood numbers studied are ranging from intermediate to very high values which are 0.5, 0.6, 0.7, 0.8 and 0.9. The corresponding density ratios are 3, 4, 5 .7, 9 and 19. Such high density ratios may be observed in , for example, ICF and astrophysical flows.
Many flow diagnostics are computed and used to analyze the complex evolution of the instability. Due to the limited avaliability of space, some of the results are not presented here. Again, some of them are computed for the first time at high Atwood numbers, such as the Reynolds stresses, and these may provide some new insights on the e↵ects of high Atwood numbers.
It is shown that the algorithm is capable of repropducing the complex physics and characteristics of the RTI at A = 0.5. The results are successfully compared with the previous experimental and numerical data where available. They also suggest that the evolution of the RTI has multiple stages. Following an early stage di↵usive growth, a non-linear interactions take place and lead the instability to a state of turbulence.
It is observed that higher Atwood numbers are characterized by increasing ratio of spike to bubble growth rates, higher speeds of bubble and especially spike fronts, faster development in instability, large scale anisotropy, higher rates of interactions around spikes, similarity in late time values of mixing and some other measures, change in bubble and spike morphologies, and mixing asymmetry. A scaling relation between h s /h b and ⇢ h /⇢ l is also provided.
Another observation is the e↵ect of solid walls on the development of the instability. This e↵ect may re-duce the penetration lengths (i.e., h s /h b ratios) and mixing, as spikes and bubbles approach the boundaries, as also pointed out in [5].
Since this is an ongoing study, the preliminary results presented here will be refined more by performing additional simulations especially at high Atwood numbers to improve the late time statistics and characteristics of RTI.