Direct numerical simulations of two-dimensional (2D) and three-dimensional (3D), single-mode and multi-mode, incompressible immiscible Rayleigh–Taylor (RT) instabilities are performed using a phase-field approach and high-order finite-difference schemes. Various combinations of Atwood number, Reynolds number, surface tension, and initial perturbation amplitude are investigated. It is found that at high Reynolds numbers, the surface tension, if significant, could prevent the formation of Kelvin–Helmholtz type instabilities within the bubble region. A relationship is proposed for the vertical distance of the bubble and spike vs the Atwood number. The spike and bubble reaccelerate after reaching a temporary plateau due to the reduction of the friction drag as a result of the formation of the spike vortices and also the formation of a momentum jet traveling upward within the bubble region. The interface for a 3D single-mode instability grows exponentially; however, a higher Reynolds number and/or a lower Atwood number could result in a noticeably larger surface area after the initial growth. It is also shown that a 3D multi-mode RT instability initially displays an exponential interface growth rate similar to single-mode RT instabilities. Due to the collapse and merging of individual single-mode instabilities, the interface area for a multi-mode RT instability is strongly dependent to the mesh resolution after the exponential growth rate. However, the ratio of kinetic energy over released potential energy exhibits an almost steady state after the initial exponential growth, with values around 0.4, independently of the mesh resolution.

The Rayleigh–Taylor (RT) instability develops at the perturbed interface of two fluids of different densities subjected to relative acceleration.1–3 It is named after Rayleigh4 and Taylor5 who were the first to study (separately) such flow patterns using classical linear stability theories. The RT instability may occur across a wide range of temporal and spatial scales in a variety of natural and engineering fluid systems with perturbed interfaces.3 For instance, the RT instability may play a key role in the dynamics of fluid flows in astrophysics (flame acceleration in supernovae6), geology (initiation and evolution of polydiapirs7), atmospheric and climate science (formation of the mammatus clouds8), inertial confinement fusion (premature fuel mixing and reduction in the heating efficacy9), and designing and fabricating soft materials (structured materials in thin polymeric liquid films10). Although the principal formation mechanism of the RT instability is buoyancy-induced, various additional factors, such as compressibility, rotation, surface tension, and viscosity, may influence its formation and evolution.3,11–14 From a theoretical point of view, such instabilities can be divided into two categories: (i) single-mode, in which the interface is perturbed by a specific perturbation wave with a known initial amplitude and (ii) multi-mode, in which the interface is perturbed in a way that more than one perturbation wave (typically many) with various amplitudes and growth rates are formed.

The evolution of an RT instability can be divided into four stages: (i) an exponential growth of the initial perturbation (given an initial perturbation amplitude greater than some critical value, normally above 0.1λ,1 where λ is the perturbation wavelength), (ii) a nonlinear growth and the formation of the spike and bubble structures, (iii) the spike rolling up and the formation of mushroom-like patterns due to significant amplification of the nonlinear growths and Helmholtz instability along the side of the spike, and (iv) spike breakup and perturbation of the bubble and consequently turbulent mixing of the two fluids.1,15 The RT instability is typically characterized by the density difference between the two fluids through the Atwood number At defined as

At=ρ̃hρ̃lρ̃h+ρ̃l,
(1)

where ρ̃ is the density and subscripts h and l refer to the heavy and light fluids, respectively.1 It should be noted that the last two stages of the evolution of a RT instability can be significantly influenced by the Atwood number. Specifically, the rolling up of the spike and its breakup might be delayed or even not occur at relatively high Atwood numbers.16–18 The Reynolds number may also play a key role in the evolution and complex flow features of the RT instabilities.11,19 It is defined as Re=ŨL̃ρr̃/μr̃, where Ũ,L̃,ρr̃, and μr̃ denote the reference velocity, length, density, and viscosity, respectively. Please note that in this paper, only variables with a tilde sign (∼) are dimensional.

There are many theoretical, computational, and experimental studies of RT instabilities available in the literature.3,20–24 Comprehensive overviews of major classical and recent studies of the RT instabilities can be found in two extensive reviews by Zhou.11,12 Such instabilities have also been widely used as benchmark problems for evaluating multiphase computational fluid dynamics (CFD) codes and/or numerical schemes, see, for example, Ref. 25. This is mainly because of the occurrence of complex interface dynamics and the formation of small-scale flow patterns.18,26–28 Although the literature about RT instabilities seems saturated, most experimental works have focused on very low Atwood numbers (At0.15) because of some technical difficulties arising due to high-density ratios as discussed in reference.21–23 Waddell et al.21 used a planar laser-induced fluorescence (PLIF) technique to experimentally study the evolution of two-dimensional (2D) single-mode miscible and immiscible RT instabilities with At = 0.155 and At = 0.336, respectively. Different interface accelerations of 0.32,0.74, and 1.34g̃ in excess of the earth's gravitational acceleration g̃ were applied to initiate the formation of the RT instabilities. The initial exponential growth rate was found to be in agreement with earlier linear stability theories, and spike and bubble eventually exhibited constant velocities. Wilkinson and Jacobs22 used the same PLIF setup used by Waddell et al.21 to study three-dimensional (3D) single-mode miscible RT instabilities with At = 0.15, finding the initial growth of the 3D instabilities were in line with the linear stability theories. However, both the spike and bubble of the 3D instabilities reaccelerated after reaching intermediate velocity plateaus which were comparable to those predicted by Goncharov's analytical solution for the spike and bubble terminal velocities.20 More recently, Malamud et al.29 studied experimentally (and also numerically) single-mode RT instabilities with At0.15 and At0.6 during the nonlinear growth regime. The authors observed a marginal reacceleration phenomenon using their laser-based experimental technique (probably for the first time). However, it was concluded that a larger laser facility was required to trigger a full-scale reacceleration for a thorough investigation. Such spike and bubble reaccelerations have also been reported in computational studies, for example, by Ramaprabhu et al.30 who studied 3D single-mode RT instabilities with 0.05At0.9. Later, Ramaprabhu et al.31 showed, based on 3D single-mode miscible RT instabilities with 0.1At0.9, that the reaccelerations of the bubble and spike stopped at very late evolution times due to the intense turbulent mixing. More recently, Liang et al.19 reported spike and bubble reaccelerations for a relatively wide Reynolds number range (10Re5000) based on computational studies of 3D single-mode immiscible RT instabilities with At = 0.15. More recently, Hu et al.13 studied spike and bubble reaccelerations in 2D single-mode miscible RT instabilities with At = 0.15 and 100Re10000. They reported a new deceleration–acceleration stage after the reacceleration stage at medium Reynolds numbers. Zanella et al.32 computationally studied 2D and 3D single-mode immiscible RT instabilities with a very low Atwood number (At0.005) and low-intermediate Reynolds numbers using a diffuse interface model. They also reported a reacceleration behavior for the low-Atwood RT instabilities studied. Lately, Ding et al.33 studied single-mode and dual-mode RT instabilities using a molecular dynamics approach. They found that the microscopic instabilities could exhibit considerable differences in terms of the growth behavior and evolution regimes compared to their macroscopic counterparts. Most previous computational works on RT instabilities focused on either 2D or 3D single-mode instabilities with a limited range and combinations of Atwood and Reynolds numbers.34 Computational studies of 3D multi-mode RT instabilities are still very limited compared to the studies based on the single-mode instabilities.24,34 Yilmaz24 used a large eddy simulation approach to computationally study multi-mode immiscible RT instabilities with 0.5At0.9 and different Reynolds numbers. It was found that interactions of individual instabilities with their surroundings were the primary mechanism for the turbulent mixing in such multi-mode instabilities. Very recently, Livescu et al.35 used direct numerical simulation (DNS) to study the evolution of 3D multi-mode RT instabilities with At=0.5,0.75, and 0.9 after the gravity is set to zero or reversed. They reported significant changes to some turbulence quantities when the acceleration was altered. Overall, the evolution of the 3D multi-mode instabilities has not been directly compared against that of the 3D single-mode instabilities in a quantitative manner.

Considering the existing literature about RT instabilities, it becomes clear that there is still a need for comprehensive comparative studies of the effects of the Atwood and Reynolds numbers on the formation and evolution of the 2D and 3D RT instabilities. Also, a better understanding of the role of the surface tension force in immiscible RT instabilities is needed.3 Additionally, more high-fidelity studies of the 3D multi-mode RT instabilities are required in order to shed light on their complex turbulent mixing characteristics. Therefore, the present study aims to contribute to the field of the incompressible immiscible RT instabilities by performing a quantitative comparative study of the formation and evolution of 2D and 3D single-mode RT instabilities with a wide range of combinations of the Atwood number (0.2At0.75), Reynolds number (300Re10000), surface tension force, and initial perturbation amplitude. The present work additionally presents quantitative comparisons between the formation and evolution of 3D single-mode and multi-mode RT instabilities. Finally, the results reported in this study provide a high-fidelity database beneficial to evaluating future CFD codes and numerical schemes. In order to achieve the aforementioned aims, a two-phase solver based on a phase-field method27,36,37 was developed and used to perform DNS of the RT instabilities. This new solver is implemented within the Xcompact3d framework38 and based on its QuasIncompact3d solver,39 which was developed previously for the variable-density Navier–Stokes equations in the low-Mach number limit. Xcompact3d is a high-order finite-difference CFD framework to study turbulent flows on a Cartesian mesh using supercomputers.40,41

The outline of this paper is as follows: in Sec. II, the governing equations, numerical methods, and problem setup are introduced. Single-mode RT instabilities are then studied in Sec. III where the effects of the Atwood number, Reynolds number, and Weber number (surface tension force) are investigated for 2D and 3D instabilities. In Sec. IV, the evolution of a 3D multi-mode RT instability is studied and compared against those of the single-mode RT instabilities. Conclusions are drawn in Sec. V. In  Appendixes A and  B, the ability of the two-phase solver to model 2D and 3D RT instabilities is evaluated and comparisons are made against relevant literature data. A mesh dependency study is also reported in  Appendix A for completeness.

1. Incompressible Navier–Stokes

The dimensionless governing equations for the flow of two incompressible and immiscible fluids, based on the one-fluid formulation, are given as

ujxj=0,
(2)
uit+12(xj(uiuj)+ujuixj)=1ρ[xj(pδij)+1Rexj(τij)+fiσ]+1Frgi,
(3)

where ui(i=1,2,3) denotes the velocity component in the ith direction (x, y, and z, respectively), ρ represents the density, and p and δij denote the pressure and the Kronecker delta, respectively,42 and the momentum equations are given in skew-symmetric form. The surface tension force and acceleration due to gravity in the ith direction are represented by fiσ and gi, respectively, while τij is the full viscous stress tensor defined as

τij=μ(uixj+ujxi)2μ3ukxkδij,
(4)

noting the second term on the right-hand side is zero due to incompressibility.

