# Numerical Simulation of Single-Mode 3D Rayleigh-Taylor Instability

Authors: Dongyung Kim
Journal: Modelling and Simulation in Engineering
 (2022)

---

## Abstract

Rayleigh-Taylor instability (RTI) is analyzed theoretically by Taylor, and 2-dimensional experimental results are obtained by Lewis in 1950. Over the 72 years, several experiments and theories are developed with the shock-driven Ritchmyer-Meshkov instability (RMI) and the shear-driven Kelvin-Helmholtz instability (KHI). Here, we emphasize the single-mode Rayleigh-Taylor instability (RTI) mixing simulation with a surface area in 3 dimensions. The simulation uses concentration equations and nonzero transport. We observed chaotic interface behavior even for this single-mode simulation, in the sense that the interface appears to have an area proportional toΔx−1, with respect to its mesh (non)convergence (i.e., rate of divergence) properties.


---

## Body

## 1. Introduction

Many important industrial problems involve flows with multiple constitutive components. Due to inherent nonlinearities and the complexity of dealing with unknown moving interfaces, multiphase flows are challenging. There are many ways to model moving interfaces. The two main approaches to simulating multiphase and multicomponent flows are interface tracking and interface capturing. In interface tracking methods (front-tracking [1], immersed interface [2], immersed boundary [3, 4]), Prometheus with PPM(the piecewise-parabolic method) [5], and CLAWPACK [6]), Lagrangian particles are used to track the interfaces. In interface capturing methods such as level-set [7, 8] and phase-field methods [9], the interface is implicitly captured by a contour of a particular scalar function. Another types of methods are Clawpack and Prometheus Method (PPM). Reference will be inserted. The numerical method we will take is Front-tracking method (FTM), which is originally developed by Glimm et al. [1] and Immersed boundary method (IBM), which is originally developed by Peskin [4]. The IBM was started to be applied to two-phase fluid flows [10–12], and the FTM was started to be applied to two-phase fluid flows. The motion of the fluid with the interface tracking method is influenced by the force generated by the interface, and the interface moves at the local fluid velocity. The strength of this method is accurate and robust with time-dependent geometry of the interface. See Figure 1 for 2-dimensional incompressible and immiscible two fluid mixing cases. Here, ρ is the variable density, and μ is the variable viscosity. Xs,t is the interface of two fluid mixing in 2 dimensions.Figure 1
Front-tracking and Immersed boundary curveΓ on a rectangular domain Ω.We extend the problems to a 3-dimensional space with compressible two fluid mixing cases. This problem is one of Rayleigh-Taylor instabilities (RTI) [13] with the density differences. RTI is applied to resolve supersonic ramjet [14] and scramjet or inertial confinement fusion (ICF) problems [15]. We consider the Navier-Stokes equations with transport for a mixture of two compressible species:
(1)∂ρ∂t+∇·ρv=0,(2)∂ρv∂t+∇·ρvv+pδij=∇·d,(3)∂E∂t+∇·E+pv=∇·κ∇T+∇·d·v,(4)∂ρψ∂t+∇·ρψv=∇·ρμ∇ψ.The following are the dependent variables:ρ is the total mass density, v is the velocity, p is the pressure, E is the total energy, ψ is the species mass fraction, κ is the coefficient of heat conductivity, and μ is the molecular mass diffusion.Equation (1) is the continuum equation, equation (2) is the momentum equation, equation (3) is the energy equation [16], and equation (4) is the mass fraction equation. In equation (4), we added and developed the volume fraction formulas to preserve the mass in the interface because it is not conservative in the interface, grid cell crossings and it is conservative without interface. After we add and develop the volume fraction formulas, it almost becomes conservative in the interface and grid cell crossing [17].The equation of state (EOS) is defined for each of the species as a gamma law gas, and according to [18], the mixture EOS is a gamma law gas also. d is the viscous stress d=2νS−1/3trSI, and Sij is the strain rate tensor:
(5)Sij=12∂vi∂xj+∂vj∂xi.In the reference book of Williams [19], there is a more detailed approximation theory of multifluid viscosity. We implement Neumann boundary conditions with reflecting boundary state at the walls for the x-velocity component. We consider compressible flows coupling the concentration equation actively into the flow dynamics. The mixing problems which we study are driven by acceleration-driven forces. The classical Rayleigh-Taylor instability with mixing regime is defined by steady acceleration of a density discontinuities. The Atwood number A=ρ2−ρ1/ρ2+ρ1, with A>0, is an important parameter to measure the effective buoyancy and thus the acceleration of the flow. See the overview in [20]. In this single-mode Rayleigh-Taylor 3D simulation, we update the numerical tool FronTier which is based on a front tracking algorithm with a high-quality treatment of a fluid interface with volume fraction formulas. This code is validated [21]. In detail, we developed a local grid-based method, one of the front tracking methods with the volume fraction formulas. A front grid (a codimension 1 grid) specifies a fluid discontinuity location. Through a regular rectangular grid, the front moves freely. In the interface and grid cells, the volume fraction formula is very effective for conserving the mass fractions. At the front, Riemann solutions, which are constructed in a normal direction, provide the physics-based dynamics to move the front at each time step.


## 2. Multiscale Modeling of Single-Mode 3D RT Instability

