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.

## I. INTRODUCTION

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 Rayleigh^{4} and Taylor^{5} 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 supernovae^{6}), geology (initiation and evolution of polydiapirs^{7}), atmospheric and climate science (formation of the mammatus clouds^{8}), inertial confinement fusion (premature fuel mixing and reduction in the heating efficacy^{9}), and designing and fabricating soft materials (structured materials in thin polymeric liquid films^{10}). 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\lambda $,^{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

where $\rho \u0303$ 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=U\u0303L\u0303\rho r\u0303/\mu r\u0303$, where $U\u0303,\u2009L\u0303,\u2009\rho r\u0303$, and $\mu r\u0303$ 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 ($At\u22640.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,\u20090.74$, and $1.34g\u0303$ in excess of the earth's gravitational acceleration $g\u0303$ 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 Jacobs^{22} 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 $At\u22480.15$ and $At\u22480.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.05\u2264At\u22640.9$. Later, Ramaprabhu *et al.*^{31} showed, based on 3D single-mode miscible RT instabilities with $0.1\u2264At\u22640.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 ($10\u2264Re\u22645000$) 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 $100\u2264Re\u226410\u2009000$. 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 ($At\u22480.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} Yilmaz^{24} used a large eddy simulation approach to computationally study multi-mode immiscible RT instabilities with $0.5\u2264At\u22640.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.2\u2264At\u22640.75$), Reynolds number ($300\u2264Re\u226410\u2009000$), 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 method^{27,36,37} was developed and used to perform DNS of the RT instabilities. This new solver is implemented within the Xcompact3d framework^{38} 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.

## II. COMPUTATIONAL METHODOLOGY

### A. Governing equations

#### 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

where $ui(i=1,2,3)$ denotes the velocity component in the *i*th 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

*i*th direction are represented by $fi\sigma $ and

*g*, respectively, while

_{i}*τ*is the full viscous stress tensor defined as

_{ij}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=U\u03032/(g\u0303L\u0303)$. As mentioned earlier, only variables accented with the tilde sign (∼) are dimensional and $U\u0303,\u2009L\u0303,\u2009\rho r\u0303,\u2009\mu r\u0303$, and $g\u0303$ 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 $U\u0303=g\u0303L\u0303$, which when inserted into the Reynolds number yields the expression $Re=\rho \u0303rL\u03033g\u0303/\mu \u0303r$ which is used to define the Reynolds number in the present study. The density of the heavy fluid is chosen as the reference density $\rho r\u0303=\rho h\u0303$, and both fluids are assumed to have a similar viscosity ($\mu r\u0303=\mu h\u0303=\mu l\u0303$) 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 domain^{18} and the boundary conditions at the fluid interface are implicitly imposed using the continuous surface force (CSF) method^{44} or the ghost fluid method.^{45} The location of the interface in the one-fluid approach can either be tracked using front-tracking methods^{46} 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

where $\varphi $ is the phase-field function, *ε* is a measure of the interface thickness, and *γ* is a positive coefficient. The phase-field function takes values $0\u2264\varphi \u22641$ with $\varphi =0,1$ representing regions of solely heavy or light fluid, respectively; it follows that the local averaged density and viscosity are given by $\rho =\varphi \rho l+(1\u2212\varphi )\rho h$ and $\mu =\varphi \mu l+(1\u2212\varphi )\mu h$, respectively. Here, following Mirjalili *et al.*,^{27} $\epsilon =\Delta $, where Δ is the spatial resolution and $\gamma =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 Kothe^{44} and followed by, for example, Refs. 36, 37, and 53 is implemented in the present work given as

where the surface tension coefficient is defined as $\sigma =1/We$ in which *We* is the Weber number defined as $We=\rho 1\u0303U\u03032L\u0303/\sigma \u0303$.

### B. Numerical methods

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 accuracy^{41} is employed. Specifically, for an arbitrary variable *f* the $1st$ and $2nd$ spatial derivatives at the *j*th mesh node in the *x _{i}* direction are approximated as follows:

with $\alpha =1/3$, *a* = 14/9, and *b* = 1/9, and

with $\alpha =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=(j\u22121)\Delta xi$ for $1\u2264j\u2264nxi$ where $nxi$ denotes the number of mesh nodes in the *x _{i}* 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) scheme^{54} 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\u0302$ at the mesh nodes *j* can be estimated (in a specific coordinate direction) using an upwind WENO reconstruction as follows:^{54,57}

where the nonlinear weights *ω _{k}* on the

*k*th stencil candidate (smoothness indicator,

*β*) defined as

_{k}^{18,54,57}

with

where $\epsilon w=10\u22126$ to avoid zero at the denominator. For a five-point WENO reconstruction with *k* = 3 and a $5th$-order accuracy, $\gamma 1=0.3,\u2009\gamma 2=0.6$, and $\gamma 3=0.1$, and for a positive velocity, the reconstruction is given as^{18,54,57}

with the smoothness indicators given as

It should be noted that the gradients of the phase function $\u2202\varphi /\u2202xi$ 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 form^{39,58–60}

where $ui\u22c6$ 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 $\rho \u0302$ appearing in Eq. (14) has the effect of an under-relaxation factor with lower values stabilizing the iterative solution process. In this work, $\rho \u0302=min(\rho )$, 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

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}