The dimensionless scaling terms Re and Fr are the Reynolds number as defined in the introduction, and the Froude number defined as Fr=Ũ2/(g̃L̃). As mentioned earlier, only variables accented with the tilde sign (∼) are dimensional and Ũ,L̃,ρr̃,μr̃, and g̃ denote the reference dimensional velocity, length, density, viscosity, and gravitational acceleration, respectively. The reference velocity scale is determined based on the acceleration due to gravity as Ũ=g̃L̃, which when inserted into the Reynolds number yields the expression Re=ρ̃rL̃3g̃/μ̃r which is used to define the Reynolds number in the present study. The density of the heavy fluid is chosen as the reference density ρr̃=ρh̃, and both fluids are assumed to have a similar viscosity (μr̃=μh̃=μl̃) in the present study.

2. Phase transport equation

There are two main approaches to model two-phase systems based on the Navier–Stokes equations, i.e., the two-fluid and the one-fluid (or mixture) approach. The former is typically beneficial when there is significant mechanical or thermal nonequilibrium between phases and/or compressibility effects.43 For the isothermal incompressible regime considered here the one-fluid approach is appropriate. The one-fluid approach applies a single set of governing equations across the entire domain18 and the boundary conditions at the fluid interface are implicitly imposed using the continuous surface force (CSF) method44 or the ghost fluid method.45 The location of the interface in the one-fluid approach can either be tracked using front-tracking methods46 or interface capturing methods including the level set method,47 the volume-of-fluid (VOF) method,48 the tangent of hyperbola for interface capturing (THINC) method,49 and the diffuse interface (DI) or phase-field methods.50 The phase-field method, in which the interface has a small but finite thickness, provides an interesting alternative that does not suffer from issues related to mass conservation, or the accurate computation of the interface, critical for evaluating the surface tension force, as may be present in level set or VOF methods, respectively.26 It has been reported that the phase-field approach can be potentially more computationally efficient compared to its sharp interface capturing counterparts.51 Various formulations have been used for the phase-field approach; however, they can be categorized into two main categories: the Cahn–Hilliard and Allen–Cahn formulations.51 The former formulation exhibits good robustness and solution stability but may suffer issues stemming from dissipation, over/under-shoots, and coarsening effects.51 Moreover, the transport equation of the Cahn–Hilliard requires calculating fourth-order spatial derivatives which is a quite challenging task and requires complex treatments.18,26,36,52 On the other hand, the Allen–Cahn formulation, which is essentially a convection-diffusion equation with a source term,51 only involves second-order spatial derivatives and is easier to implement numerically.

Following the works of Mirjalili et al.,27 Chiu and Lin,36 Mirjalili and Mani,51 a modified Allen–Cahn formulation was implemented in Xcompact3d as

ϕt+ujϕxj=xj[εγ(εϕxjϕ(1ϕ)1|ϕxi|ϕxj)],
(5)

where ϕ is the phase-field function, ε is a measure of the interface thickness, and γ is a positive coefficient. The phase-field function takes values 0ϕ1 with ϕ=0,1 representing regions of solely heavy or light fluid, respectively; it follows that the local averaged density and viscosity are given by ρ=ϕρl+(1ϕ)ρh and μ=ϕμl+(1ϕ)μh, respectively. Here, following Mirjalili et al.,27 ε=Δ, where Δ is the spatial resolution and γ=max(u12+u22+u32), the maximum velocity magnitude in the domain (evaluated at every time step).

A surface tension force formula based on a CSF model, initially proposed by Brackbill and Kothe44 and followed by, for example, Refs. 36, 37, and 53 is implemented in the present work given as

fiσ=σxj[(1|ϕxk|ϕxj)]ϕxi,
(6)

where the surface tension coefficient is defined as σ=1/We in which We is the Weber number defined as We=ρ1̃Ũ2L̃/σ̃.

The governing equations are spatially discretized on a Cartesian mesh (with uniformly distributed mesh nodes) using high-order finite-difference schemes. In Eqs. (2) and (3), a 6th-order compact finite-difference method with quasi-spectral accuracy41 is employed. Specifically, for an arbitrary variable f the 1st and 2nd spatial derivatives at the jth mesh node in the xi direction are approximated as follows:

αfxi|j1+fxi|j+αfxi|j+1=afj+1fj12Δxi+bfj+2fj24Δxi,
(7)

with α=1/3, a = 14/9, and b = 1/9, and

α2f2xi|j1+2f2xi|j+α2f2xi|j+1=afj+12fj+fj1Δxi2+bfj+22fj+fj24Δxi2,
(8)

with α=2/11, a = 12/11, and b = 3/11. In the present paper, for the uniformly distributed mesh nodes used, fj=f(xj), where xj=(j1)Δxi for 1jnxi where nxi denotes the number of mesh nodes in the xi direction.

Special treatment is needed to discretize the convective term of Eq. (5), the phase-field transport equation, to avoid spurious oscillations and solution instability due to the existence of two-phase interfaces and sudden changes in the fluid properties (density here). For this equation, a 5th-order weighted essentially nonoscillatory (WENO) scheme54 is used, following the strategy employed by McSherry et al.55 in a solver based on the level set method. Various WENO schemes, particularly the one developed by Jiang and Shu,54 have been widely used to solve phase functions (or other phase-related discontinuous functions) in different multiphase modeling approaches.17,26,56 These schemes help to obtain a robust and accurate solution across the interface by suppressing spurious numerical oscillations. Specifically, given a pointwise function f (here, for example, the gradient of the phase function) at a set of equally spaced mesh nodes, the value f̂ at the mesh nodes j can be estimated (in a specific coordinate direction) using an upwind WENO reconstruction as follows:54,57

f̂k=1NkωkIk,
(9)

where the nonlinear weights ωk on the kth stencil candidate (smoothness indicator, βk) defined as18,54,57

ωk=ωk̂k=1Nkωk̂,
(10)

with

ωk̂=γk(βk+εw)2,
(11)

where εw=106 to avoid zero at the denominator. For a five-point WENO reconstruction with k = 3 and a 5th-order accuracy, γ1=0.3,γ2=0.6, and γ3=0.1, and for a positive velocity, the reconstruction is given as18,54,57

I1=13fj+56fj+116fj+2,
(12a)
I2=16fj1+56fj+13fj+1,
(12b)
I3=13fj276fj1+116fj,
(12c)

with the smoothness indicators given as

β1=1312(fj2fj+1+fj+2)2+14(3fj4fj+1+fj+2)2,
(13a)
β2=1312(fj12fj+fj+1)2+14(fj1fj+1)2,
(13b)
β3=1312(fj22fj1+fj)2+14(fj24fj1+3fj)2.
(13c)

It should be noted that the gradients of the phase function ϕ/xi in the right-hand side term of Eq. (5), the surface tension force Eq. (6) and the gradient of the viscosity required to evaluate the full viscous stress tensor are also obtained using the WENO reconstruction. It is worth noting that, as discussed comprehensively by Huang et al.,18 a 5th-order accuracy is only obtainable in smooth regions and the overall accuracy order of a two-phase system can be well below 1st-order due to the discontinuous nature of the phase function.

At each time step, the solver first advances the phase function by solving Eq. (5) and updates the density and viscosity fields. The discretized momentum equation is advanced using a conventional fractional step method. Due to the varying density, a variable-coefficient Poisson equation for the pressure arises when enforcing the continuity constraint Eq. (2). This is reformulated in terms of a correction to a constant-coefficient Poisson equation of the form39,58–60

2pν+1xi2=2pνxi2+ρ̂[1Δtxi(uiui)xi(1ρpνxi)],
(14)

where ui is the intermediate velocity field obtained by solving Eq. (3) without the pressure term. Equation (14) is solved iteratively, with superscript ν the iteration counter, applying a direct spectral Poisson solver to update the pressure field at each iteration. The ρ̂ appearing in Eq. (14) has the effect of an under-relaxation factor with lower values stabilizing the iterative solution process. In this work, ρ̂=min(ρ), as proposed in Ref. 60 due to the large density differences in the problems considered here. After the first two time–steps the pressure field is extrapolated forward in time to form the right hand side of Eq. (14) and a single iteration is performed following the approach of Dodd and Ferrante.42 A separate mesh is used for the pressure field which is staggered by a half mesh spacing with respect to the velocity mesh. This is to avoid the issues related to the pressure-velocity decoupling.41 All governing equations are integrated in time using a two-step explicit Adams–Bashforth method.

Discretizing the divergence of viscous stress tensor Eq. (4), appearing in Eq. (3) (the momentum equations) directly in conservative form can give rise to oscillations in the solution when the viscosity varies.60 Instead, by applying the chain rule, this term can be expressed as

τijxj=μ2uixj2+μxj(uixj+ujxi),
(15)

where incompressibility has again been applied to simplify the expression. The first term is equivalent to the viscous stress tensor in the incompressible case and is likewise discretized using the 6th-order compact scheme for 2nd derivatives. With the exception of the viscosity gradient all other terms are discretized using the 6th-order compact scheme for 1st derivatives. The viscosity is a function of the phase field; therefore the same WENO scheme is used for computing the gradients of the viscosity for consistency with the calculation of the phase-field gradients in the phase-field transport equation.

Similarly to the other Xcompact3d solvers, the new two-phase solver is combined with an efficient highly scalable pencil (or 2D) domain decomposition library for distributed-memory high performance computing platforms, namely, 2Decomp&FFT.61,62 This 2D domain decomposition library is also a highly scalable tool to perform fast Fourier transforms needed for the pressure Poisson equation which is solved in the spectral space for efficiency.41 

1. 2D single-mode

All 2D RT instabilities studied here are modeled in a rectangular domain of size [0,d]×[0,4d] where d is the reference length and the interface is initially located at y0, defined as

y0=2d+Acos(2πxd).
(16)

Equation (16) represents a perturbation with wave number k=2π/d and an amplitude of A=0.1d superimposed on a planar interface. Slip boundary conditions for all flow quantities are applied to all boundaries.27 The gravitational acceleration points downward with a magnitude of g = 1.0. The phase-field function is initialized using

ϕ=12[1.0+tanh(y0y2ε)].
(17)

Various combinations of the Atwood number, Reynolds number, and surface tension are studied, see Table I for the 2D RT instabilities simulations. Moreover, a mesh dependence study is provided in  Appendix A. Finally, the time step for the simulations is given by Δt=2.5×104/At, similar to that used by Ding et al.26 A schematic showing the computational setup, domain dimensions and some key features and parameters of 2D single-mode RT instabilities is provided in Fig. 1 (left).

TABLE I.

Single-mode 2D Rayleigh–Taylor instability cases.