We simulate 3D single-mode Rayleigh-Taylor with nonzero transport and obtain the surface area with averaged volume fraction formulas. It was already observed [22] without averaged volume fraction formulas; the interface for a related but distinct 2D Richtmyer-Meshkov instability has a length proportional to Δx−1, with respect to its mesh (non)convergence (i.e., rate of divergence) properties.The mechanism of the interfacial instabilities is the following. The amplitude of mode is 0.1, the phase of mode is 270 degrees, and the boundary type of perturbation in direction is periodic. We show in Figures2–5 several time steps in the evolving flow. Here, we used transport coefficients for water (Navier-Stokes viscosity: 0.00085105 cm2/ms, Navier-Stokes mass diffusion: 0.00016366 cm2/ms, and Navier-Stokes thermal conductivity: 0.00112 cm2/ms).Figure 2
Plot of the interface between the heavy (above) and light (below) fluids at initial time for Rayleigh-Taylor fluid instability. For this study, density ratio is 2 : 1, the Atwood number is 1/3, the peak to peak amplitude of initial disturbance is 0.06, grid size is20×20×200, and computational domain is 1cm×1cm×10cm.Figure 3
Plot of early time for Rayleigh-Taylor fluid instability. Flow and grid parameters as in Figure2.Figure 4
Plot of middle time for Rayleigh-Taylor fluid instability. Flow and grid parameters as in Figure2.Figure 5
Plot of late time for Rayleigh-Taylor fluid instability. Flow and grid parameters as in Figure2.In the interface surface area in Figure6, we observed that the surface area is divergent with time (expressed unphysical units) in three mesh refinements. We convert the surface area to mesh units following the formula
(6)interfacesurfacemixingzonevolume=physicalsurfacearea/ΔX∗ΔYhmax−hmin∗domainX∗domainY/ΔX∗ΔY∗ΔZ=physicalsurfacearea∗ΔZhmax−hmin∗domainX∗domainY.Figure 6
Plot of the interface surface area vs. time using physical units. Results for three mesh levels are displayed for the identical 3D single-mode RT instability.Here,hmax is the maximum height, and hmin is the minimum height in the mixing zone. We observe the chaotic interface behavior in Figure 7. By (6), at late time, the mesh fraction of surface area per mixing zone mesh block ranges from 24% to 14%. We interpreted these members as indicating that the flow is a mixture of chaotic patches and nonchaotic patches.Figure 7
Plot of the interface area divided by the mixing zone volume vs. time. Both the area and the volume are measured in mesh units (thequantityphysicalarea/physicalvolume×Δxisplotted). Results for three mesh levels are displayed for a 3D single-mode RT instability.The flow morphology of a single-mode Rayleigh-Taylor disturbance is far less chaotic than a comparison multimode 2D Richtmyer-Meshkov flow, but there is still vortex shedding from the mushroom caps. The flow is more or less chaotic in these vortex shedding regions at late time. So, we study the chaotic flow near the bottom of the mushroom caps and the top of the spike tips in more detail [23].The mesh level surface fraction is time independent. In the later time, the mesh level surface fraction is about 20% relative to the mixing zone itself. The divergent nature prototypical error analysis of the interface without transport physics is proposed as a formula form [23]:
(7)Error=C1×Δx×interfacearea=C1C2.Ci are O1 constants independent of Δx. C1 is related to numerical mass diffusion and might be taken as 3.0 for a typical numerical algorithm [24]. In the single-mode 3D problem and for grid levels considered here, C2 has a value in the range 0.14 to 0.24.To study the convergence and mesh refinement, we defineRe=VL/<ν>, ν is the kinematic viscosity, and <·> is an ensemble average. V is the turbulent fluid velocity, V=<δV2>. Here, δVz=Vz−<Vz> is the fluctuating part of the velocity, and L is the mixing zone height (hmax−hmin). In the simulation, we used a constant dynamic viscosity. We also use the mesh Reynolds number Remesh=VΔx/<ν> and the Schmidt number Sc=ν/μ as a new dimensionless parameter, where μ is the coefficient of molecular mass diffusion. The Kolmogorov length scale λK=ν3/ε1/4 or the viscous inner scale, approximately 50λK, is a measure of the lengths at which viscosity plays a role, and this is related to the level of mesh refinement needed for a numerical simulation. For comparison, Δx for the finest grid is comparable to λK and well below the inner viscous scale, indicating that the calculation is in close to fully resolved simulation. We need to do more mesh refinement to get fully resolved simulation. Here, ε is the dissipation rate:
(8)ε=ν2Sij22,where Sij is (5). See Table 1 relating Re and λK to Remesh.Table 1
Table relatingRe and λK to Remesh at t=11.0 for several mesh levels. Kolmogorov scale λKmesh=λK/Δx. Schmidt number Sc=5.2.
Grid sizeReRemeshλKmesh10×10×100731.5032.511.2220×20×200937.0421.091.7140×40×4001212.3612.552.80


## 3. Conclusions

In this paper, we investigate the simulation of a single-mode Rayleigh-Taylor instability with the volume fraction formulas in 3D. By the implementation of surface area function, we can precisely observe the chaotic behavior of the interface flow. The flow morphology is far less chaotic than a comparison 2D multimode Richtmyer-Meshkov flow. There is still vortex shedding from the mushroom caps, at which locations of the flow are more or less chaotic. So, we can study the chaotic flow of the bottom of the mushroom caps and the top of the spike tips in more detail. Here, we studied the Kolmogorov length scale which is relative to the level of mesh refinement for numerical simulation. Finally, in the future work, we will apply Prandtl and Batchelor scales which are need to resolve more fully for numerical simulation.



---

*Source: 6544820-2022-11-09.xml*