### C. Problem definition

#### 1. 2D single-mode

All 2D RT instabilities studied here are modeled in a rectangular domain of size $[0,d]\xd7[0,4d]$ where *d* is the reference length and the interface is initially located at *y*_{0}, defined as

Equation (16) represents a perturbation with wave number $k=2\pi /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

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 $\Delta t=2.5\xd710\u22124/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).

Case . | At
. | $\rho h\rho l$ . | Re
. | $\sigma =1We$ . | Mesh size . |
---|---|---|---|---|---|

1 | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 65 × 257 |

2 | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 129 × 513 |

3 | 0.5 | 3 | 3000 | 0.0 | 129 × 513 |

4^{a} | 0.5 | 3 | 3000 | 0.0 | 129 × 513 |

5 | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 257 × 1025 |

6 | 0.5 | 3 | 1000 | $1\xd710\u22125$ | 129 × 513 |

7 | 0.5 | 3 | 300 | $1\xd710\u22125$ | 129 × 513 |

8 | 0.5 | 3 | 10 000 | $1\xd710\u22125$ | 129 × 513 |

9 | 0.2 | 1.5 | 3000 | $1\xd710\u22125$ | 129 × 513 |

10 | 0.333 | 2 | 3000 | $1\xd710\u22125$ | 129 × 513 |

11 | 0.666 | 5 | 3000 | $1\xd710\u22125$ | 129 × 513 |

12 | 0.75 | 7 | 3000 | $1\xd710\u22125$ | 129 × 513 |

13 | 0.5 | 3 | 3000 | $1\xd710\u22124$ | 129 × 513 |

14 | 0.5 | 3 | 3000 | $5\xd710\u22124$ | 129 × 513 |

15 | 0.5 | 3 | 10 000 | $5\xd710\u22124$ | 129 × 513 |

Case . | At
. | $\rho h\rho l$ . | Re
. | $\sigma =1We$ . | Mesh size . |
---|---|---|---|---|---|

1 | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 65 × 257 |

2 | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 129 × 513 |

3 | 0.5 | 3 | 3000 | 0.0 | 129 × 513 |

4^{a} | 0.5 | 3 | 3000 | 0.0 | 129 × 513 |

5 | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 257 × 1025 |

6 | 0.5 | 3 | 1000 | $1\xd710\u22125$ | 129 × 513 |

7 | 0.5 | 3 | 300 | $1\xd710\u22125$ | 129 × 513 |

8 | 0.5 | 3 | 10 000 | $1\xd710\u22125$ | 129 × 513 |

9 | 0.2 | 1.5 | 3000 | $1\xd710\u22125$ | 129 × 513 |

10 | 0.333 | 2 | 3000 | $1\xd710\u22125$ | 129 × 513 |

11 | 0.666 | 5 | 3000 | $1\xd710\u22125$ | 129 × 513 |

12 | 0.75 | 7 | 3000 | $1\xd710\u22125$ | 129 × 513 |

13 | 0.5 | 3 | 3000 | $1\xd710\u22124$ | 129 × 513 |

14 | 0.5 | 3 | 3000 | $5\xd710\u22124$ | 129 × 513 |

15 | 0.5 | 3 | 10 000 | $5\xd710\u22124$ | 129 × 513 |

^{a}

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

Case . | Mode . | At
. | $\rho h\rho l$ . | Re
. | $\sigma =1We$ . | A
. | Domain size $(\xd7d)$ . | Mesh size . |
---|---|---|---|---|---|---|---|---|

1 | Single | 0.5 | 3 | 1024 | $1\xd710\u22125$ | 0.05 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

2 | Single | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.05 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

3 | Single | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

4 | Single | 0.2 | 1.5 | 1024 | $1\xd710\u22125$ | 0.05 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

5 | Single | 0.2 | 1.5 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