CaseAtρhρlReσ=1WeMesh size
0.5 3000 1×105 65 × 257 
0.5 3000 1×105 129 × 513 
0.5 3000 0.0 129 × 513 
4a 0.5 3000 0.0 129 × 513 
0.5 3000 1×105 257 × 1025 
0.5 1000 1×105 129 × 513 
0.5 300 1×105 129 × 513 
0.5 10 000 1×105 129 × 513 
0.2 1.5 3000 1×105 129 × 513 
10 0.333 3000 1×105 129 × 513 
11 0.666 3000 1×105 129 × 513 
12 0.75 3000 1×105 129 × 513 
13 0.5 3000 1×104 129 × 513 
14 0.5 3000 5×104 129 × 513 
15 0.5 10 000 5×104 129 × 513 
CaseAtρhρlReσ=1WeMesh size
0.5 3000 1×105 65 × 257 
0.5 3000 1×105 129 × 513 
0.5 3000 0.0 129 × 513 
4a 0.5 3000 0.0 129 × 513 
0.5 3000 1×105 257 × 1025 
0.5 1000 1×105 129 × 513 
0.5 300 1×105 129 × 513 
0.5 10 000 1×105 129 × 513 
0.2 1.5 3000 1×105 129 × 513 
10 0.333 3000 1×105 129 × 513 
11 0.666 3000 1×105 129 × 513 
12 0.75 3000 1×105 129 × 513 
13 0.5 3000 1×104 129 × 513 
14 0.5 3000 5×104 129 × 513 
15 0.5 10 000 5×104 129 × 513 
a

RHS of Eq. (5) is neglected for this case.

FIG. 1.

Schematics showing the computational domains, dimensions, and key features of the 2D (left) and 3D (right) single-mode instabilities studied here. All boundaries (highlighted by the red lines) are slip conditions of all quantities. Also, y*=y2d, where y = 0 is at the bottom of the domains. It should be noted that the 3D multi-mode instabilities have a similar domain setup to that shown here for the 3D single-mode instabilities but with a different domain size as provided in Table II.

FIG. 1.

Schematics showing the computational domains, dimensions, and key features of the 2D (left) and 3D (right) single-mode instabilities studied here. All boundaries (highlighted by the red lines) are slip conditions of all quantities. Also, y*=y2d, where y = 0 is at the bottom of the domains. It should be noted that the 3D multi-mode instabilities have a similar domain setup to that shown here for the 3D single-mode instabilities but with a different domain size as provided in Table II.

Close modal
TABLE II.

Single-mode and multi-mode 3D Rayleigh–Taylor instability cases. For the multi-mode RT instabilities, the amplitude corresponds to the velocity field perturbation.

CaseModeAtρhρlReσ=1WeADomain size (×d)Mesh size
Single 0.5 1024 1×105 0.05 1×4×1 129×513×129 
Single 0.5 3000 1×105 0.05 1×4×1 129×513×129 
Single 0.5 3000 1×105 0.1 1×4×1 129×513×129 
Single 0.2 1.5 1024 1×105 0.05 1×4×1 129×513×129 
Single 0.2 1.5 1000 1×105 0.1 1×4×1 129×513×129 
Single 0.333 1000 1×105 0.1 1×4×1 129×513×129 
Single 0.5 1000 1×105 0.1 1×4×1 129×513×129 
Single 0.666 1000 1×105 0.1 1×4×1 129×513×129 
Single 0.75 1000 1×105 0.1 1×4×1 129×513×129 
10 Single 0.2 1.5 3000 1×105 0.1 1×4×1 129×513×129 
11 Single 0.333 3000 1×105 0.1 1×4×1 129×513×129 
12 Single 0.666 3000 1×105 0.1 1×4×1 129×513×129 
13 Single 0.75 3000 1×105 0.1 1×4×1 129×513×129 
14 Single 0.5 3000 1×105 0.005 1×4×1 129×513×129 
15 Multi 0.5 3000 1×105 0.005 π2×3π×π2 301×1801×301 
16 Multi 0.5 3000 1×105 0.005 π2×3π×π2 401×2401×401 
17 Multi 0.5 3000 1×105 0.005 π2×3π×π2 501×3001×501 
CaseModeAtρhρlReσ=1WeADomain size (×d)Mesh size
Single 0.5 1024 1×105 0.05 1×4×1 129×513×129 
Single 0.5 3000 1×105 0.05 1×4×1 129×513×129 
Single 0.5 3000 1×105 0.1 1×4×1 129×513×129 
Single 0.2 1.5 1024 1×105 0.05 1×4×1 129×513×129 
Single 0.2 1.5 1000 1×105 0.1 1×4×1 129×513×129 
Single 0.333 1000 1×105 0.1 1×4×1 129×513×129 
Single 0.5 1000 1×105 0.1 1×4×1 129×513×129 
Single 0.666 1000 1×105 0.1 1×4×1 129×513×129 
Single 0.75 1000 1×105 0.1 1×4×1 129×513×129 
10 Single 0.2 1.5 3000 1×105 0.1 1×4×1 129×513×129 
11 Single 0.333 3000 1×105 0.1 1×4×1 129×513×129 
12 Single 0.666 3000 1×105 0.1 1×4×1 129×513×129 
13 Single 0.75 3000 1×105 0.1 1×4×1 129×513×129 
14 Single 0.5 3000 1×105 0.005 1×4×1 129×513×129 
15 Multi 0.5 3000 1×105 0.005 π2×3π×π2 301×1801×301 
16 Multi 0.5 3000 1×105 0.005 π2×3π×π2 401×2401×401 
17 Multi 0.5 3000 1×105 0.005 π2×3π×π2 501×3001×501 

2. 3D single-mode and multi-mode

The single-mode 3D RT instabilities studied here are modeled in a [0,d]×[0,4d]×[0,d] domain where the interface is initially located at y0, defined as

y0=2d+A[cos(2πxd)+cos(2πzd)].
(18)

With A = 0.05, Eq. (18) provides an initial maximum perturbation amplitude comparable to that produced by Eq. (16) with A = 0.1 in 2D. The phase-field function is initialized using Eq. (17). A slip boundary condition is applied to all boundaries of the domain with a 129×513×129 mesh size. The gravitational acceleration points downward with a magnitude of g = 1.0. The parameters of the simulations are presented in Table II. Similarly to the 2D cases a time step of Δt=2.5×104/At is used for the 3D cases. A schematic showing the computational setup, domain dimensions, and key features of 3D single-mode RT instabilities is provided in Fig. 1 (right).

In addition to the single-mode RT instabilities, one 3D multi-mode RT instability with At = 0.5 and Re = 3000 is also studied in a domain size of [0,πd/2]×[0,3πd]×[0,πd/2] discretized with three mesh resolutions, with 301×1801×301,401×2401×401, and 501×3001×501 mesh nodes (cases 15–17 of Table II, respectively). The interface is initially located at y0=3πd/2 and, unlike the 2D cases, an initial perturbation is applied to the velocity field. Specifically, the vertical component of the velocity vector is initialized as24 

u2=Aβ(1+cos(2πxd)),
(19)

where β is a random number between −1 and 1, and the amplitude A = 0.005. The random numbers are created based on the random_number() routine of FORTRAN and applied to all mesh nodes. By using Eq. (19), the amplitude of the perturbations is smoothed toward the vertical boundaries. Directly perturbing the velocity field to develop 3D multi-mode RT instabilities is superior to perturbing the interface amplitude as the former approach produces resolution-independent perturbations without additional “mesh smoothing” requirements (for when the amplitude of the interface perturbation is smaller than the mesh resolution).24,63 The perturbation spectrum is a white noise down to the mesh scale, meaning that, at very early times when high wavenumber modes are dominant, there may be slight differences between different simulations. However, soon after when the nonlinear mode interactions start, the footprint of the initial perturbations is no longer impacting the flow.63 

1. General features

In a system of two fluids with a perturbed interface, as shown in Fig. 2 (and Fig. 16 in  Appendix A), the heavier fluid on the top penetrates the lighter fluid due to the force of gravity, a motion which consequently moves the lighter fluid upward. This behavior results in the formation of the spike and bubble features which are the tip of the falling and rising fluids, respectively. As the solution evolves, the penetration length of the spike increases, and it eventually rolls up and forms two transient counter-rotating vortices due to shear forces and Kelvin–Helmholtz type instabilities.15 At around t*=tAt=2.25, a pair of secondary vortices appear at the tails of the spike vortices which then evolve into relatively complex chaotic flow features and promote turbulent mixing at later times.15,17,26,64 The interface remains symmetric despite its complex evolution.

FIG. 2.

Top: effect of the Reynolds number on the evolution of the 2D single-mode RT instability with At = 0.5. Bottom: temporal evolution of the vertical coordinates of the tip of the falling (spike) and rising (bubble) fluid of the 2D single-mode RT instability with At = 0.5 (the order of the lines is similar to Fig. 16 in  Appendix A).

FIG. 2.

Top: effect of the Reynolds number on the evolution of the 2D single-mode RT instability with At = 0.5. Bottom: temporal evolution of the vertical coordinates of the tip of the falling (spike) and rising (bubble) fluid of the 2D single-mode RT instability with At = 0.5 (the order of the lines is similar to Fig. 16 in  Appendix A).

Close modal

2. Reynolds number effect

In this subsection, the effect of the Reynolds number on the evolution of the 2D single-mode RT instabilities with At = 0.5 is investigated. In addition to the case with Re = 3000, three more test cases with Re=300,1000, and 10 000 are studied. Figure 2 (top) shows the temporal evolution of the 2D RT cases with different Reynolds numbers up to tAt=3.75. Figure 2 (bottom) shows the coordinates of the spike and bubble vs tAt. It should be noted that y*=y2 in this study. From Fig. 2, it is evident that increasing the Reynolds number does not noticeably affect the penetration lengths of the spike and bubble features except for the penetration length of the spike of the case with Re = 300 which is slightly lower compared to those of the other cases. A comparable trend for the sensitivity to the Reynolds number was also reported in Liang et al.19 for 3D single-mode RT instabilities with At = 0.15. The Reynolds number is found to have a substantial influence on the topology of the transient spike vortices and the chaotic region created due to the spike rolling up. The counter-rotating vortices appear earlier as the Reynolds number increases which can be seen for example in the snapshots tAt=1.75 of Fig. 2 (left). The secondary vortices at the tail of the initial vortices also appear earlier as the Reynolds number increases which can be seen in snapshots tAt=2.25 of Fig. 2. The cases with Re = 300 and 1000 do not develop the secondary vortices until very late, only becoming apparent at tAt=3.75. With Re = 10 000, the bubble is perturbed and vortical structures are formed within the bubble region and in the vicinity of the stem of the spike, which grows over time creating complex features; despite this, the RT instability remains symmetric at increased Reynolds numbers.

3. Atwood number effect

In this subsection, the effect of the Atwood number on the evolution of the 2D single-mode RT instabilities with Re = 3000 is investigated. In addition to the initial case with At = 0.5, four more test cases with At=0.2,0.333,0.666, and 0.75 are also considered, corresponding to a density ratio of ρh/ρl=1.5,2,5, and 7, respectively. Figure 3 shows the temporal evolution of the RT cases with different Atwood numbers up to tAt=3.0 while Fig. 4 shows the temporal evolution of the coordinates of the spike and bubble. By increasing the Atwood number, the penetration lengths of the spike and bubble features increase and decrease, respectively, as seen in Figs. 3 and 4. From Fig. 4, which also shows the absolute penetration lengths of the spike |ys*| and bubble |yb*| vs the Atwood number at tAt=3.0, it can be seen that the spike exhibits a significant nonlinear growth as the Atwood number increases, while the bubble behaves almost linearly. Increasing the Atwood number reduces the significance of the transient counter-rotating vortices at the spike tip and consequently reduces the intensity of the chaotic flow features at the tail of the secondary vortices. This trend is consistent with some previous observations, for example, in Ref. 17. Specifically, the latter authors reported that for a high Atwood number (ρh/ρl=30) the spike roll-up was very limited and for a very high Atwood number (ρh/ρl=1000) the spike did not roll up. This is mainly attributed to significant increases in the magnitude of the spike tip velocity as the Atwood number increases.

FIG. 3.

Effect of the Atwood number for the 2D single-mode RT instability at Re = 3000.

FIG. 3.

Effect of the Atwood number for the 2D single-mode RT instability at Re = 3000.

Close modal
FIG. 4.

Top left: evolution of the vertical coordinates of the spike and bubble for the 2D RT instability at Re = 3000. Top right: evolution of the velocities of the spike and bubble (with a zoomed-in view mainly for the bubble velocities shown at the very top of the figure). Bottom left: |ys*|,|yb*|, and |ysb*| as a function of At. Bottom right: |ysb*| as a function of At.

FIG. 4.

Top left: evolution of the vertical coordinates of the spike and bubble for the 2D RT instability at Re = 3000. Top right: evolution of the velocities of the spike and bubble (with a zoomed-in view mainly for the bubble velocities shown at the very top of the figure). Bottom left: |ys*|,|yb*|, and |ysb*| as a function of At. Bottom right: |ysb*| as a function of At.

Close modal

Figure 4 also shows the velocity U=d|y*|/dt of the spike and bubble in the 2D RT instabilities with different Atwood numbers. During the initial exponential growth stage, the velocities of the spike and bubble increase for all Atwood numbers as the solution evolves. Based on the linear Rayleigh–Taylor instability theory of two irrotational incompressible fluids, the initial growth of the RT amplitude can be estimated as19,21,22

|ysb*|=|ys*|+|yb*|2=a1eUt+a2eUt,
(20)

where |ys*| and |yb*| are the absolute vertical distance of the spike and bubble from y*=0 and a1 and a2 are the fitting coefficients. For very low surface tension values, the growth rate (velocity) of the initial stage is estimated as U=Atgk. However, after the initial growth, the spike and bubble velocities almost level off at semi-terminal velocities before they reaccelerate. The bubble terminal velocity can be estimated using a correlation derived by Goncharov20 as

Ub=2At1+AtgCk,
(21)

where C = 3 for 2D RT instabilities and C = 1 for 3D instabilities (k=2π is the perturbation wavenumber). The spike terminal velocity can be estimated (less accurately) using a similar correlation but replacing 1+At with 1At in the denominator of Eq. (21).20 The horizontal dashed and dotted lines in the U graph of Fig. 4 show the spike and bubble terminal velocities, estimated by the appropriate form of Eq. (21), respectively. A good agreement is seen between the terminal velocities obtained in the present simulations and Eq. (21), particularly for the bubble. Similar trends for the bubble and spike velocities have been reported previously.19,20

As mentioned previously, the bubble (and spike) reacceleration in RT instabilities has been reported in the literature19,22,30,65 where it is suggested that the formation of the spike tip vortices plays a key role in such reaccelerations. Specifically, such vortices act to reduce the friction drag on the spike. They could also produce secondary vortices and consequently a momentum jet that can propel the bubble forward.30 From Fig. 4, it is clear that the reacceleration of both the spike and bubble starts at around tAt=2. On the other hand, from Fig. 3, the spike tip vortices are fully developed in all RT instabilities at tAt=2 and for later times, the secondary chaotic structures (or the momentum jets) move upward toward the bubble region. Ramaprabhu et al.30 reported such reaccelerations for RT instabilities with At0.75 and showed that RT instabilities with At = 0.8 and At = 1.0 did not exhibit such flow behavior. RT instabilities with At0.8 do not form strong spike tip vortices, something that again points to the key role of such vortical structures in the reacceleration mechanism. Therefore, based on the current observations, it can be concluded that the vortical structures of the spike tip do play a key role in the reacceleration of the spike and bubble in the 2D RT instabilities with At0.75.

Figure 4 also shows the vertical distance of the bubble and spike |ysb*| vs the Atwood number following the initial growth of the 2D single-mode RT instabilities (tAt1.5). By considering the exponential relationship between |ysb*| and At as seen from Eq. (20), a relationship between the vertical distance of the bubble and spike, and the Atwood number for tAt1.5 is proposed as follows:

|ysb*|=η[θ1e(ψkAtP)+θ2e(ψkAtP)],
(22)

where θ1=0.653,θ2=0.362, k is the initial perturbation wavenumber (2π here), P = 2, and ψ is the time coefficient defined as ψ=0.25/At|0.5. The adjustment factor η=0.47,0.635,0.805,1 for tAt=1.5,2,2.5, and 3, respectively. The θ parameters are selected based on the best fit at tAt=3, and the adjustment parameter is introduced to shift the curves and avoid introducing new θ parameters for other times. The mathematical reason for having a time coefficient can be hinted from Eq. (20) where the time is inside the exponential function. Here, a global time coefficient ψ is defined based on the time step used for the RT instability with At = 0.5 with ψ=1000Δt. Figure 4 shows that Eq. (22) predicts the vertical distance of the bubble and spike quite accurately compared to the present DNS data particularity at later times [the blue dashed lines superimposed on the |ysb*| curves represent the predictions obtained by Eq. (22) with the parameters given in Table III]. It should be noted that η = 1 for tAt=3 since its curve is chosen as the reference. Note that other combinations of θ1, θ2, ψ, and η may exist and the focus here is the general exponential form of Eq. (22).

TABLE III.

Parameters of Eq. (22) used for the 2D and 3D single-mode RT instabilities.

DimensiontAtReηθ1θ2ψP
2D 1.5 3000 0.47 0.653 0.362 0.353 
2D 3000 0.635 0.653 0.362 0.353 
2D 2.5 3000 0.805 0.653 0.362 0.353 
2D 3000 0.653 0.362 0.353 
3D 1000 and 3000 1.12 0.653 0.362 0.353 
3D 2.25 1000 and 3000 1.28 0.653 0.362 0.353 
DimensiontAtReηθ1θ2ψP
2D 1.5 3000 0.47 0.653 0.362 0.353 
2D 3000 0.635 0.653 0.362 0.353 
2D 2.5 3000 0.805 0.653 0.362 0.353 
2D 3000 0.653 0.362 0.353 
3D 1000 and 3000 1.12 0.653 0.362 0.353 
3D 2.25 1000 and 3000 1.28 0.653 0.362 0.353 

4. Weber number (surface tension) effect

In this subsection, the effect of the surface tension on the evolution of the 2D single-mode RT instabilities with various Atwood and Reynolds numbers are investigated. The effect of surface tension is studied by keeping the Atwood and Reynolds numbers constant with At = 0.5 and Re = 3000. Surface tension values of σ=1/We=0,1×105,1×104, and 5×104 are investigated. For σ = 0 an additional case is studied in which the right hand side term of Eq. (5) is neglected (case 4 of Table I). As mentioned earlier, the phase-field method can be based on the Allen–Cahn or Cahn–Hilliard formulations which mainly differ in the diffusion term of their phase transport equation.16,51 Unlike the Cahn–Hilliard formulation, the Allen–Cahn transport formulation is independent of chemical potential parameters including the surface tension.16 The chemical potential Φ of the Cahn–Hilliard approach can be defined as Φ=(1/ε)σαζεσαΔϕ,26 where ζ is the bulk energy density and α is a constant. If the surface tension is zero, the chemical potential term and consequently the right-hand side of a Cahn–Hilliard equation becomes zero yielding a transient advection equation. The reason for introducing test case 4 in this study is to examine the effect of the right-hand side term of the Allen–Cahn method in the absence of the surface tension force.

Figure 5 (left) shows the temporal evolution of the 2D RT instabilities with different values of the surface tension coefficient as mentioned above for times up to tAt=3.5. No noticeable differences are apparent in the evolution of the spike and bubble features and their absolute penetration lengths. The only visible difference is in the structure of the secondary vortices and the associated chaotic tail for the case with σ=5×104 after tAt=2.0. It confirms that, unless the surface tension is relatively significant (low Weber numbers), the RT instability with At = 0.5 and Re = 3000 can be safely modeled with σ = 0 which has been commonly done in most previous studies of this specific RT instability. Another important conclusion is that, although the diffusion term of the Allen–Cahn formulation is independent of surface tension, the right-hand side can be safely neglected in the absence of the surface tension or when its value is very low (high We values), which can noticeably reduce the computational time per time step. Neglecting the right-hand side of Eq. (5) cuts the cost per time step by 7.5% (comparing cases 3 and 4 of Table I). In addition to the above studies, the impact of surface tension on the case with At = 0.75 has also been studied although graphical results are not shown for brevity. Similarly to the case with At = 0.5, the surface tension mainly affects the secondary vortices and their associated chaotic tail, and it does not have a noticeable effect on the evolution and absolute penetration lengths of the spike and bubble features.

FIG. 5.

Left: effect of the surface tension (cases 3, 4, 2, 13, and 14 of Table I) on the evolution of the 2D RT instability for At = 0.5 and Re = 3000. Right: effect of the surface tension on the evolution of the 2D RT instability for At = 0.5 and Re = 10 000.

FIG. 5.

Left: effect of the surface tension (cases 3, 4, 2, 13, and 14 of Table I) on the evolution of the 2D RT instability for At = 0.5 and Re = 3000. Right: effect of the surface tension on the evolution of the 2D RT instability for At = 0.5 and Re = 10 000.

Close modal

Extending the study to a higher Reynolds number regime, the effect of surface tension in cases with Re = 10 000 is investigated. Figure 5 (right) shows a direct comparison between the temporal evolution of the RT instabilities with Re = 10 000 and σ=1×105 and 5×104. A relatively high surface tension is found to have a significant effect on the evolution of the RT instability. In the case with σ=5×104, it is observed that the shape of the bubble region and the stem of the spike are much closer to those with the same Atwood number but lower Reynolds numbers. This is because a high surface tension force could resist the shear forces and consequently cancels the formation of the Kelvin–Helmholtz type instabilities that arise in the bubble region and also the stem section of the spike.

1. General features

To illustrate the evolution of a 3D single-mode RT instability, Fig. 6 shows snapshots of the interface (defined as the iso-surface ϕ=0.5) for the case At = 0.5 and Re = 1024 at several time instants. The temporal evolution of the same 3D RT instability on the z = 0.5 and z = x planes are shown in Fig. 7. Flow features shown in the latter two figures compare very well with those reported by He et al.66 and Lee and Kim.67 As seen in Fig. 6, the spike tip vortices discussed previously in the context of 2D RT instabilities form as a vortex ring in the 3D case. The saddle point, as highlighted in Fig. 6, is a specific feature of 3D RT instabilities. Another key difference between the 2D and 3D RT instabilities is the formation of a secondary (two-layer) roll-up that forms a “girdle”66 around the spike inside the bubble which is seen clearly in the x = z planes of Fig. 7 at tAt=2.5 and 3. Similarly to the vortex ring, the secondary roll-up is also formed due to the existence of Kelvin–Helmholtz type instabilities.66 The interface of the 3D single-mode RT instability, like the 2D single-mode RT instabilities, maintains its symmetries despite the evolution of complex features. Quantitative comparisons between the present DNS study and those reported by He et al.,66 Lee and Kim67 based on the coordinates of the spike, bubble, and saddle point of the RT instability shown in Fig. 6 and another 3D single-mode RT instability with At = 0.2 and Re = 1024 are provided in  Appendix B for completeness.