6 | Single | 0.333 | 2 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

7 | Single | 0.5 | 3 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

8 | Single | 0.666 | 5 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

9 | Single | 0.75 | 7 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

10 | Single | 0.2 | 1.5 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

11 | Single | 0.333 | 2 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

12 | Single | 0.666 | 5 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

13 | Single | 0.75 | 7 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

14 | Single | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.005 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

15 | Multi | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.005 | $\pi 2\xd73\pi \xd7\pi 2$ | $301\xd71801\xd7301$ |

16 | Multi | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.005 | $\pi 2\xd73\pi \xd7\pi 2$ | $401\xd72401\xd7401$ |

17 | Multi | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.005 | $\pi 2\xd73\pi \xd7\pi 2$ | $501\xd73001\xd7501$ |

Case . | Mode . | At
. | $\rho h\rho l$ . | Re
. | $\sigma =1We$ . | A
. | Domain size $(\xd7d)$ . | Mesh size . |
---|---|---|---|---|---|---|---|---|

1 | Single | 0.5 | 3 | 1024 | $1\xd710\u22125$ | 0.05 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

2 | Single | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.05 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

3 | Single | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

4 | Single | 0.2 | 1.5 | 1024 | $1\xd710\u22125$ | 0.05 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

5 | Single | 0.2 | 1.5 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

6 | Single | 0.333 | 2 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

7 | Single | 0.5 | 3 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

8 | Single | 0.666 | 5 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

9 | Single | 0.75 | 7 | 1000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

10 | Single | 0.2 | 1.5 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

11 | Single | 0.333 | 2 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

12 | Single | 0.666 | 5 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

13 | Single | 0.75 | 7 | 3000 | $1\xd710\u22125$ | 0.1 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

14 | Single | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.005 | $1\xd74\xd71$ | $129\xd7513\xd7129$ |

15 | Multi | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.005 | $\pi 2\xd73\pi \xd7\pi 2$ | $301\xd71801\xd7301$ |

16 | Multi | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.005 | $\pi 2\xd73\pi \xd7\pi 2$ | $401\xd72401\xd7401$ |

17 | Multi | 0.5 | 3 | 3000 | $1\xd710\u22125$ | 0.005 | $\pi 2\xd73\pi \xd7\pi 2$ | $501\xd73001\xd7501$ |

#### 2. 3D single-mode and multi-mode

The single-mode 3D RT instabilities studied here are modeled in a $[0,d]\xd7[0,4d]\xd7[0,d]$ domain where the interface is initially located at *y*_{0}, defined as

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\xd7513\xd7129$ 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 $\Delta t=2.5\xd710\u22124/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,\pi d/2]\xd7[0,3\pi d]\xd7[0,\pi d/2]$ discretized with three mesh resolutions, with $301\xd71801\xd7301,\u2009401\xd72401\xd7401$, and $501\xd73001\xd7501$ mesh nodes (cases 15–17 of Table II, respectively). The interface is initially located at $y0=3\pi 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 as^{24}

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}

## III. SINGLE-MODE RAYLEIGH–TAYLOR INSTABILITIES

### A. 2D instabilities

#### 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.

#### 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,\u20091000$, 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*=y\u22122$ 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 $\rho h/\rho 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 ($\rho h/\rho l=30$) the spike roll-up was very limited and for a very high Atwood number ($\rho h/\rho 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.

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 as^{19,21,22}

where $|ys*|$ and $|yb*|$ are the absolute vertical distance of the spike and bubble from $y*=0$ and *a*_{1} and *a*_{2} are the fitting coefficients. For very low surface tension values, the growth rate (velocity) of the initial stage is estimated as $U=At\u2009g\u2009k$. 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 Goncharov^{20} as