FIG. 6.

Temporal evolution of the interface of the 3D single-mode RT instability for At = 0.5 and Re = 1024. From left to right, the snapshots correspond to tAt=1,2, and 2.5.

FIG. 6.

Temporal evolution of the interface of the 3D single-mode RT instability for At = 0.5 and Re = 1024. From left to right, the snapshots correspond to tAt=1,2, and 2.5.

Close modal
FIG. 7.

Temporal evolution of the 3D single-mode RT instability for At = 0.5 and Re = 1024 on the z = 0.5 (left) and z = x (right) planes.

FIG. 7.

Temporal evolution of the 3D single-mode RT instability for At = 0.5 and Re = 1024 on the z = 0.5 (left) and z = x (right) planes.

Close modal

2. Reynolds number and amplitude effects

In this subsection, the effects of the Reynolds number and the initial perturbation amplitude on 3D single-mode RT instabilities are investigated for instabilities with At = 0.5, Re=1000,1024, and 3000, and A=0.005,0.5, and 0.1 (cases 1–4 and 14 of Table II). Figure 8 shows the coordinates of the spike, bubble, and saddle point of these 3D instabilities as a function of tAt. Similarly to the 2D single-mode RT instabilities, with a specific initial perturbation amplitude, varying the Reynolds number does not significantly affect the absolute penetration lengths of the spike and bubble (and also the saddle point) of the 3D single-mode RT instabilities. On the other hand, changing the initial amplitude has a noticeable effect as shown in Fig. 8. For instance, with Re = 3000 and at tAt=2.5, the RT instability with A = 0.1 exhibits 15% and 25% higher absolute penetration lengths of the spike and bubble, respectively, compared with the RT instability with A = 0.05. The RT instability with A = 0.005, unlike the other RT instabilities with much higher initial amplitudes, does not show an exponential initial growth, as also predicted by linear stability theories.1 However, after the initial growth and after around tAt=1.5, all RT instabilities show comparable growth rates for the spike and bubble as shown in Fig. 8.

FIG. 8.

Left: temporal evolution of the vertical coordinates of the spike, bubble, and saddle point for the 3D single-mode RT instabilities with At = 0.5 (the order of the lines is similar to Fig. 18 in  Appendix B). Right: temporal evolution of the velocities of the spike and bubble for the 3D single-mode RT instabilities with At = 0.5 (the horizontal dashed and dotted red lines show the predictions of Eq. (21) for the spike and bubble terminal velocities, respectively).

FIG. 8.

Left: temporal evolution of the vertical coordinates of the spike, bubble, and saddle point for the 3D single-mode RT instabilities with At = 0.5 (the order of the lines is similar to Fig. 18 in  Appendix B). Right: temporal evolution of the velocities of the spike and bubble for the 3D single-mode RT instabilities with At = 0.5 (the horizontal dashed and dotted red lines show the predictions of Eq. (21) for the spike and bubble terminal velocities, respectively).

Close modal

Figure 9 shows a direct comparison between the penetration lengths of the spike and bubble of the 2D and 3D RT instabilities with At = 0.5 and Re = 3000. As mentioned earlier, the 3D case with A = 0.05 has the same maximum initial perturbation amplitude as the 2D case with A = 0.1. From Fig. 9 it is clear that with similar Atwood and Reynolds numbers, a 3D RT instability exhibits larger absolute penetration lengths of the spike and bubble compared to a 2D RT instability with the same initial perturbation magnitude. For example, at tAt=2.5, the 3D case with Re = 3000 has 43% and 51% higher absolute penetration lengths of the spike and bubble, respectively, compared with the 2D case. Such differences in the absolute penetration lengths between the 2D and 3D cases have been reported previously, for example, in Bian et al.65 A similar prediction can also be made following Goncharov's equation for the terminal velocity [Eq. (21)] which suggests that for a 2D case the bubble terminal velocity is around 43% lower than that of a comparable 3D case.

FIG. 9.

Temporal evolution of the vertical coordinates of the spike and bubble. Direct comparison between 2D and 3D RT instabilities with At = 0.5, Re = 3000, and the same maximum initial perturbation amplitude.

FIG. 9.

Temporal evolution of the vertical coordinates of the spike and bubble. Direct comparison between 2D and 3D RT instabilities with At = 0.5, Re = 3000, and the same maximum initial perturbation amplitude.

Close modal

The velocity of the spike and bubble in the 3D single-mode RT instabilities with different Reynolds numbers and perturbation amplitudes studied in this subsection are shown in Fig. 8 (right). Unlike the 2D RT instabilities, no semi-terminal velocity is observed for the spikes of the 3D instabilities; however, the bubbles do reach a semi-terminal velocity which are in good agreement with the value obtained using Eq. (21) indicated by the horizontal lines in Fig. 8. With the exception of the case with A = 0.005, which has a significantly different initial growth rate, the spike velocities of the cases with different initial perturbation amplitudes approach the same value at the end of the simulation. The terminal bubble velocity of the case with A = 0.005 is 10%15% lower than those of the cases with A = 0.05 and 0.1. Similarly to the 2D instabilities, the bubble of the 3D instabilities with A = 0.1 reaccelerate while those of the 3D instabilities with A = 0.005 and 0.05 do not exhibit such behavior. This suggests an increase in the nonlinearity of the RT instabilities as the perturbation amplitude increases as also predicted by linear stability studies.1 As mentioned previously, the 3D RT instability with A = 0.05 has a maximum initial perturbation amplitude equal to that of the 2D instability with A = 0.1. Therefore, based on the data in Fig. 8, it can be concluded that a 2D RT instability exhibits bubble (and spike) reacceleration (i.e., higher nonlinearity) with lower initial perturbation amplitudes compared to a 3D RT instability with similar Atwood and Reynolds numbers.

3. Atwood number effect

In this subsection the effect of the Atwood number on the evolution of the 3D single-mode RT instabilities at Re = 1000 and Re = 3000 is investigated with an initial perturbation amplitude A = 0.1 to include a strong nonlinear growth feature. For these cases, the Atwood numbers At=0.2,0.333,0.5,0.666, and 0.75 are studied. Figure 10 shows the temporal evolution of the coordinates of the spike, bubble, and saddle point. At both Re = 1000 and Re = 3000, the instabilities shown in Fig. 10 have very similar absolute spike, bubble, and saddle penetration lengths, in line with the previous observations for the 2D RT instabilities. As for the 2D cases, by increasing the Atwood number, the absolute penetration lengths of the spike and bubble features increase and decrease, respectively. An increased penetration length of the 3D saddle feature with increasing Atwood number is also observed. Figure 10 also shows the velocities of the spike and bubble in the 3D RT instabilities with Re = 1000 and different Atwood numbers. It can be observed that the velocities of the spike (for cases with At0.333) and bubble increase as the solution evolves. Then the spike and bubble velocities almost level off (at semi-terminal velocities) before they reaccelerate again at around tAt=2 like the 2D RT instabilities. However, the spike of the 3D cases with At=0.5,0.666, and 0.75 do not reach an initial plateau and show an increasing trend (except for the case with At = 0.666 which decelerates after tAt=2 since the spike reaches the lower boundary very soon afterward).

FIG. 10.

Top: temporal evolution of the vertical coordinates of the spike, bubble, and saddle point for the 3D single-mode RT instabilities with Re = 1000 and 3000 (the line groupings from top to bottom are: bubble, saddle, and spike coordinates). Bottom left: temporal evolution of the velocities of the spike and bubble fluid for the 3D single-mode RT instabilities with Re = 1000. Bottom right: vertical distance of the spike and bubble as a function of At for Re = 1000 and 3000 at tAt=1.5,1.75,2, and 2.25. The blue dashed lines superimposed on the |ysb*| curves represent the predictions obtained by Eq. (22) with parameters given in Table III.

FIG. 10.

Top: temporal evolution of the vertical coordinates of the spike, bubble, and saddle point for the 3D single-mode RT instabilities with Re = 1000 and 3000 (the line groupings from top to bottom are: bubble, saddle, and spike coordinates). Bottom left: temporal evolution of the velocities of the spike and bubble fluid for the 3D single-mode RT instabilities with Re = 1000. Bottom right: vertical distance of the spike and bubble as a function of At for Re = 1000 and 3000 at tAt=1.5,1.75,2, and 2.25. The blue dashed lines superimposed on the |ysb*| curves represent the predictions obtained by Eq. (22) with parameters given in Table III.

Close modal

Figure 10 also shows the vertical distance of the bubble and spike ysb* as a function of Atwood number for the 3D single-mode RT instabilities. It is clear that Eq. (22), with similar values of θ1, θ2, and ψ used for the 2D RT instabilities, can also be used for the 3D instabilities with Re = 1000 and Re = 3000 studied here particularly for tAt2 as shown in Fig. 10. However, for the 3D instabilities P = 4 instead of P = 2 in the 2D instabilities. This is related to higher spike and bubble velocities of a 3D instability compared to its 2D counterpart as mentioned earlier. Specifically, within its very initial growth stage, a 3D instability penetrates a significant portion of the domain length, and as a result, it exhibits less exponential behavior compared to a 2D case [please note that the Atwood number is less than unity and a higher P value results in a lower exponential value in Eq. (22)].

Cross sections of the interface of the 3D single-mode RT instabilities with Re = 1000 and 3000 and different Atwood numbers at tAt=2 are shown on the z = x planes in Fig. 11 and also in 3D in Fig. 12. While the penetration lengths of the spike and bubble were observed to be very similar, the 3D single-mode RT instabilities with different Atwood number values exhibit quite different flow features, particularly with respect to the spike vortex ring and the secondary roll-ups. As with the 2D instabilities, in the 3D instabilities, the vortical structures become weaker as the Atwood number increases. The interface area Aint permits a more quantitative means of comparison, its evolution as a function of both tAt and t is shown in Fig. 13. For Re = 1000, all lines almost collapse on top of each other for the Aint vs tAt graph. However, the Aint vs t graph of Fig. 13 provides clearer trends between the interface areas of the instabilities with different Reynolds and Atwood numbers. It is observed that, at a specific time, a higher Reynolds number results in the formation of a slightly larger interface area. Also, at a specific time, a lower Atwood number results in a larger interface area. This behavior is attributed to the formation of larger and wider vortical structures in the cases with higher Reynolds number and/or lower Atwood number as can be seen in Figs. 11 and 12.

FIG. 11.

Location of the interface (ϕ=0.5) for the 3D single-mode RT instabilities with different Atwood numbers for Re = 1000 (red) and Re = 3000 (black) on the z = x plane at tAt=2.

FIG. 11.

Location of the interface (ϕ=0.5) for the 3D single-mode RT instabilities with different Atwood numbers for Re = 1000 (red) and Re = 3000 (black) on the z = x plane at tAt=2.

Close modal
FIG. 12.

Interface (ϕ=0.5) of the 3D single-mode RT instabilities for Re = 1000 (top) and Re = 3000 (bottom) at tAt=2.

FIG. 12.

Interface (ϕ=0.5) of the 3D single-mode RT instabilities for Re = 1000 (top) and Re = 3000 (bottom) at tAt=2.

Close modal
FIG. 13.

Evolution of the interface area of the 3D single-mode RT instabilities for Re = 1000 and Re = 3000 for different Atwood numbers as a function of tAt (left) and t (right).

FIG. 13.

Evolution of the interface area of the 3D single-mode RT instabilities for Re = 1000 and Re = 3000 for different Atwood numbers as a function of tAt (left) and t (right).

Close modal

In this section, a 3D multi-mode RT instability is investigated thanks to highly resolved DNS with up to 750 ×106 mesh nodes. The mesh resolutions have been chosen based on previous numerical works24,35,68 with the aim to resolve as many turbulent scales as possible for the chosen combination of Reynolds and Atwood numbers. The aim here is to compare 3D single-mode and multi-mode RT instabilities but also to discuss the influence of the mesh resolution when individual instabilities of a multi-mode RT instability merge and collapse as they develop. Three DNS are performed with three different resolutions as seen in Table II. As discussed in Sec. II C 2, unlike the single-mode RT instabilities, the initial perturbation is applied to the velocity field. This perturbation results in the formation of many, almost uniformly distributed, single-mode type RT instabilities as can be seen in Fig. 14 where the temporal evolution of the interface (defined as ϕ=0.5 iso-surface) of the 3D multi-mode RT instability with A = 0.005, At = 0.5, and Re = 3000 is plotted for three simulations with 301×1801×301,401×2401×401, and 501×3001×501 mesh nodes, respectively. As the solution evolves, the initial individual RT instabilities grow rapidly, forming elongated instabilities with very irregular shapes. At later times, the individual irregular RT instabilities start merging and collapsing forming a chaotic (turbulent) field as seen in the last few snapshots of Fig. 14. Qualitatively comparable flow patterns to those shown in Fig. 14 have been previously reported for multi-mode RT instabilities.24,35,68,69 A visual inspection of the interface seems to suggest that it is not affected by the mesh resolution up to t*=3.

FIG. 14.

Time evolution of the interface (iso-surface of ϕ=0.5) of the 3D multi-mode RT instability with Re = 3000 and At = 0.5 for different mesh resolutions.

FIG. 14.

Time evolution of the interface (iso-surface of ϕ=0.5) of the 3D multi-mode RT instability with Re = 3000 and At = 0.5 for different mesh resolutions.

Close modal

To provide a more quantitative comparison, the evolution of the interface area Aint normalized by the domain cross section area Ad as a function of t*=tAt for the 3D multi-mode RT instability is compared with those of the 3D single-mode cases with A = 0.005 and 0.05 in Fig. 15. The interface growth of the multi-mode RT instability shows a similar exponential trend to the single-mode cases until t*=3 when it slows down and departs the single-mode behavior. This behavior is a consequence of the individual building block RT instabilities of the multi-mode RT instability. As these building block instabilities grow, the surface of the multi-mode RT instability increases, however, as shown in Fig. 14, the individual RT instabilities start merging and collapsing, resulting in a reduction of the growth rate of the total interface area. Specifically, the irregular chaotic instabilities, formed by the initial individual single-mode RT instabilities, continue to grow, but relatively bulkier with less wrinkled surfaces irregular instabilities, as can be seen in t*>3 snapshots of Fig. 14. At later times when individual instabilities start merging and collapsing, the interface area of the multi-mode RT instability increases as the mesh becomes finer as shown in Fig. 15. This is attributed to the formation of more small-scale irregular instabilities when the mesh resolution is increased. Such small instabilities merge at much later times compared to their larger counterparts and therefore help the instability maintain a more wrinkled surface, hence a larger surface area. The interface area graph of Fig. 15 suggests that a much finer mesh than the finest one used here would lead to an even bigger interface area. It could suggest that the mesh resolutions in the present study might not be enough to capture properly the interface area after the initial exponential growth. In order to investigate the influence of the mesh resolution on the 3D multi-mode RT instability, the ratio of the kinetic energy KE to released potential energy PE has been used previously, see, for instance, Ref. 24. Figure 15 provides a comparison of the temporal evolution of KE/PE with KE=1/2ΩuiuidΩ and PE=Ω[ρ(xi,0)ρ(xi,t)]x2dΩ, where Ω denotes the computational domain volume.24 It can be seen that for the present study, KE/PE has a similar trend for the three mesh resolutions used in the present study, suggesting that even the coarsest mesh used here is able to capture the kinetic and released potential energies to reasonable levels. It is worth noting that after t*=3, KE/PE reaches values around 0.4 (for all mesh resolutions) which is comparable to the values previously reported for similar multi-mode instabilities, for example, to the one studied by Yilmaz.24 The behavior observed in Figs. 14 and 15 clearly shows the significant influence of the growth mechanism of the single-mode RT instabilities and of the mesh resolution on the overall behavior and the evolution of the interface area of multi-mode RT instabilities.

FIG. 15.

Left: temporal evolution of the interface area for the 3D single- and multi-mode RT instabilities for Re = 3000 and At = 0.5. Right: evolution of the ratio of the kinetic energy to released potential energy for the 3D multi-mode instability.

FIG. 15.

Left: temporal evolution of the interface area for the 3D single- and multi-mode RT instabilities for Re = 3000 and At = 0.5. Right: evolution of the ratio of the kinetic energy to released potential energy for the 3D multi-mode instability.

Close modal

Direct numerical simulations of 2D and 3D, single-mode and multi-mode, incompressible immiscible Rayleigh–Taylor instabilities with low to medium Atwood numbers were successfully performed using a newly developed two-phase solver within Xcompact3d, a high-order finite-difference CFD framework. A phase-field approach based on the Allen–Cahn formulation in conjunction with a CSF surface tension model was implemented to perform the simulations. For the single-mode 2D and 3D instabilities, various combinations of Atwood number (0.2At0.75), Reynolds number (Re=300,1000,3000, and 10 000 for the 2D instabilities, and Re=1000,1024, and 3000 for the 3D instabilities), surface tension force, and initial perturbation amplitude were investigated. To the authors' best knowledge, this is one of the very few studies of RT instabilities using high-order schemes (6th-order compact schemes in conjunction with 5th-order WENO schemes).

The overall topology of the single-mode instabilities, including the spike roll-up and vortex ring formation, secondary roll-up of the 3D instabilities, the turbulent flow patterns as the result of the formation of the secondary vortices at the tail of the spike vortices, and the saddle points of the 3D instabilities, was captured accurately. The transient movements of the coordinates of the spike, bubble, and saddle point were also predicted accurately. It was found that, at low and medium Atwood numbers, unless the surface tension effect is relatively significant, RT instabilities could be safely modeled without considering the surface tension force. Moreover, it was concluded that, in the absence of the surface tension or when its value is very low, the diffusion term of the Allen–Cahn phase-field transport equation could be safely neglected, with a significant positive impact on the performance of the flow solver. The effect of the surface tension was found to be more dominant at relatively high Reynolds numbers. Specifically, it was found to prevent the formation of Kelvin–Helmholtz type instabilities within the bubble region (i.e., bubble perturbation) and also the stem section of the spike in a 2D RT instability with Re = 10 000. For both the 2D and 3D single-mode RT instabilities, the Reynolds number was found to have significant effects on the topology of the transient vortices and the chaotic region created due to the spike rolling up while it did not noticeably influence the absolute penetration lengths of the spike and bubble (and the saddle point in 3D RT instabilities) features. On the other hand, the latter absolute penetration lengths were greatly affected by increasing the Atwood number from At = 0.2 to At = 0.75 in both 2D and 3D single-mode RT instabilities. Specifically, a higher Atwood number, resulted in a longer absolute penetration lengths of the spike and the saddle point and a lower absolute penetration length of the bubble. However, the spike exhibited a significant nonlinear growth with respect to the Atwood number while the rate of change of the bubble coordinate was almost linear. It was shown that with similar Atwood and Reynolds numbers a 3D RT instability exhibited larger absolute penetration lengths of the spike and bubble compared to a 2D RT instability (by 40%50%). For both 2D and 3D single-mode instabilities, the spike and bubble reaccelerated after reaching a temporary plateau. Such behavior was observed for all Atwood numbers studied here (except for the spike of the 3D cases with At0.5 which showed an increasing velocity trend for the majority of the simulation time). Confirming some earlier observations, it was concluded that such flow reacceleration was triggered by the reduction of the friction drag due to the formation of the spike vortices and also the formation of a momentum jet (due to the secondary vortices at the tail of the spike vortices) traveling upward within the bubble region. It was concluded that a 2D RT instability could exhibit bubble (and spike) reacceleration (i.e., higher nonlinearity) with lower initial perturbation amplitudes compared to a 3D RT instability with similar Atwood and Reynolds. A relationship between the vertical distance of the bubble and spike |ysb*| and the Atwood number was proposed for the 2D and 3D single-mode RT instabilities as |ysb*|=η[θ1exp(ψkAtP)+θ2exp(ψkAtP)], with θ1=0.653,θ2=0.362, P = 2, and ψ=0.25/At|0.5 for the 2D single-mode RT instabilities. It was shown that for the 3D single-mode instabilities, the latter correlation could be used with similar θ1, θ2, and ψ values to those used for the 2D instabilities but with P = 4 and different adjustment coefficients η.

The Reynolds number was found to have a noticeable influence on the area of the interface between the heavy and light fluids in a 3D single-mode instability. The interface grew exponentially with time for all cases; however, a higher Reynolds number resulted in a noticeably larger surface area at any specific time after the initial growth of the instability. On the other hand, at a specific tAt, a lower Atwood number exhibited a slightly larger interface surface area, more evident in cases with a higher Reynolds number (Re = 3000). This is suggested to be related to the formation of relatively larger vortical structures and consequently large turbulent features at the tail of the spike roll-up at smaller Atwood numbers as mentioned earlier. A direct comparison between a 3D multi-mode RT instability at At = 0.5 and Re = 3000 with comparable 3D single-mode instabilities showed that the multi-mode instability exhibited an exponential growth rate for the interface surface area similarly to those of the single-mode instabilities. However, unlike the single-mode instabilities, the growth rate of the interface area eventually slows down because of the individual single-mode building block RT instabilities of the multi-mode RT instability. Specifically, the individual instabilities started merging and collapsing after an initial growth and formed relatively bulkier but less wrinkled irregular instabilities. The three different mesh resolutions exhibit a similar interface growth rate prior to the formation of the late-time elongated irregular instabilities and similar values for the ratio of the kinetic energy to released potential energy.