where *C* = 3 for 2D RT instabilities and *C* = 1 for 3D instabilities ($k=2\pi $ is the perturbation wavenumber). The spike terminal velocity can be estimated (less accurately) using a similar correlation but replacing $1+At$ with $1\u2212At$ 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 literature^{19,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 $At\u22640.75$ and showed that RT instabilities with *At* = 0.8 and *At* = 1.0 did not exhibit such flow behavior. RT instabilities with $At\u22650.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 $At\u22640.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 ($tAt\u22651.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 $tAt\u22651.5$ is proposed as follows:

where $\theta 1=0.653,\u2009\theta 2=0.362$, *k* is the initial perturbation wavenumber ($2\pi $ here), *P* = 2, and *ψ* is the time coefficient defined as $\psi =0.25/At|0.5$. The adjustment factor $\eta =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 $\psi =1000\Delta 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).

Dimension . | $tAt$ . | Re
. | η
. | θ_{1}
. | θ_{2}
. | ψ
. | P
. |
---|---|---|---|---|---|---|---|

2D | 1.5 | 3000 | 0.47 | 0.653 | 0.362 | 0.353 | 2 |

2D | 2 | 3000 | 0.635 | 0.653 | 0.362 | 0.353 | 2 |

2D | 2.5 | 3000 | 0.805 | 0.653 | 0.362 | 0.353 | 2 |

2D | 3 | 3000 | 1 | 0.653 | 0.362 | 0.353 | 2 |

3D | 2 | 1000 and 3000 | 1.12 | 0.653 | 0.362 | 0.353 | 4 |

3D | 2.25 | 1000 and 3000 | 1.28 | 0.653 | 0.362 | 0.353 | 4 |

Dimension . | $tAt$ . | Re
. | η
. | θ_{1}
. | θ_{2}
. | ψ
. | P
. |
---|---|---|---|---|---|---|---|

2D | 1.5 | 3000 | 0.47 | 0.653 | 0.362 | 0.353 | 2 |

2D | 2 | 3000 | 0.635 | 0.653 | 0.362 | 0.353 | 2 |

2D | 2.5 | 3000 | 0.805 | 0.653 | 0.362 | 0.353 | 2 |

2D | 3 | 3000 | 1 | 0.653 | 0.362 | 0.353 | 2 |

3D | 2 | 1000 and 3000 | 1.12 | 0.653 | 0.362 | 0.353 | 4 |

3D | 2.25 | 1000 and 3000 | 1.28 | 0.653 | 0.362 | 0.353 | 4 |

#### 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 $\sigma =1/We=0,1\xd710\u22125,1\xd710\u22124$, and $5\xd710\u22124$ 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 $\Phi $ of the Cahn–Hilliard approach can be defined as $\Phi =(1/\epsilon )\sigma \alpha \zeta \u2032\u2212\epsilon \sigma \alpha \Delta \varphi $,^{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 $\sigma =5\xd710\u22124$ 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 $\u223c7.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.

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 $\sigma =1\xd710\u22125$ and $5\xd710\u22124$. A relatively high surface tension is found to have a significant effect on the evolution of the RT instability. In the case with $\sigma =5\xd710\u22124$, 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.

### B. 3D instabilities

#### 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 $\varphi =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 Kim^{67} 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.

#### 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,\u20091024$, 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 $\u223c15%$ and $\u223c25%$ 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.

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 $\u223c43%$ and $\u223c51%$ 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.

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 $\u223c10%\u201315%$ 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 $At\u22640.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).

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 $tAt\u22652$ 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 *A _{int}* 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

*A*vs $tAt$ graph. However, the

_{int}*A*vs

_{int}*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.

## IV. 3D MULTI-MODE RT INSTABILITY

In this section, a 3D multi-mode RT instability is investigated thanks to highly resolved DNS with up to 750 ×10^{6} mesh nodes. The mesh resolutions have been chosen based on previous numerical works^{24,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 $\varphi =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\xd71801\xd7301,\u2009401\xd72401\xd7401$, and $501\xd73001\xd7501$ 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$.

To provide a more quantitative comparison, the evolution of the interface area *A _{int}* normalized by the domain cross section area

*A*as a function of $t*=tAt$ for the 3D multi-mode RT instability is compared with those of the 3D single-mode cases with

_{d}*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\u222b\Omega uiui\u2009d\Omega $ and $PE=\u222b\Omega [\rho (xi,0)\u2212\rho (xi,t)]x2\u2009d\Omega $, 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.

## V. CONCLUSIONS

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.2\u2264At\u22640.75$), Reynolds number ($Re=300,\u20091000,\u20093000$, and 10 000 for the 2D instabilities, and $Re=1000,\u20091024$, 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 $\u223c40%\u201350%$). 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 $At\u22650.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*|=\eta [\theta 1\u2009exp\u2009(\psi kAtP)+\theta 2\u2009exp\u2009(\u2212\psi kAtP)]$, with $\theta 1=0.653,\theta 2=0.362$, *P* = 2, and $\psi =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.

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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

### APPENDIX A: VALIDATION AND MESH STUDIES (2D RT INSTABILITIES)

##### 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 $\sigma =1\xd710\u22125$ 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 $\Delta t=2.5\xd710\u22124/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.

##### 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 ($\varphi =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 $At\u22480.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.

### APPENDIX B: VALIDATION (3D SINGLE-MODE RT INSTABILITIES)

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.