The authors would like to acknowledge the financial support and computational time on ARCHER2, the UK supercomputing facility, from the UK Turbulence Consortium (UKTC) under the EPSRC Grant No. EP/R029326/1. The authors also thank the financial support from the UK Supercomputing Facility under Grant Nos. ARCHER-eCSE13-3 and ARCHER2-eCSE01-6. The use of Imperial College London's High Performance Computing Facility and associated support services to conduct part of simulations presented in this paper is gratefully acknowledged.

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.4722736, Ref. 70.

1. Comparison with reference data

The ability of the two-phase solver to accurately capture a 2D RT instability benchmark with At = 0.5 and Re = 3000 as previously studied by various researchers with different surface-capturing methods, including the phase-field method, among others, Tryggvason,2 Liang et al.,15 Huang et al.,17,18 Ding et al.,26 Chiu and Lin.36 An ultra-low surface tension value of σ=1×105 is used to mimic experimental circumstances where a liquid surfactant is added to reduce the surface tension.15 Note that this 2D RT instability has been widely modeled without the surface tension force, for example, by Huang et al.17 and Ding et al.26 Similarly to the simulations discussed in the main body of this paper, the time step is determined by Δt=2.5×104/At for the simulations discussed in  Appendixes A and  B.

Figure 16 (left) shows the temporal evolution of the interface on a 257 × 1025 mesh (case 5 of Table I). The complex dynamics of the interface compares well with those reported previously for a similar Atwood number. The coordinates of the spike and bubble of the RT instability vs tAt are presented in Fig. 16 (right) with a direct comparison with the data of Huang et al.17 and a very good agreement is found between the present DNS results and reference data.

FIG. 16.

Left: evolution of the interface of the 2D single-mode RT instability with At = 0.5 with 257 × 1025 mesh nodes. Right: direct comparison between the present DNS study and the reference data of Huang et al.17 for the vertical coordinates of the tip of the falling (spike) and rising (bubble) fluid in the 2D single-mode RT instability with At = 0.5 and with 257 × 1025 mesh nodes.

FIG. 16.

Left: evolution of the interface of the 2D single-mode RT instability with At = 0.5 with 257 × 1025 mesh nodes. Right: direct comparison between the present DNS study and the reference data of Huang et al.17 for the vertical coordinates of the tip of the falling (spike) and rising (bubble) fluid in the 2D single-mode RT instability with At = 0.5 and with 257 × 1025 mesh nodes.

Close modal
2. Mesh study

The ability of the two-phase solver to capture the finest features of a 2D RT instability at different spatial resolutions is examined by performing a mesh independence study. The results provided in this subsection belong to the cases with At = 0.5 and Re = 3000. In Fig. 17, the interface (ϕ=0.5) at tAt=2, tAt=2.5, and tAt=3 is plotted for three mesh resolutions: 65 × 257, 129 × 513, and 257 × 1025 mesh nodes. Except for the 65 × 257 simulation, which under-predicts the peak locations of the spike and bubble features, the main difference is in the details of the roll-ups and their associated chaotic (turbulent) features. Also, from Fig. 17, it is evident that the differences between the different mesh resolutions grow as the solution evolves and becomes more apparent at later times. Coordinates predicted for the 129 × 513 and 257 × 1025 simulations are almost on top of each other while those of the 65 × 257 simulation are slightly under-predicted (not shown here). It is clear that for the present instability, 129 × 513 mesh nodes can successfully capture the fine details of the RT instability. High-order schemes such as the ones used in the present study are able to capture more small scales at a given resolution than low-order schemes.41 Huang et al.17,18 also used a comparable mesh resolution for a similar RT instability. Furthermore, Mirjalili et al.27 showed that a reasonable solution can be obtained on a 128 × 512 mesh for a RT instability with At0.75. Therefore, a mesh resolution of 129 × 513 is used to evaluate the effects of the Atwood and Reynolds numbers and the surface tension force in the present study.

FIG. 17.

Direct comparison for the location of the interface (ϕ=0.5) for the 2D single-mode RT instabilities with At = 0.5, Re = 3000, and 65 × 257 mesh nodes (green), 129 × 513 mesh nodes (blue), and 257 × 1025 mesh nodes (red). From left to right the snapshots correspond to tAt=2,2.5, and 3.

FIG. 17.

Direct comparison for the location of the interface (ϕ=0.5) for the 2D single-mode RT instabilities with At = 0.5, Re = 3000, and 65 × 257 mesh nodes (green), 129 × 513 mesh nodes (blue), and 257 × 1025 mesh nodes (red). From left to right the snapshots correspond to tAt=2,2.5, and 3.

Close modal

The ability of the two-phase solver to model 3D single-mode RT instabilities is evaluated in this section. Specifically, single-mode instabilities with A = 0.05, Re = 1024, and At = 0.2 and 0.5 (cases 1 and 4 of Table II) are modeled and compared against the results of similar RT instabilities reported by He et al.,66 Lee and Kim.67 The coordinates of the spike, bubble, and saddle point of the 3D RT instabilities (A = 0.05, Re = 1024) with At = 0.5 and 0.2 vs tAt are presented in Fig. 18. The latter figure provides a direct comparison with the simulation results of He et al.66 and Lee and Kim.67 A very good agreement is observed between the present DNS results and the reference data for the coordinates of the spike, bubble, and saddle point.

FIG. 18.

Direct comparison between the present DNS study and the reference data of He et al.,66 Lee and Kim67 for the temporal evolution of the vertical coordinates of the tip of the falling (spike) and rising (bubble) fluid and vertical coordinate of the saddle point in the 3D single-mode RRT instability for Re = 1024, and At = 0.5 (left) and At = 0.2 (right).

FIG. 18.

Direct comparison between the present DNS study and the reference data of He et al.,66 Lee and Kim67 for the temporal evolution of the vertical coordinates of the tip of the falling (spike) and rising (bubble) fluid and vertical coordinate of the saddle point in the 3D single-mode RRT instability for Re = 1024, and At = 0.5 (left) and At = 0.2 (right).

Close modal
1.
D. H.
Sharp
, “
An overview of Rayleigh-Taylor instability
,”
Physica D: Nonlinear Phenomena
,
12
(
1–3
),
3
18
(
1984
).
2.
G.
Tryggvason
, “
Numerical simulations of the Rayleigh–Taylor instability
,”
J. Comput. Phys.
75
,
253
282
(
1988
).
3.
G.
Boffetta
and
A.
Mazzino
, “
Incompressible Rayleigh–Taylor turbulence
,”
Annu. Rev. Fluid Mech.
49
,
119
143
(
2017
).
4.
R.
Lord
, “
Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density
,” in
Proceedings of the London Mathematical Society
(
1882
), pp.
170
177
.
5.
G. I.
Taylor
, “
The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. I
,”
Proc. R. Soc. London, Ser. A
201
,
192
196
(
1950
).
6.
J.
Bell
,
M.
Day
,
C.
Rendleman
,
S.
Woosley
, and
M.
Zingale
, “
Direct numerical simulations of type Ia supernovae flames. II. The Rayleigh–Taylor instability
,”
Astrophys. J.
608
,
883
(
2004
).
7.
G.
Zulauf
,
J.
Zulauf
,
A.
Thiessen
, and
E.
Hattingen
, “
Formation of dome-in-dome structures: Results from experimental studies and comparison with natural examples
,”
J. Struct. Geol.
118
,
324
339
(
2019
).
8.
S.
Ravichandran
,
E.
Meiburg
, and
R.
Govindarajan
, “
Mammatus cloud formation by settling and evaporation
,”
J. Fluid Mech.
899
,
A27
(
2020
).
9.
G.
Zhu
,
P.
Shi
,
Z.
Yang
,
J.
Zheng
,
M.
Luo
,
J.
Ying
, and
X.
Sun
, “
A new method to suppress the Rayleigh–Taylor instability in a linear device
,”
Phys. Plasmas
26
,
042107
(
2019
).
10.
J.
Marthelot
,
E.
Strong
,
P. M.
Reis
, and
P.-T.
Brun
, “
Designing soft materials with interfacial instabilities in liquid films
,”
Nat. Commun.
9
,
1
7
(
2018
).
11.
Y.
Zhou
, “
Rayleigh–Taylor and Richtmyer–Meshkov instability induced flow, turbulence, and mixing. I
,”
Phys. Rep.
720–722
,
1
136
(
2017
).
12.
Y.
Zhou
, “
Rayleigh–Taylor and Richtmyer–Meshkov instability induced flow, turbulence, and mixing. II
,”
Phys. Rep.
723
,
1
160
(
2017
).
13.
Z.-X.
Hu
,
Y.-S.
Zhang
,
B.
Tian
,
Z.
He
, and
L.
Li
, “
Effect of viscosity on two-dimensional single-mode Rayleigh–Taylor instability during and after the reacceleration stage
,”
Phys. Fluids
31
,
104108
(
2019
).
14.
T.
Luo
,
J.
Wang
,
C.
Xie
,
M.
Wan
, and
S.
Chen
, “
Effects of compressibility and Atwood number on the single-mode Rayleigh–Taylor instability
,”
Phys. Fluids
32
,
012110
(
2020
).
15.
H.
Liang
,
X.
Hu
,
X.
Huang
, and
J.
Xu
, “
Direct numerical simulations of multi-mode immiscible Rayleigh–Taylor instability with high Reynolds numbers
,”
Phys. Fluids
31
,
112104
(
2019
).
16.
F.
Ren
,
B.
Song
,
M. C.
Sukop
, and
H.
Hu
, “
Improved lattice Boltzmann modeling of binary flow based on the conservative Allen–Cahn equation
,”
Phys. Rev. E
94
,
023311
(
2016
).
17.
Z.
Huang
,
G.
Lin
, and
A. M.
Ardekani
, “
Consistent, essentially conservative and balanced-force Phase-field method to model incompressible two-phase flows
,”
J. Comput. Phys.
406
,
109192
(
2020
).
18.
Z.
Huang
,
G.
Lin
, and
A. M.
Ardekani
, “
A mixed upwind/central WENO scheme for incompressible two-phase flows
,”
J. Comput. Phys.
387
,
455
480
(
2019
).
19.
H.
Liang
,
Q.
Li
,
B.
Shi
, and
Z.
Chai
, “
Lattice Boltzmann simulation of three-dimensional Rayleigh–Taylor instability
,”
Phys. Rev. E
93
,
033113
(
2016
).
20.
V.
Goncharov
, “
Analytical model of nonlinear, single-mode, classical Rayleigh–Taylor instability at arbitrary Atwood numbers
,”
Phys. Rev. Lett.
88
,
134502
(
2002
).
21.
J.
Waddell
,
C.
Niederhaus
, and
J. W.
Jacobs
, “
Experimental study of Rayleigh–Taylor instability: Low Atwood number liquid systems with single-mode initial perturbations
,”
Phys. Fluids
13
,
1263
1273
(
2001
).
22.
J.
Wilkinson
and
J. W.
Jacobs
, “
Experimental study of the single-mode three-dimensional Rayleigh–Taylor instability
,”
Phys. Fluids
19
,
124102
(
2007
).
23.
M. J.
Andrews
and
S. B.
Dalziel
, “
Small Atwood number Rayleigh–Taylor experiments
,”
Philos. Trans. R. Soc., A
368
,
1663
1679
(
2010
).
24.
I.
Yilmaz
, “
Analysis of Rayleigh–Taylor instability at high Atwood numbers using fully implicit, non-dissipative, energy-conserving large eddy simulation algorithm
,”
Phys. Fluids
32
,
054101
(
2020
).
25.
V. R.
Gopala
and
B. G.
van Wachem
, “
Volume of fluid methods for immiscible-fluid and free-surface flows
,”
Chem. Eng. J.
141
,
204
221
(
2008
).
26.
H.
Ding
,
P. D.
Spelt
, and
C.
Shu
, “
Diffuse interface model for incompressible two-phase flows with large density ratios
,”
J. Comput. Phys.
226
,
2078
2095
(
2007
).
27.
S.
Mirjalili
,
C. B.
Ivey
, and
A.
Mani
, “
Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows
,”
Int. J. Multiphase Flow
116
,
221
238
(
2019
).
28.
L.
Zheng
,
S.
Zheng
, and
Q.
Zhai
, “
Multiphase flows of N immiscible incompressible fluids: Conservative Allen–Cahn equation and lattice Boltzmann equation method
,”
Phys. Rev. E
101
,
013305
(
2020
).
29.
G.
Malamud
,
L.
Elgin
,
T.
Handy
,
C.
Huntington
,
R.
Drake
,
D.
Shvarts
,
A.
Shimony
, and
C.
Kuranz
, “
Design of a single-mode Rayleigh–Taylor instability experiment in the highly nonlinear regime
,”
High Energy Density Phys.
32
,
18
30
(
2019
).
30.
P.
Ramaprabhu
,
G.
Dimonte
,
Y.-N.
Young
,
A.
Calder
, and
B.
Fryxell
, “
Limits of the potential flow approach to the single-mode Rayleigh–Taylor problem
,”
Phys. Rev. E
74
,
066308
(
2006
).
31.
P.
Ramaprabhu
,
G.
Dimonte
,
P.
Woodward
,
C.
Fryer
,
G.
Rockefeller
,
K.
Muthuraman
,
P.-H.
Lin
, and
J.
Jayaraj
, “
The late-time dynamics of the single-mode Rayleigh–Taylor instability
,”
Phys. Fluids
24
,
074107
(
2012
).
32.
R.
Zanella
,
G.
Tegze
,
R. L.
Tellier
, and
H.
Henry
, “
Two-and three-dimensional simulations of Rayleigh–Taylor instabilities using a coupled Cahn–Hilliard/Navier–Stokes model
,”
Phys. Fluids
32
,
124115
(
2020
).
33.
J.
Ding
,
P.
Sun
,
S.
Huang
, and
X.
Luo
, “
Single-and dual-mode Rayleigh–Taylor instability at microscopic scale
,”
Phys. Fluids
33
,
042102
(
2021
).
34.
D. L.
Youngs
, “
Rayleigh–Taylor mixing: Direct numerical simulation and implicit large eddy simulation
,”
Phys. Scr.
92
,
074006
(
2017
).
35.
D.
Livescu
,
T.
Wei
, and
P.
Brady
, “
Rayleigh–Taylor instability with gravity reversal
,”
Phys. D
417
,
132832
(
2021
).
36.
P.-H.
Chiu
and
Y.-T.
Lin
, “
A conservative phase field method for solving incompressible two-phase flows
,”
J. Comput. Phys.
230
,
185
204
(
2011
).
37.
S.
Mirjalili
,
S. S.
Jain
, and
M.
Dodd
, “
Interface-capturing methods for two-phase flows: An overview and recent developments
,”
Cent. Turbul. Res. Annu. Res. Briefs
2017
,
117
135
; available at https://web.stanford.edu/group/ctr/ResBriefs/2017/09_Dodd_117_135.pdf.
38.
P.
Bartholomew
,
G.
Deskos
,
R. A.
Frantz
,
F. N.
Schuch
,
E.
Lamballais
, and
S.
Laizet
, “
Xcompact3D: An open-source framework for solving turbulence problems on a Cartesian mesh
,”
SoftwareX
12
,
100550
(
2020
).
39.
P.
Bartholomew
and
S.
Laizet
, “
A new highly scalable, high-order accurate framework for variable-density flows: Application to non-Boussinesq gravity currents
,”
Comput. Phys. Commun.
242
,
83
94
(
2019
).
40.
S. K.
Lele
, “
Compact finite difference schemes with spectral-like resolution
,”
J. Comput. Phys.
103
,
16
42
(
1992
).
41.
S.
Laizet
and
E.
Lamballais
, “
High-order compact schemes for incompressible flows: A simple and efficient method with quasi-spectral accuracy
,”
J. Comput. Phys.
228
,
5989
6015
(
2009
).
42.
M. S.
Dodd
and
A.
Ferrante
, “
A fast pressure-correction method for incompressible two-fluid flows
,”
J. Comput. Phys.
273
,
416
434
(
2014
).
43.
M.
Ishii
, “
Two-fluid model for two-phase flow
,”
Multiphase Sci. Technol.
5
,
1
(
1990
).
44.
J. U.
Brackbill
,
D. B.
Kothe
, and
C.
Zemach
, “
A continuum method for modeling surface tension
,”
J. Comput. Phys.
100
,
335
354
(
1992
).
45.
X.-D.
Liu
,
R. P.
Fedkiw
, and
M.
Kang
, “
A boundary condition capturing method for Poisson's equation on irregular domains
,”
J. Comput. Phys.
160
,
151
178
(
2000
).
46.
T. E.
Tezduyar
, “
Interface-tracking and interface-capturing techniques for finite element computation of moving boundaries and interfaces
,”
Comput. Methods Appl. Mech. Eng.
195
,
2983
3000
(
2006
).
47.
M.
Sussman
,
P.
Smereka
, and
S.
Osher
, “
A level set approach for computing solutions to incompressible two-phase flow
,”
J. Comput. Phys.
114
,
146
159
(
1994
).
48.
C. W.
Hirt
and
B. D.
Nichols
, “
Volume of fluid (VOF) method for the dynamics of free boundaries
,”
J. Comput. Phys.
39
,
201
225
(
1981
).
49.
F.
Xiao
,
Y.
Honma
, and
T.
Kono
, “
A simple algebraic interface capturing scheme using hyperbolic tangent function
,”
Int. J. Numer. Methods Fluids
48
,
1023
1040
(
2005
).
50.
D. M.
Anderson
,
G. B.
McFadden
, and
A. A.
Wheeler
, “
Diffuse-interface methods in fluid mechanics
,”
Annu. Rev. Fluid Mech.
30
,
139
165
(
1998
).
51.
S.
Mirjalili
and
A.
Mani
, “
Consistent, energy-conserving momentum transport for simulations of two-phase flows using the phase field equations
,”
J. Comput. Phys.
426
,
109918
(
2021
).
52.
S.
Dong
and
J.
Shen
, “
A time-stepping scheme involving constant coefficient matrices for phase-field simulations of two-phase incompressible flows with large density ratios
,”
J. Comput. Phys.
231
,
5788
5804
(
2012
).
53.
E.
Olsson
and
G.
Kreiss
, “
A conservative level set method for two phase flow
,”
J. Comput. Phys.
210
,
225
246
(
2005
).
54.
G.-S.
Jiang
and
C.-W.
Shu
, “
Efficient implementation of weighted ENO schemes
,”
J. Comput. Phys.
126
,
202
228
(
1996
).
55.
R. J.
McSherry
,
K. V.
Chua
, and
T.
Stoesser
, “
Large eddy simulation of free-surface flows
,”
J. Hydrodyn.
29
,
1
12
(
2017
).
56.
M.
Sussman
and
P.
Smereka
, “
Axisymmetric free boundary problems
,”
J. Fluid Mech.
341
,
269
294
(
1997
).
57.
J.
Zhang
and
T. L.
Jackson
, “
A high-order incompressible flow solver with WENO
,”
J. Comput. Phys.
228
,
2426
2442
(
2009
).
58.
F.
Nicoud
, “
Numerical study of a channel flow with variable properties
,”
Cent. Turbul. Res. Annu. Res. Briefs
1988
,
289
(
1999
); available at https://web.stanford.edu/group/ctr/ResBriefs98/nicoud.pdf.
59.
F.
Nicoud
, “
Conservative high-order finite-difference schemes for low-Mach number flows
,”
J. Comput. Phys.
158
,
71
97
(
2000
).
60.
E.
Motheau
and
J.
Abraham
, “
A high-order numerical algorithm for DNS of low-Mach-number reactive flows with detailed chemistry and quasi-spectral accuracy
,”
J. Comput. Phys.
313
,
430
454
(
2016
).
61.
N.
Li
and
S.
Laizet
, “
2DECOMP&FFT—A highly scalable 2D decomposition library and FFT interface
,” in
Cray User Group 2010 Conference
(
2010
), pp.
1
13
.
62.
S.
Laizet
and
N.
Li
, “
Incompact3d: A powerful tool to tackle turbulence problems with up to O(105) computational cores
,”
Int. J. Numer. Methods Fluids
67
,
1735
1757
(
2011
).
63.
J. M.
Stone
and
T.
Gardiner
, “
Nonlinear evolution of the magnetohydrodynamic Rayleigh–Taylor instability
,”
Phys. Fluids
19
,
094104
(
2007
).
64.
T.
Zhang
,
J.
Wu
, and
X.
Lin
, “
An interface-compressed diffuse interface method and its application for multiphase flows
,”
Phys. Fluids
31
,
122102
(
2019
).
65.
X.
Bian
,
H.
Aluie
,
D.
Zhao
,
H.
Zhang
, and
D.
Livescu
, “
Revisiting the late-time growth of single-mode Rayleigh–Taylor instability and the role of vorticity
,”
Phys. D
403
,
132250
(
2020
).
66.
X.
He
,
R.
Zhang
,
S.
Chen
, and
G. D.
Doolen
, “
On the three-dimensional Rayleigh–Taylor instability
,”
Phys. Fluids
11
,
1143
1152
(
1999
).
67.
H. G.
Lee
and
J.
Kim
, “
Numerical simulation of the three-dimensional Rayleigh–Taylor instability
,”
Comput. Math. Appl.
66
,
1466
1474
(
2013
).
68.
G. C.
Burton
, “
Study of ultrahigh Atwood-number Rayleigh–Taylor mixing dynamics using the nonlinear large-eddy simulation method
,”
Phys. Fluids
23
,
045106
(
2011
).
69.
G.
Dimonte
,
D.
Youngs
,
A.
Dimits
,
S.
Weber
,
M.
Marinak
,
S.
Wunsch
,
C.
Garasi
,
A.
Robinson
,
M.
Andrews
,
P.
Ramaprabhu
 et al, “
A comparative study of the turbulent Rayleigh–Taylor instability using high-resolution three-dimensional numerical simulations: The Alpha-Group collaboration
,”
Phys. Fluids
16
,
1668
1693
(
2004
).
70.
A.
Hamzehloo
,
P.
Bartholomew
, and
S.
Laizet
(
2021
). “Direct numerical simulations of incompressible Rayleigh–Taylor instabilities at low and medium Atwood numbers,”
Zenodo
.