Scale-dependent Rayleigh–Taylor dynamics with variable acceleration in a finite-sized domain for three-dimensional flows

Hwang, Hanul; Chan, Wai Hong Ronald; Jain, Suhas S.; Abarzhi, Snezhana I.

doi:10.1063/5.0059898

Three-dimensional Rayleigh–Taylor instability (RTI) with the time-varying acceleration in a finite domain is investigated in a systematic framework. The acceleration magnitude follows a power law in time with an exponent greater than −2. Applying the group theory, the instabilities are demonstrated considering the irreducible representations for observable periodic structures with a square symmetry in the plane normal to the acceleration. We derive the dynamical system and illustrate the universal form of the solutions in the linear and nonlinear regimes. The scale-dependent dynamics are shown to be single scale and multiscale in the two regimes, respectively. For the nonlinear regime solutions, fundamental scales are derived bridging the solutions in the finite- and infinite-sized domains. Special solutions for bubbles and spikes are identified from a one-parameter family of solutions. The effect of domain confinement is that the velocity and curvature decreases and shear increases as the domain size is reduced. The theory provides predictions for the flow field and demonstrates the interfacial behavior of RTI. Our results are in good agreement with the prior studies and also provide new benchmarks for experiments and simulations.

I. INTRODUCTION

The Rayleigh–Taylor instability (RTI) occurs at the interface between two fluids with different densities when accelerated against their density gradient. RTI leads to interpenetration and mixing of the fluids, which is a critical element to the flow evolution. RTI prompts the growth of interface perturbations along the acceleration direction depending on the initial condition and eventually leads to Rayleigh–Taylor (RT) mixing. RT mixing, which manifests distinct characteristics from the classical Kolmogorov turbulence, is driven by the transport of mass, momentum, and energy: RT mixing is a heterogeneous and anisotropic process with nonlocal interactions coupling a wide range of scales.¹ RTI and RT mixing are long-standing paradigm-shifting problems,^2–5 which received widespread attention due to their applications in a broad range of natural, technological, and industrial processes. Examples include inertial confinement fusion (ICF),⁶ supernovae,⁷ and even galactic processes,⁸ which are multiscale phenomena and which also often involve multiphysics. Therefore, understanding the fundamental mechanism and properties of RTI and RT mixing is crucial to further investigate these complex natural phenomena and to control the engineering applications. In the present study, we analyze three-dimensional (3D) RTI dynamics with time-varying acceleration in a finite-sized domain.

As observed in numerous laboratory studies and simulations and in applications,^6,8,9 RTI exhibits several common features in its evolution.¹⁰ It undergoes a series of stages, namely, the scale-dependent linear regime, the scale-dependent nonlinear regime, and the scale-invariant mixing regime. In the linear regime, small amplitudes of an initial perturbation from an equilibrium state initiate RTI. For constant acceleration, perturbations in the linear regime grow exponentially in time. As RTI evolves, the nonlinear regime follows, and the growth of the perturbation amplitude obeys a power-law in time. Thus, the growth rate is relatively slower than in the linear regime.^11,12 At this stage, the interface transforms into a composition of large-scale coherent structures and small-scale irregular structures. The large-scale structures consist of bubbles and spikes, which correspond respectively to the portion of lighter fluid penetrating the heavier fluid and vice versa. Along with the growth of bubbles and spikes, small-scale vortical structures emerge at the interface due to interfacial shear. The resulting enhanced interaction and coupling of scales lead to a transition from the nonlinear regime to the mixing regime where the perturbation amplitude is believed to grow self-similarly with time.^13,14

RTI and RT mixing are nonequilibrium phenomena.¹⁵ Thus, investigation of the subject is challenging due to its inherent complexity. This complicated nature of RTI and RT mixing poses a severe limitation on the study of theory, numerical simulations, and experiments as well as data analysis. From a theoretical standpoint, one should provide a comprehensive explanation incorporating a number of important features of RTI, such as universal properties of asymptotic solutions, multiscale and nonlocal dynamics, and symmetry of RT flows, along with chaotic properties of the mixing. In the computational simulations, an accurate representation of small-scale dynamics near the interface in tandem with the accurate resolution of large-scale bubbles and spikes is challenging. Due to the statistical unsteadiness of interfacial RT mixing, wide spatial and temporal distributions of flow quantities make the experimental measurements extremely difficult in terms of accuracy and repeatability. Therefore, data analysis requires succinct accuracy, adequate precision, and a high acquisition rate of data over a considerable span of scales in space and in time.

Despite the associated difficulties, successful accomplishments are reported in theory, experiments, and simulations.^12,16–18 Particularly, group theory is shown to successfully identify symmetries and invariants of RT dynamics and to find that coherent structures are spatially periodic patterns determined by the criteria of isotropy and structural stability.^{1,7,11,12,19,20}

In the experiments, advanced diagnostics and design facilitate more accurate measurements of RTI and RT mixing. Robey et al.²¹ access high-energy-density conditions in RT-unstable flows, with high values of pressure and temperature. Meshkov²² identifies a concentration jump at the interface, which is an essential structure of the developed mixing zone. Delicate features of interfacial dynamics with complex initial conditions have also been investigated.^14,19 Meshkov and Abarzhi¹ report heterogeneity, anisotropy, and sensitivity of RT mixing in a wide range of setups achieving excellent agreement of theory and experiments. The complicated character of scale coupling, such as the ordered nature and relaminarization of accelerated RT mixing, is revealed.

From a computational standpoint, high-resolution numerical simulations have become feasible, and the cutting-edge data assimilation techniques with error estimation and uncertainty quantification have been developed.^13,23–27 These developments incorporate the numerical modeling of miscible interface with mass transfer,²⁸ studies of the instabilities with magnetohydrodynamics,²⁹ In addition, increasingly accurate interface tracking and capturing methods,^25,27 as well as molecular-level simulations,^7,30,31 have been proven to successfully capture initial condition sensitivities and small-scale dynamics of RT-related phenomena.^27,30,32,33

RTI and RT mixing are a subject of active research. Examples of very recent studies include the effects of viscosity and compressibility on single-mode Rayleigh–Taylor instability,^34,35 RTI at spherical interfaces between viscous fluids,³⁶ dynamic stabilization of the RTI in miscible liquids,³⁷ and numerical modeling of RTI and RT interfacial mixing by using coupled Cahn–Hilliard/Navier–Stokes models.^38,39

Most RTI and RT mixing phenomena involve time-varying acceleration rather than constant or impulsive acceleration.^7,8,40 Moreover, RT experiments and simulations are often conducted in finite-sized domains rather than infinite-sized domains.^6,14,41 In order to describe the problem in more realistic conditions, a systematic way to address the effects of the aforementioned conditions is of crucial importance.¹⁷ More specifically, a better understanding is required on the following: the effect of the domain length on RTI evolution and on the coupling of dominant length scales, the effect of variable acceleration on RTI evolution and on scale coupling, and the interplay of the finite domain sizes and variable acceleration.

In this work, we consider variable accelerations whose magnitudes change with time and whose directions remain the same. The time dependence of the acceleration magnitude is represented by a power-law function. On the side of fundamental science, power-law functions are important to study because they can lead to new invariants and scaling properties of the dynamics.¹⁵ On the side of engineering applications, power-law functions can be used to adjust the acceleration's time dependence in realistic environments and thus ensure practicality of our results.^6,8,15 Examples of RTI induced by an acceleration with power-law time dependence include blast-wave-driven RTI in core-collapse supernovae, RT/Richtmyer–Meshkov (RM)-unstable plasma irregularities in Earth's ionosphere; RTI induced by unsteady shocks in inertial confinement fusion; and RT instabilities in the fossil fuel industry.¹⁶

Accelerations with power-law time dependence can naturally link our present work to earlier studies:¹⁶ The classical RTI corresponds to the case of zero exponent of the acceleration's power-law. The classical Richtmyer–Meshkov (RM) instability, where shock induced impulsive acceleration leads to a constant initial growth rate, corresponds to the case of acceleration exponents smaller than −2.^1–5,15 Furthermore, for integer exponent values the acceleration can have a polynomial time dependence, and for an exponent equal to unity, the acceleration changes linearly with time. In this work, we focus on RTI subject to accelerations whose magnitudes vary as power laws with time and have exponents greater than −2.^1,15,40 We shall consider in other works the cases of RM instability with the exponent of the acceleration's power-law smaller than −2, and the RT-to-RM transition when the exponent of the acceleration's power-law is equal to −2.

The present study investigates RTI dynamics with variable acceleration in a finite-sized domain for 3D flows with square symmetry. We apply and generalize group theory to illustrate the interplay of length scales and the time-varying acceleration. Solutions in the linear and nonlinear regimes are derived for a range of Atwood numbers. A broad range of sizes of the finite domain is explored to discover new invariants and scaling properties as well as the practicality of the theory.

Our work presents the analytical study of 3D highly symmetric RTI with variable acceleration in a finite-sized domain. The novelties of our work are in the systematic and meticulous study on the coupling of the effects of variable acceleration and the finite-sized domain in RTI over a wide range of acceleration parameters, the identification of universal scaling relations and parameterized solutions in linear and nonlinear regimes, and in identifying benchmarks for experiments and simulations. RTI and interfacial mixing are a sample case of nonequilibrium dynamics of interface and interfacial mixing, including the fundamentals of fluid dynamics and applications in nature and technology. Hence, the present work is directly relevant to the scope and the topic of the Collection on “Interfaces and Mixing and Beyond” of Physics of Fluids.

The paper is structured as follows. Section II describes the governing equations and boundary conditions as well as the theory to form the solutions. Also, a power-law-dependent time acceleration is introduced. In Secs. III A and III B, the dynamical system and shear function are derived. Then, solutions in the linear regime and the effect of the domain sizes are discussed in Sec. III C. Section III D highlights a continuous family of one-parameter asymptotic solutions in the nonlinear regime. We propose universal scaling relations in the nonlinear regimes bridging the solutions in a finite-sized domain to those in an infinite-sized domain in Sec. III E. Section III F distinguishes the singularity of the spike solutions. Special solutions for bubbles and spikes are demonstrated in Secs. III G and III H. Based on the identified multiscale character of the dynamics, we explain transitions from scale-dependent dynamics to self-similar mixing for RT bubbles and spikes in Sec. III I. We conclude our results in Sec. IV and discuss new benchmarks for numerical simulations and experiments.

II. METHODS

In this section, the governing equations and boundary conditions are introduced. Then, the fluid potentials satisfying the governing equations and boundary conditions are sought by employing group theory.

A. Boundary value problem setup

The equations governing RTI of inviscid fluids are the conservation laws for mass, momentum, and energy,

(1)

\frac{\partial ρ}{\partial t} + \frac{\partial ρ v_{i}}{\partial x_{i}} = 0,

(1a)

\frac{\partial ρ v_{i}}{\partial t} + \frac{\partial ρ v_{i} v_{j}}{\partial x_{j}} + \frac{\partial P}{\partial x_{i}} = 0,

(1b)

\frac{\partial E}{\partial t} + \frac{\partial (E + P) v_{i}}{\partial x_{i}} = 0,

(1c)

in a spatial coordinate system where

(x_{1}, x_{2}, x_{3}) = (x, y, z)

such that x and y constitute the horizontal plane, which is perpendicular to the acceleration and z represents the vertical dimension. A schematic of the problem setup is depicted in Fig. 1. Here, t denotes the time, and density, velocity, pressure, and energy are introduced as

(ρ, v, P, E)

⁠, respectively, with

v = (v_{1}, v_{2}, v_{3})

⁠. The total energy E is of the form

E = ρ (e + | v |^{2} / 2)

⁠, where e is the specific internal energy. Note that the conservation equations above describe each of the two fluids. The subscripts h and l represent the quantities of the heavy and light fluids, respectively. The effect of surface tension, viscosity, and compressibility can be considered. They are beyond the scope of the present study.

FIG. 1.

View large Download slide

Schematic of the problem. Two fluids with different densities ρ_h and ρ_l are present under time-varying acceleration in the direction opposite to the density gradient. The interface is slightly perturbed from the equilibrium state.

In the laboratory frame, the boundary conditions at the ends of the domain in the vertical direction are

v_{3} (z = Z) = 0, v_{3} (z = - Z) = 0,

(2)

where Z is half the vertical separation of the two ends.

The boundary conditions at the interface are derived from the governing equations in the bulk by introducing a continuous local scalar function, $θ (x, y, z, t)$ ⁠. The continuous scalar function $θ (x, y, z, t)$ is given in the form

θ (x, y, z, t) = {\begin{matrix} θ (x, y, z, t) > 0, z > z^{*} (x, y, t), \\ θ (x, y, z, t) = 0, z = z^{*} (x, y, t), \\ θ (x, y, z, t) < 0, z < z^{*} (x, y, t) . \end{matrix}

(3)

Here, $z^{*} (x, y, t)$ represents the location of the freely evolving unstable interface. Then, the variables $(ρ, v, P, E)$ over the entire flow field are expressed employing the Heaviside step function $H (θ)$ as

(ρ, v, P, E) = H (θ) {(ρ, v, P, E)}_{h} + H (- θ) {(ρ, v, P, E)}_{l} .

(4)

We define an interface unit normal vector $n_{n}$ and an interface unit tangential vector $n_{t}$ using the scalar function $θ (x, y, z, t)$ as

n_{n} = \frac{\nabla θ}{| \nabla θ |}, n_{n} \cdot n_{t} = 0.

(5)

In order to reformulate the governing equations in Eq. (1) to interface boundary conditions, we further introduce a mass flux vector $\tilde{j}$ across the moving interface $z^{*} (x, y, t)$ as

\tilde{j} = ρ (n_{n} \frac{\dot{θ}}{| \nabla θ |} + v),

(6)

where $\nabla θ$ is the gradient of θ and $\dot{θ}$ is its time derivative.

After substituting Eq. (4) into Eq. (1) and manipulating terms using the interface normal and tangential vectors, $n_{n}$ and $n_{t}$ ⁠, and a mass flux $\tilde{j}$ ⁠, Eq. (6) yields the following interface boundary conditions:

(7)

[\tilde{j} \cdot n_{n}] = 0,

(7a)

[(P + \frac{{(\tilde{j} \cdot n_{n})}^{2}}{ρ})] = 0,

(7b)

[(\tilde{j} \cdot n_{n}) \frac{(\tilde{j} \cdot n_{t})}{ρ}] = 0,

(7c)

[(\tilde{j} \cdot n_{n}) (W + \frac{{\tilde{j}}^{2}}{2 ρ^{2}})] = 0,

(7d)

where

W = e + P / ρ

is the specific enthalpy. The square brackets on a variable χ denote the difference of χ across the interface,

[χ] = χ_{h} - χ_{l}

⁠. Note that the conditions in Eq. are derived from the conservation of mass, vertical momentum, horizontal momentum, and energy, respectively, in the inertial frame of reference and are exact.

The immiscibility of the fluids at the interface further simplifies Eq. (7). The mass flux $\tilde{j}$ across the interface in Eq. (7) should satisfy $\tilde{j} \cdot n_{n} |_{θ = 0^{\pm}} = 0$ ⁠. We note that this condition is identical to the contact discontinuity.⁴⁰ Due to this contact discontinuity, the conditions in Eq. (7) may be written as

(8)

[v \cdot n_{n}] = 0,

(8a)

[P] = 0,

(8b)

[v \cdot n_{t}] = arbitrary,

(8c)

[W] = arbitrary .

(8d)

The first two conditions in Eq. represent the continuity of the normal velocity and the pressure, whereas the last two illustrate the discontinuities of the tangential velocity and the enthalpy across the interface, respectively.

The governing equations in Eq. (1), together with the boundary conditions in Eqs. (2) and (8) and with initial conditions, pose the RTI problem.

B. Time-dependent acceleration

The acceleration, $g (t)$ ⁠, is directed from the heavy fluid to the light fluid and varies with time as a power-law function,

g (t) = - G t^{a} {\hat{e}}_{z}, t > 0,

(9)

where G is a prefactor with units $[m / s^{a + 2}], {\hat{e}}_{z}$ is a unit vector in the vertical direction, and a is a dimensionless acceleration exponent. It is worth emphasizing that for time-varying acceleration as a power-law function, the timescale of the early time dynamics is distinguished by the acceleration exponent.^15,19,42 If $a > - 2$ ⁠, the acceleration $g (t)$ determines the timescale and early time dynamics, and this type of instability is termed RT-type instability. If $a < - 2$ ⁠, the initial condition, such as initial growth rate, controls the timescale and early time dynamics, and this type of instability is termed RM-type instability. The RT-to-RM transition occurs at the acceleration exponent a = −2. Here, we focus only on RT-type instabilities, $a > - 2$ ⁠.

C. Group theory

Despite the complexity of RTI and RT mixing, observations demonstrate common features in their large-scale coherent motions.¹ The solutions that satisfy the conservation laws in Eq. (1) and the boundary conditions in Eq. (8) are expected to address the observed universality and symmetric characteristics of RTI and RT mixing. For example, periodicity of the large-scale coherent motions of RTI permits symmetry in a solution describing RT dynamics. This symmetry implies invariance under certain transformations. To serve that purpose, group theory is applied.

For scale-dependent RT dynamics, the group theory approach exploits space groups. Due to the spatial periodicity of RTI, RT dynamics are invariant with respect to the Fedorov and/or Schoenflies group, G. The generators of this Fedorov group G are translations in the horizontal plane, rotations, and reflections.¹¹ The space groups are classified with respect to their dimensionality: there are seven one-dimensional and 17 two-dimensional crystallographic groups.⁴³ However, not all space groups are considered here because structural stability and isotropy significantly limit the types of space groups. The space groups that satisfy the following requirements are selected. First, the space groups must have anisotropy in the acceleration direction. Second, the space group must have inversion in the plane normal to the periodic plane. These requirements ensure that the coherent structures are structurally stable and thus are observable and repeatable. The two requirements physically imply that large-scale fluctuations do not alter the space group G and the corresponding dominant wave vector. Consideration of structural stability and isotropy leaves, for 3D flow, p4 mm, p6 mm, p2 mm, cmm, and p2 among the crystallographic groups following international classification and Fedorov's notation for space groups.⁴³

D. Fourier series expansions

Due to the aforementioned requirements, the irreducible representation of the symmetry group p4 mm allows the application of the Fourier series in combination with the generators of G to solve the governing equations. We consider the square lattice whose spatial period $a_{1 (2)}$ has equal lengths $| a_{1 (2)} | = λ$ and are orthogonal to each other. The wavevectors $k_{1 (2)}$ of the reciprocal lattice are defined as $k_{1 (2)} a_{2 (1)} = 0$ and $k_{1 (2)} a_{1 (2)} = 2 π$ ⁠, thus resulting in $k_{1 (2)} = 2 π / λ$ ⁠.

The spatial period λ serves as a characteristic length scale in the early time and late-time nonlinear dynamics of RTI. Then, λ is the scale that describes large-scale coherent motion. It is also known that a secondary instability, due to the interfacial shear, develops the vortical structures. However, the length scale of the vortical structures is assumed to be small compared to the length scale of the large-scale coherent structures λ. Therefore, a fluid potential is employed to describe the large-scale coherent motion of the fluids as

v_{h, l} = \nabla Φ_{h, l},

(10)

where $Φ_{h, l}$ is the fluid potential of the heavy and light fluid, respectively. With the continuity equation, $Φ_{h, l}$ follows the Laplace equation,

Δ Φ_{h, l} = 0.

(11)

The fluid potentials for the heavy and light fluids for the space group p4 mm may be written as

(12)

Φ_{h} (x, y, z, t) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} Φ_{m n} (t) (z \sinh (α_{m n} k Z) + \frac{cos (mkx) cos (nky)}{α_{m n} k} \cosh [α_{m n} k (z - Z)]) + f_{h} (t),

(12a)

Φ_{l} (x, y, z, t) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} {\tilde{Φ}}_{m n} (t) (- z \sinh (α_{m n} k Z) + \frac{cos (mkx) cos (nky)}{α_{m n} k} \cosh [α_{m n} k (z + Z)]) + f_{l} (t) .

(12b)

Here, α_mn is defined as

α_{m n}^{2} = m^{2} + n^{2}

with integers

m, n = 1, 2, 3 \dots

⁠. Here,

Φ_{m n}

and

{\tilde{Φ}}_{m n}

are the matrices of Fourier amplitudes for the heavier and lighter fluids, respectively, where they satisfy

Φ_{m n} = Φ_{n m}

and

Φ_{00} = 0

for the heavy fluid and similar requirements for the light fluid.

f_{h} (t)

and

f_{l} (t)

are time-dependent functions.

Recall that the space group p4 mm is applied to describe 3D periodicity of the flow. The transformations in the x and y directions are $x \to - x$ and $y \to - y$ as well as $(x, y) \to (- x, - y)$ ⁠. Moreover, not only reflections but also rotations by angles $\pm π / 2, \pm π$ must be applied to the Fourier series expansion. In the vicinity of the tip of the bubble or spike, the spatial expansions of the interface $z^{*} (x, y, t)$ with the aforementioned symmetry constraints only permit

z^{*} (x, y, t) = \sum_{i = 0}^{\infty} \sum_{j = 0}^{\infty} ζ_{i, j} x^{2 i} y^{2 j} .

(13)

Here, we note that $ζ_{i, j}$ is that coefficient of the expansion satisfying $ζ_{i, j} = ζ_{j, i}$ ⁠.

III. RESULTS

In this section, we identify the dependence of the RTI growth rate on the acceleration parameter. For the late-time nonlinear dynamics, we find a continuous family of asymptotic solutions and identify a special family of solutions for bubbles and spikes. Also, the effect of the finite-sized domain kZ on the solutions is examined.

A. Dynamical system

For a given space group, the irreducible representation and project operators with infinite Fourier series enable the velocity potentials $Φ_{h}$ and $Φ_{l}$ to be expanded as forms of Eq. (12). In addition, the interface in Eq. (13) can be expanded at the vicinity of the tips of bubbles and spikes in a given unit cell. Further substitution of the fluid potentials in Eq. (12) and the interface in Eq. (13) into the governing equation in Eq. (8) leads to a dynamical system in terms of moments, M or N, and surface variables ζ.

The moment expressions for M and N are introduced as

(14)

M_{(a, b, c)} = \sum_{a = 1}^{\infty} \sum_{b = 1}^{\infty} \sum_{c = 1}^{\infty} Φ_{m n} (t) {(m k)}^{a} {(n k)}^{b} {(α_{m n} k)}^{c} \sinh (m α_{m n} Z),

(14a)

{\tilde{M}}_{(a, b, c)} = \sum_{a = 1}^{\infty} \sum_{b = 1}^{\infty} \sum_{c = 1}^{\infty} {\tilde{Φ}}_{m n} (t) {(m k)}^{a} {(n k)}^{b} {(α_{m n} k)}^{c} \sinh (m α_{m n} Z),

(14b)

N_{(a, b, c)} = \sum_{a = 1}^{\infty} \sum_{b = 1}^{\infty} \sum_{c = 1}^{\infty} Φ_{m n} (t) {(m k)}^{a} {(n k)}^{b} {(α_{m n} k)}^{c} \cosh (m α_{m n} Z),

(14c)

{\tilde{N}}_{(a, b, c)} = \sum_{a = 1}^{\infty} \sum_{b = 1}^{\infty} \sum_{c = 1}^{\infty} {\tilde{Φ}}_{m n} (t) {(m k)}^{a} {(n k)}^{b} {(α_{m n} k)}^{c} cos h (m α_{m n} Z),

(14d)

where m and n are the integers

m, n = 1, 2, 3 \dots

⁠. Note that the moments with tilde,

\tilde{M}

and

\tilde{N}

⁠, represent the moments for the light fluid; otherwise, the moments are for the heavy fluid. We note that the Fourier matrix

Φ_{m n}

for the heavy fluids (and, similarly, for the light fluid) appears due to the independent translational directions. Yet, wavevectors

| k_{1} |

and

| k_{2} |

have the same magnitude.

To proceed further, the pressure term in the first two nontrivial conditions in Eq. (8) is replaced by the pressure term in the Bernoulli equation in the frame accelerating at $d v / d t$ ⁠, which is

\frac{\partial Φ_{h, l}}{\partial t} + P + (g + \frac{\partial v}{\partial t}) z + \frac{1}{2} | v |^{2} = C (t),

(15)

where C(t) is a function of time. In addition, the boundary conditions in the accelerating frame are modified as

v_{z} (z = Z) = - v (t) \cdot {\hat{e}}_{z}, v_{z} (z = - Z) = - v (t) \cdot {\hat{e}}_{z},

(16)

where ${\hat{e}}_{z}$ is a unit vector in z direction.

Then, the substitution of moment expressions in Eq. (14) into Eq. (8) with Eqs. (15) and (16) leads to the dynamical system with the equations

(17)

\dot{ζ_{1}} - (2 N_{(2, 0, - 1)} + N_{(0, 0, 1)}) ζ_{1} - \frac{M_{(2, 0, 0)}}{2} = 0,

(17a)

\dot{ζ_{1}} - (2 {\tilde{N}}_{(2, 0, - 1)} + {\tilde{N}}_{(0, 0, 1)}) ζ_{1} + \frac{{\tilde{M}}_{(2, 0, 0)}}{2} = 0,

(17b)

\begin{matrix} (1 + A) (\frac{{\dot{N}}_{(2, 0, 0)}}{2} + ζ_{1} {\dot{M}}_{0} - \frac{N_{(2, 0, 0)}^{2}}{4} - g (t) ζ_{1}) \\ = (1 - A) (\frac{{\dot{\tilde{N}}}_{(2, 0, 0)}}{2} - ζ_{1} {\dot{\tilde{M}}}_{0} - \frac{{\tilde{N}}_{(2, 0, 0)}^{2}}{4} - g (t) ζ_{1}), \end{matrix}

(17c)

M_{0} + {\tilde{M}}_{0} = 0.

(17d)

Here,

ζ_{1} = ζ_{0, 1} = ζ_{1, 0}

is the curvature and

A = (ρ_{h} - ρ_{l}) / (ρ_{h} + ρ_{l})

is the Atwood number. Dot notation denotes the time derivative. Also, note that

M_{0} = - v (t) \cdot {\hat{e}}_{z} = - v_{3}

and

{\tilde{M}}_{0} = v (t) \cdot {\hat{e}}_{z} = v_{3}

⁠. It is worth to emphasize that the dynamical system is composed of an infinite power series in terms of x with time-dependent coefficients

ζ_{m} (t)

⁠. Thus, the dynamical system can be spatially expanded to the desired order of accuracy.

Also, we note that various choices of parameter sets (a, b, c) are possible for the moments of the space group p4 mm. For conciseness, we replace the moments in the dynamical system in Eq. (17) with independent moments. From the definition of the moment expression for p4 mm in Eq. (14), we easily notice that for any kind of moments, such as $M, \tilde{M}, N, \tilde{N}$ ⁠, the following is satisfied:

M_{(a, b, c)} = M_{(b, a, c)} .

(18)

Applying Eq. (18) to each of the moments M, $\tilde{M}$ ⁠, N, and $\tilde{N}$ in Eq. (17) results in sets of two independent moments. For example, $M_{1} = M_{(0, 0, 1)}$ and $M_{2} = M_{(0, 0, 2)}$ for the moment M, where

M_{(0, 0, 1)} = 2 M_{(2, 0, - 1)} = M_{(2, 0, - 1)} + M_{(0, 2, - 1)},

(19)

M_{(0, 0, 2)} = 2 M_{(2, 0, 0)} = M_{(2, 0, 0)} + M_{(0, 2, 0)} .

(20)

Exploiting the independent set of moments, the dynamical system of p4 mm can be written in simplified form as

(21)

\dot{ζ_{1}} - 2 N_{1} ζ_{1} - \frac{M_{2}}{4} = 0,

(21a)

\dot{ζ_{1}} - 2 {\tilde{N}}_{1} ζ_{1} + \frac{{\tilde{M}}_{2}}{4} = 0,

(21b)

\begin{matrix} (1 + A) (\frac{{\dot{N}}_{1}}{4} + ζ_{1} {\dot{M}}_{0} - \frac{N_{1}^{2}}{8} - g (t) ζ_{1}) \\ = (1 - A) (\frac{{\dot{\tilde{N}}}_{1}}{4} - ζ_{1} {\dot{\tilde{M}}}_{0} - \frac{{\tilde{N}}_{1}^{2}}{8} - g (t) ζ_{1}), \end{matrix}

(21c)

M_{0} + {\tilde{M}}_{0} = 0.

(21d)

The dynamical system in Eq. (21) is characterized by one length scale and two natural time scales. The length scale is determined by the wavevector $k = 2 π / λ$ ⁠, where λ is the spatial period set by the initial conditions. The two time scales are the time scale related to accelerations, τ_G, and the time scale related to initial growth rate, τ₀. The acceleration-related time scale, τ_G, is set by the acceleration parameters a and G in Eq. (9) as $τ_{G} = {(k G)}^{- 1 / (a + 2)}$ ⁠, whereas the initial-growth-rate-related time scale τ₀ is determined as $τ_{0} = {(k v_{0})}^{- 1}$ where v₀ is the initial growth-rate. For $a > - 2$ ⁠, the time scales relate as $τ_{G} ≪ τ_{0}$ ⁠, and the fastest process with the smallest time scale is defined by the acceleration, leading to RT-type acceleration-driven dynamics. For $a < - 2$ ⁠, the time scales relate as $τ_{G} ≫ τ_{0}$ ⁠, and the fastest process with the smallest time scale is defined by the constant initial growth-rate, leading to RM-type initial-growth-rate-driven dynamics.^15,19 RTI considered in this study is acceleration-driven and of RT-type such that $τ_{G} ≪ τ_{0}$ and $a > - 2$ ⁠, respectively. Therefore, the time scale of RTI here is $τ = τ_{G}$ ⁠. Moreover, we assume that the initial time t₀ follows $t_{0} ≫ τ_{G}, τ_{0}$ meaning the RT dynamics is independent of the choice of time origin t₀.

B. The shear function

In the nonlinear regime of RTI, the flow features consist of not only large-scale coherent structures but also secondary instabilities caused by the small-scale vortices. It is understood that the source of these vortices is shear that develops across the interface. Therefore, shear serves as an important quantity to be parameterized and is of interest to RTI dynamics. The shear function Γ is defined as and to first order is

Γ = \lim_{x \to 0} \frac{1}{x} (v_{h} (x, y) \cdot n_{t} - v_{l} (x, y) \cdot n_{t}) = {\tilde{N}}_{1} - N_{1} .

(22)

Note that the working fluids are inviscid in the current study. Thus, no viscous layer exists at the vicinity of the interface. Even though this inviscid framework limits our study to large-scale dynamics, the introduction of the shear function allows us to describe the qualitative properties of vortical structures and their rotation direction. For example, positive Γ indicates the rotation from the heavy fluid to the light fluid.

Hereafter, we introduce symbolic letters for clarity in describing quantities and solutions of bubbles and spikes. Any quantity χ evaluated in a finite-sized domain is denoted with an underline as $\underline{χ}$ ⁠. Quantities without underlines indicate those evaluated in an infinite-sized domain (e.g., χ). The finite quantity $\underline{χ}$ obviously approaches the infinite quantity χ as the domain size increases $k Z \to \infty$ ⁠. Therefore, solutions in the finite-sized domain can be extended to those of infinite domains without losing generality. Also, properties related to spikes are marked with hats as $\hat{χ}$ ⁠, whereas those of bubbles are simply represented as χ. Moreover, the curvature ζ₁ to the first approximation $N = i + j = 1$ in Eq. (13) is replaced by ζ for brevity.

C. Early time dynamics

The early time dynamics describe the linear regime where $(t - t_{0}) ≪ τ$ ⁠. In the linear regime, the initial conditions for the curvature $ζ (t_{0})$ and the velocity $v (t_{0})$ are assumed to be small, leading to $0 < | ζ (t_{0}) k | ≪ 1$ and $0 < | v (t_{0}) {(k / g)}^{1 / 2} | ≪ 1$ ⁠. Thus, the moments in Eq. (14) are approximated to first order.

Applying the aforementioned assumptions to the dynamical system in Eq. (21) leads to the solutions in the linear regime as

\begin{matrix} \underline{ξ} = C_{1} \sqrt{\frac{t}{τ}} I_{\frac{1}{2 s}} (\sqrt{A \tanh (k Z)} \frac{{(t / τ)}^{s}}{s}) \\ + C_{2} \sqrt{\frac{t}{τ}} I_{- \frac{1}{2 s}} (\sqrt{A \tanh (k Z)} \frac{{(t / τ)}^{s}}{s}), \end{matrix}

(23)

\underline{v} = - \frac{4}{k} \frac{d \underline{ξ}}{d t},

(24)

where $s = (a + 2) / 2$ and I_p is the modified Bessel function of the pth order. Also, ξ is a dimensionless curvature defined as $ξ = ζ / k$ ⁠. The curvature ζ has a negative value for bubbles as it is concave down (⁠ $ζ < 0$ ⁠) in the considered configuration (see Fig. 1) and has a positive value for spikes. We note that bubbles move upward (v >0) while spikes move downward (⁠ $\hat{v} < 0$ ⁠). The magnitudes of solutions for bubbles and spikes in the linear regime are identical.

The dimensionless curvature ξ and velocity $v / \sqrt{g / k}$ for bubbles and spikes with different acceleration parameters a and Atwood numbers A are shown in Fig. 2 for two finite-sized domains. As the acceleration exponent a increases from a = −2 to a = 2, bubbles and spikes are more curved and move faster. Both bubbles and spikes have larger curvature and faster velocity for larger domain size (kZ = 5) compared to the smaller domain size (kZ = 0.5). Fixing the acceleration parameter as a = 2, velocity and curvature increase as the Atwood number increases.

FIG. 2.

View large Download slide

The linear regime solutions for the curvature ζ and velocity v for bubbles and spikes for the subspace group p4 mm. The linear regime solutions are shown with different acceleration exponents (a) and (b). The acceleration exponents are a = −2 (dash-dotted lines), a = 0 (dashed line), and a = 2 (solid line). The solutions are shown with different Atwood numbers (c) and (d). The Atwood numbers are A = 0.15 (dash-dotted lines), A = 0.5 (dashed lines), and A = 0.95 (solid line). The effect of different finite-sized domains is also illustrated in all panels, kZ = 5 (black lines) and kZ = 0.5 (gray lines).

D. Late-time dynamics

The late-time dynamics with a power-law acceleration in time in Eq. (9) show complicated features compared to the linear regime. For the late-time dynamics, we seek a one-parameter family of asymptotic solutions in the form of a power-law function hinted by the same form of acceleration in Eq. (9).

1. One-parameter family of solutions

Garabedian⁴⁴ first suggested the existence of one-parameter family of asymptotic solutions for two-dimensional periodic RTI for a one-fluid system. A rigorous analytical analysis was given to the stationary solutions for three-dimensional periodic RTI with constant acceleration for a one-fuid system.⁴⁵ The importance of the flow symmetry and the validity of the parameterized family of solutions was investigated in terms of uniqueness, dimensionality, and universality.⁴⁶

Previous studies demonstrated that the symmetry of the flow is a crucial characteristic in order to pose a constraint and permit a unique symmetric expression [e.g., see Eq. (12)] among the general solutions satisfying the Laplace equation in Eq. (11). In other words, the flow can be described by a one-parameter family of asymptotic solutions under a symmetric condition p4 mm; both the fluid potentials $Φ$ in Eq. (12) and the free surface $z^{*} (x, y, t)$ in Eq. (13) are invariant with respect to the spatial group.⁴⁶

To illustrate, the number of equations N_e is related to the number of unknown variables N_t as $N_{t} = N_{e} + 2$ for the order N approximation of the Fourier series expansions in Eq. (12). Manipulating $N_{e} - 1$ equations yields an expression $f (v, R_{1}, R_{2}) = 0$ ⁠, where R₁ and R₂ are the curvature of the tip in each of the independent horizontal directions, from which we find the velocity v and the harmonics $Φ_{m n}$ as functions of the parameters R₁ and R₂, i.e., a two-parameter solution set in the Nth approximation.⁴⁶ However, if the flow is symmetric, the radii of the curvature R₁ and R₂ are dependent on the flow symmetry. Thus, this dependence removes one undetermined variable leading to $N_{t} = N_{e} + 1$ (cf., for zeroth-order solutions N_t = N_e and the curvature R₁ and R₂ are fixed).

Following this analogy, space groups containing a fourth-order axis (p4 mm) are shown to have a one-parameter family of asymptotic solutions for 3D periodic RTI. On the other hand, for the rectangular symmetry group p2 mm, the flow permits a two-parameter family of solutions due to the two independent translations, $k_{1}$ and $k_{2}$ ⁠, in the xy plane.

Here, we retain four unknowns for $Φ_{m n}$ in order to obtain a one-parameter family of solutions. The number of equations N_e is N_e = 4 as provided in Eq. (21), and also $N_{t} = N_{e} + 1 = 4 + 1 = 5$ since the curvatures in x and y directions are dependent due to the symmetry condition in Eq. (13) with $i + j = 1$ ⁠. Then, we derive a solution, for example, the velocity at the tip v, as a function of the other unknown variable, such as the curvature ζ.

Recall that the Fourier matrix $Φ_{m n}$ has two-dimensional elements depending on the row and column index m and n. As the matrix indices increases, more unknowns in the dynamical system are introduced. Hence, different sets of Fourier amplitudes can be considered for $Φ_{m n}$ $m + n < 3$ ⁠, for instance, $(Φ_{10}, Φ_{20}, {\tilde{Φ}}_{10}, {\tilde{Φ}}_{20})$ and $(Φ_{10}, Φ_{11}, {\tilde{Φ}}_{10}, {\tilde{Φ}}_{11})$ ⁠. Therefore, the choice of the elements in the Fourier matrix requires all possible choices to be compared and to be justified. Earlier works address this justification for the ideal family of solutions for spatially periodic solutions of RTI.^45,46 Suppose such a family of solutions exists and is unique. Then, there must exist a functional limit with respect to the free parameter (e.g., the curvature ζ) of the Fourier amplitudes $Φ_{m n}$ as the order N of the approximation increases

Φ (ζ) = \lim_{N \to \infty} Φ_{N} (ζ) .

(25)

The smoothness of the family of solutions was also demonstrated.^45,46 This smoothness implies the convergence of the Fourier harmonics as the order increases $N \to \infty$ within the functional limits set by the free parameter

| Φ_{m n} |_{m + n = l + 1} ≪ | Φ_{m n} |_{m + n = l} .

(26)

Thus, the absolute value of the amplitude $Φ_{m n}$ decays as l increases. Moreover, the solutions corresponding to each approximation should be similar. In other words, no significant difference is observed between $Φ_{11}$ and $Φ_{20}$ ⁠.⁴⁵ As a result, we justify our choice of Fourier moments for p4 mm as $Φ_{10}, Φ_{20}, {\tilde{Φ}}_{10}$ ⁠, and ${\tilde{Φ}}_{20}$ in the present study.

2. Asymptotic solutions

In order to derive the asymptotic solution, the following scalings for the terms in Eq. (21) are assumed:

\begin{matrix} ζ \sim k {(\frac{t}{τ})}^{β}, Φ_{m}, {\tilde{Φ}}_{m} \sim \frac{1}{k τ} {(\frac{t}{τ})}^{α}, \\ M_{m}, N_{m}, {\tilde{M}}_{m}, {\tilde{N}}_{m} \sim \frac{k^{n}}{k τ} {(\frac{t}{τ})}^{α}, v \sim \frac{1}{k τ} {(\frac{t}{τ})}^{α}, \end{matrix}

(27)

where α and β are determined by balancing the orders of magnitude of the terms in Eq. (21). Substituting Eq. (27) to Eq. (21) and balancing the exponents of the terms lead to two possible pairs of exponents $(α, β)$ depending on the acceleration exponent a.

The power-law exponents of the first pair for RT-type instability (⁠ $a > - 2$ ⁠) are $α = a / 2$ ⁠, β = 0. The power-law exponents of the second pair for RM-type instability (⁠ $a < - 2$ ⁠) are $α = - 1, β = 0$ ⁠. We note that for the acceleration exponent a = −2, the transition from the first case to the second case occurs. As stated earlier, the present study considers RT-type instabilities. Thus, the acceleration exponent satisfies $a > - 2$ ⁠. Substituting Eq. (27) and $α = a / 2$ ⁠, β = 0 into Eq. (21) and retaining only the dominant terms lead to a dynamical system in the nonlinear regime,

(28)

- 2 N_{1} ζ - \frac{M_{2}}{4} = 0,

(28a)

- 2 {\tilde{N}}_{1} ζ + \frac{{\tilde{M}}_{2}}{4} = 0,

(28b)

\begin{matrix} (1 + A) (\frac{1}{8} N_{1}^{2} - g (t) ζ) = (1 - A) (- \frac{1}{8} {\tilde{N}}_{1}^{2} - g (t) ζ), \end{matrix}

(28c)

M_{0} + {\tilde{M}}_{0} = 0.

(28d)

Note again that the acceleration g(t) follows

g (t) = \frac{1}{k τ^{2}} {(\frac{t}{τ})}^{a}

⁠.

Then, the moments in Eq. (14) are truncated to the second order with $Φ_{10}, Φ_{20}, {\tilde{Φ}}_{10}$ ⁠, and ${\tilde{Φ}}_{20}$ as previously discussed. As a result, we arrive at the solutions for Eq. (28), which leads to the parameterized family of solutions,

\underline{v} = \frac{\sqrt{2}}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} \frac{{(A \underline{ξ})}^{\frac{1}{2}} \tanh (k Z) (3 + 8 \underline{ξ} \tanh (k Z)) (3 - 8 \underline{ξ} \tanh (k Z))}{{[{(- 3 + \tanh {(k Z)}^{2})}^{2} {48 \tanh (k Z) \underline{ξ} - A (9 + 64 \tanh {(k Z)}^{2} {\underline{ξ}}^{2})}]}^{\frac{1}{2}}},

(29)

\underline{Γ} = \frac{12 \sqrt{2}}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} \frac{{(A \underline{ξ})}^{\frac{1}{2}}}{{48 \tanh (k Z) \underline{ξ} - A (9 + 64 \tanh {(k Z)}^{2} {\underline{ξ}}^{2})}^{\frac{1}{2}}} .

(30)

Recall that bubbles move upward (v >0) and concave down (⁠ $ζ < 0$ ⁠), whereas spikes move downward (v <0) and concave up (⁠ $ζ > 0$ ⁠).

E. Universal nonlinear scaling

The effect of the finite domain size kZ appears in the solutions throughout the linear regime, Eqs. (23) and (24), and the nonlinear regime, Eqs. (29) and (30), in the form of a hyperbolic function, such as $\tanh (k Z)$ ⁠. However, how the finite domain sizes are related to the solutions is distinct for the two different regimes. This distinct effect of the finite domain size kZ enables a universal scaling in the nonlinear regimes for the parameters of interest.

The domain extent kZ can be related to the original length scale of RTI, which is the horizontal length scale λ in the nonlinear regime. In the linear regime, there exists a scale separation between the heights of the large coherent structure h and the finite domain kZ since the initial velocity and amplitude are assumed to be small. Thus, they are independent of each other. On the other hand, in the nonlinear regime, the vertical length scale h of bubbles and spikes increases to a value comparable to horizontal length scale of the system, $h \sim O (1 / k)$ ⁠. This similar order of magnitude in length scale implies the possibility of dependence between the characteristic scales, such as k and g(t), and the domain size kZ.

Following this analogy, the universal scaling bridging the solution in the finite domain to the infinite domain is established for each characteristic scale in the RTI system. In order to demonstrate this scaling, we start from the definition of the Atwood bubble, which is the fastest solution of the parameterized family of solution in Eqs. (29) and (30). The Atwood bubble satisfies

{\frac{\partial \underline{v}}{\partial \underline{ζ}} |}_{\underline{ζ} = {\underline{ζ}}_{A}} = 0, {\frac{\partial^{2} \underline{v}}{\partial {\underline{ζ}}^{2}} |}_{\underline{ζ} = {\underline{ζ}}_{A}} < 0,

(31)

where ζ_A is the curvature ζ for the Atwood bubble and $ζ_{A} \in [0, ζ_{c r}]$ ⁠, where ζ_cr is the maximum possible curvature for the nonlinear bubble.

In order to compute ζ_A, differentiation of Eq. (29) with respect to the curvature ζ leads to the following characteristic equation:

\begin{matrix} 4096 A {({\underline{ξ}}_{A} \tanh (k Z))}^{4} + 4096 {({\underline{ξ}}_{A} \tanh (k Z))}^{3} \\ + 1152 A {({\underline{ξ}}_{A} \tanh (k Z))}^{2} - 27 A = 0, \end{matrix}

(32)

where ${\underline{ξ}}_{A} = {\underline{ζ}}_{A} / k$ ⁠. Equation (32) implies the form of transformation for the length scale in an infinite-sized domain k and the length scale in finite-sized domain $\underline{k}$ as

\underline{k} = S_{k} k, S_{k} = \tanh (k Z) .

(33)

In addition, we define the following fundamental scales:

(34)

S_{v} = \sqrt{\frac{S_{g}}{S_{k}}},

(34a)

S_{τ} = \frac{1}{\sqrt{S_{k} S_{g}}},

(34b)

S_{ω} = \frac{1}{S_{τ}},

(34c)

S_{λ} = \frac{2 π}{S_{k}},

(34d)

where

S_{χ}

denotes the scaling factor linking

\underline{χ}

and χ. Then,

S_{g}

is derived by substituting Eq. into Eq. and employing Eq. as

S_{g} = \frac{{(- 3 + \tanh^{2} (k Z))}^{2}}{4} .

(35)

Thus, Eqs. (33)–(35) define a complete set of fundamental scales in the RTI system. Universal scales in Eqs. (33)–(35) are illustrated in Fig. 3. All scales for the length and acceleration asymptotically reach unity as the finite domain size increases. As the domain size kZ is reduced, the scales for the curvature $S_{k}$ and the velocity $S_{v}$ decreases. On the other hand, the scales for the gravity $S_{g}$ and characteristic time $S_{τ}$ increases. The scale for the shear, which is equivalent with $S_{ω}$ ⁠, is not monotonic. However, we note that our analysis in the nonlinear regime is valid for $k Z ⪆ 1$ ⁠. Thus, the scale for shear increases as the domain size is decreased to kZ = 1. These trends of scale variations indicate the modulations of the corresponding quantities in the same tendency. For example, the curvature ζ and velocity v decrease, whereas the gravity g, shear Γ, and time scale τ increase. Beyond kZ > 3, less than 1% of the domain size effect is observed.

FIG. 3.

View large Download slide

Universal scales as functions of the finite domain size kZ: $S_{k}$ (black solid line), $S_{g}$ (black dashed line), $S_{v}$ (black dash-dotted line), $S_{ω}$ (gray solid line), and $S_{τ}$ (gray dashed line).

The universal nonlinear scales not only provide the physical insights on the effect of the finite domain size on different characteristic parameters but also enable a significant efficiency when describing solutions in the finite domain. For example, the velocity v in Eq. (29) and the shear Γ in Eq. (30) are rewritten as

(36)

v = \frac{\sqrt{2}}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(A ξ)}^{\frac{1}{2}} (3 - 8 ξ) (3 + 8 ξ) D,

(36a)

Γ = \frac{12 \sqrt{2}}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(A ξ)}^{\frac{1}{2}} D,

(36b)

where

D

is defined as

D = {48 ξ - A (9 + 64 ξ^{2})}^{- 1 / 2}

⁠.

Moreover, the Fourier amplitudes are

(37)

Φ_{1} = \frac{- 4 \sqrt{2}}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(A ξ)}^{\frac{1}{2}} (3 - 8 ξ) (1 + 4 ξ) D,

(37a)

Φ_{2} = \frac{\sqrt{2}}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(A ξ)}^{\frac{1}{2}} (3 - 8 ξ) (1 + 8 ξ) D,

(37b)

{\tilde{Φ}}_{1} = \frac{4 \sqrt{2}}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(A ξ)}^{\frac{1}{2}} (3 + 8 ξ) (1 - 4 ξ) D,

(37c)

{\tilde{Φ}}_{2} = \frac{- \sqrt{2}}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(A ξ)}^{\frac{1}{2}} (3 + 8 ξ) (1 - 8 ξ) D .

(37d)

The effects of the Atwood number on the family of the solutions in the nonlinear regimes are illustrated in Fig. 4. The curvature ζ is normalized by the critical curvature $ζ_{c r} = ζ_{max}$ ⁠, which is the maximum curvature of bubbles and spikes in nonlinear regime. As the Atwood number varies, fluids with a higher contrasting density (A = 0.95) show faster velocity and larger shear than fluids with a similar density (A = 0.15). In Figs. 4(c) and 4(d), the velocity v and shear Γ for spikes diverge at a certain curvature. This blow-up of the solutions indicates that the solutions for spikes of the form in Eqs. (29) and (30) are applicable only for a finite range of curvatures. For higher Atwood numbers, this range of curvature narrows.

FIG. 4.

View large Download slide

The nonlinear-regime family of asymptotic solutions for bubbles (a) and (b) and spikes (c) and (d) for the p4 mm subspace group. The velocity $v / \sqrt{g / k}$ ⁠, shear $Γ / \sqrt{g k}$ ⁠, and curvature ξ are shown for three different Atwood numbers: A = 0.15 (dash-dotted lines), A = 0.5 (dashed lines), and A = 0.95 (solid lines).

It is also informative to see how the maximum value of the solutions is achieved depending on the Atwood number. For this purpose, the properties of bubbles are re-scaled by their maximum values. For example, the velocity v is normalized by the Atwood velocity, which has the fastest solution, $v_{A} = v_{max}$ ⁠. Similarly, the shear Γ is normalized by the critical shear $Γ_{c r} = Γ_{max}$ ⁠.

We take a closer look at the solutions for bubbles and spikes separately. First, the solutions for bubbles normalized by their maximum values are shown in Fig. 5. Figure 5(a) shows the maximum velocity and the corresponding curvature. Note that in the theoretical family of solutions, bubbles are less curved when the density contrast of the fluids is small. On the other hand, more shear is exerted on bubbles of lower Atwood number. This larger shear is in agreement with slower velocity since the shear acts as a drag on the bubble [Figs. 5(b) and 5(c)].

FIG. 5.

View large Download slide

One-parameter family of asymptotic solutions for bubbles in the nonlinear regime for p4 mm. The solutions are normalized by their maximum values. The velocity $v / v_{max}$ and shear $Γ / Γ_{max}$ as functions of the curvature $ζ / ζ_{c r}$ are given in (a) and (b), respectively; (c) the one-to-one mapping of the re-scaled velocity and shear. (d) The first two Fourier amplitudes are illustrated as functions of the Atwood number A. The effect of the Atwood number A = 0.15 (dash-dotted line), A = 0.5 (dashed line), and A = 0.95 (solid line) is shown. The Fourier amplitudes for the heavy (black lines) and light (gray lines) fluids are plotted for the first and second harmonics, which are represented with a solid line and dashed line, respectively.

Moreover, the convergence behavior of the Fourier amplitude for the solutions of Atwood bubble is presented in Fig. 5(d). Since the first harmonic is large in magnitude compared to the second harmonic, the nonlinear solutions for bubbles are convergent for the entire range of Atwood numbers $0 < A \leq 1$ ⁠. Also, no singularity is found in the solutions as the bubble solutions are regular and stable. This regularity is in agreement with the experiments.⁴¹

F. Singularity of the solutions for RT spikes

Next, we investigate the nonlinear solutions of spikes. Unlike the family solutions of bubbles, the solutions of spikes exhibit a singular behavior, which is also in agreement with the experiments.^1,41 This singularity originates from the pole of the denominator of the velocity in Eq. (36), which is

- 48 \hat{ξ} + A (9 + 64 {\hat{ξ}}^{2}) = 0.

(38)

The root of Eq. (36) satisfying $\hat{ξ} > 0$ is ${\hat{ξ}}_{c r} = 3 (1 - \sqrt{1 - A^{2}}) / (8 A)$ ⁠. As the curvature reaches this singular point $\hat{ζ} \to {\hat{ζ}}_{r}$ ⁠, the velocity $\hat{v}$ approaches the negative infinity (⁠ $\hat{v} \to - \infty$ ⁠). Thus, the velocity at $\hat{ζ} = {\hat{ζ}}_{r}$ corresponds to the maximum velocity ${\hat{v}}_{max}$ for spikes (⁠ $| {\hat{v}}_{max} | = \infty$ ⁠).

Due to the existence of the singularity, it is worthwhile to clarify the domain of the curvature $\hat{ζ}$ ⁠. Recall the definition of the critical curvature $\hat{ζ} = {\hat{ζ}}_{c r}$ ⁠, at which the velocity of the spike is zero (⁠ ${\hat{v}}_{c r} = 0$ ⁠). This definition naturally implies $0 < {\hat{ζ}}_{r} < {\hat{ζ}}_{c r}$ for spikes. Therefore, the domain for the spike curvature is $\hat{ζ} \in [{\hat{ζ}}_{r}, {\hat{ζ}}_{c r}]$ and the corresponding interval for velocity is $\hat{v} \in [{\hat{v}}_{max}, 0]$ ⁠.

As a conservative measure for singularity, the solutions of spikes at $\hat{ζ} = {\hat{ζ}}_{r}$ are regularized by

{\bar{\hat{v}}}_{r} = {\hat{v}}_{r} {(\frac{\hat{ζ} - {\hat{ζ}}_{r}}{k})}^{1 / 2}, {\bar{\hat{Γ}}}_{r} = {\hat{Γ}}_{r} {(\frac{\hat{ζ} - {\hat{ζ}}_{r}}{k})}^{1 / 2} .

(39)

The solutions of spikes in the nonlinear regime regularized by quantities in Eq. (39) are shown in Fig. 11.

In contrast to the solutions for bubbles, the solutions for spikes are evaluated in the curvature domain $\hat{ζ} \in [{\hat{ζ}}_{r}, {\hat{ζ}}_{c r}]$ ⁠, as previously demonstrated. As the Atwood number increases, the corresponding curvature for the fastest solution slightly increases. Thus, spikes are more curved for the fluid pairs with higher density contrast. Due to the singularity, the curvature for the maximum velocity is limited and determined by the edge of the range [Figs. 6(a) and 6(b)]. Moreover, the shear monotonically increases as the curvature approaches the critical value. In Fig. 6(c), we observe that the velocity decreases as the shear increases for A = 0.95 while the other Atwood numbers, A = 0.5 and A = 0.15, have local maxima. In other words, the shear always decelerates spikes leading to slower velocities when the density difference is large. However, for the cases when the density contrast of the fluids is not too large, more curved spikes due to large shear can lead to faster velocities. Moreover, the convergence behavior of the Fourier amplitude for Atwood spike is presented in Fig. 6(d). We note that the convergence holds for $0 < A < 45 / 53 \approx 0.85$ since the solutions for the spikes in the nonlinear regime are singular and unstable. Within the convergence regime, the first harmonic of Fourier amplitudes $| Φ_{1} |$ and $| {\tilde{Φ}}_{1} |$ is dominant (e.g., $| Φ_{1} | ≫ | Φ_{2} |$ and $| {\tilde{Φ}}_{1} | ≫ | {\tilde{Φ}}_{2} |$ ⁠) even though the solutions are singular.

FIG. 6.

View large Download slide

One-parameter family of asymptotic solutions for spikes in the nonlinear regime with the subspace group p4 mm. The solutions are normalized by their maximum values. The velocity $v / {\hat{\underline{v}}}_{r}$ and shear $Γ / {\hat{\underline{Γ}}}_{r}$ as functions of the curvature $ζ / {\hat{\underline{ζ}}}_{c r}$ are given in (a) and (b), respectively. (c) The one-to-one mapping of the rescaled velocity and shear; (d) the first two Fourier amplitudes as a function of the Atwood number A. In (a)–(c), the effect of the Atwood number A = 0.15 (dash-dotted line), A = 0.5 (dashed line), and A = 0.95 (solid line) is shown. In (d), the Fourier amplitudes for the heavy (black lines) and light (gray lines) fluids are plotted for the first and second harmonics, which are represented by solid lines and dashed lines, respectively.

Another different feature of solutions for spikes in the nonlinear regime compared to those for bubbles is the large magnitude of the shear. The shear approaches the negative infinity (⁠ $\hat{Γ} \to - \infty$ ⁠) as $\hat{ζ} \to {\hat{ζ}}_{r}$ ⁠. Therefore, the magnitude of the shear is maximized (⁠ $\hat{Γ} = {\hat{Γ}}_{max}$ ⁠) at the singular point $\hat{ζ} = {\hat{ζ}}_{r}$ ⁠. It is known that the interfacial shear induces small-scale vortical structures at the vicinity of spikes and bubbles. Hence, these vortical structures are significantly pronounced at the interface in the vicinity of spikes.^30,32,41

The effect of the finite-sized domain kZ on the fastest solutions of bubbles and spikes are shown for three different Atwood numbers in Fig. 7. As the Atwood number is increased, the magnitude of the curvature ζ, the velocity v, and the shear Γ are increased. Analogous to Fig. 3, domain size has a negligible effect on the solutions beyond kZ = 3. The solutions of spikes are more sensitive to the Atwood number than those of bubbles. The ratios of the magnitude of the curvature, velocity, and shear for kZ = 4 are 1.43, 1.7, and 1.8 for bubbles and 6.6, 9.6 and 13.9 for spikes, respectively.

FIG. 7.

View large Download slide

The effect of the finite-sized domain kZ on the fastest solutions of one-parameter family of solutions. Figures on the left column are the curvature (a), the velocity (c), and the shear (e) for bubbles. Figures on the right column are the curvature (b), the velocity (d), and the shear (f) for spikes. Solutions for the spikes are regularized. Solutions are plotted with Atwood numbers A = 0.15 (dash-dotted line), A = 0.5 (dashed line), and A = 0.95 (solid line).

G. Special solutions: Bubbles

The family of solutions in the nonlinear regime in Eq. (36) describes all the possible values for the velocities v and shear Γ at the tips of bubbles and spikes parameterized by the curvature ξ. Thus, special solutions with specific properties for the diagnostic parameters can be identified from the general form. The specific properties include a particular value for the curvature or velocity. We note that the diagnostic parameters, velocity v, shear Γ, and curvature ξ, are related to each other in the nonlinear regime (e.g., see Fig. 5). Therefore, special solutions are characterized by one of the diagnostic parameters (e.g., ξ). The solutions highlighted here are the critical bubble, the convergence-limit bubble, the Taylor bubble, the Layzer-drag bubble, the Atwood bubble, and the flat bubble. Note that all the quantities presented herein are evaluated in the infinite domain.

1. Critical bubble

The critical bubble is defined as the bubble with the maximum curvature ζ_cr. At $ζ = ζ_{c r}$ ⁠, the bubble velocity is zero (v_cr = 0), and the shear Γ_cr achieves its maximum value Γ_max as follows:

ζ_{c r} = - \frac{3}{8} k, v_{c r} = 0, Γ_{c r} = \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{\frac{6 A}{1 + A}} .

(40)

The solutions in Eq. (40) are shown in Fig. 8.

FIG. 8.

View large Download slide

Special solutions for bubbles among the one-parameter family of asymptotic solutions with p4 mm. (a) Curvature, (b) velocity, and (c) shear are illustrated as a function of Atwood number A. A one-to-one relation between the shear and the velocity is presented in (d). The various lines correspond to the solution, the critical bubble (dash-dotted lines), the convergence-limit bubble (dotted lines), Taylor bubble (dashed lines), Layzer-drag bubble (gray solid lines), and Atwood bubble (black solid lines).

2. Convergence-limit bubble

The convergence-limit bubble is categorized as the bubble whose curvature ζ_cn is at the limit of the convergence of the Fourier amplitudes. More specifically, the amplitudes for the lighter fluid in Eq. (37) satisfy $| {\tilde{Φ}}_{1} | > | {\tilde{Φ}}_{2} |$ for $ζ \in (ζ_{c r}, 0)$ ⁠, whereas the amplitudes for the heavier fluid are $| Φ_{1} | > | Φ_{2} |$ only at $ζ \in (ζ_{c n}, 0)$ ⁠. Figure 9 shows the Fourier amplitudes of the heavier fluid on semilogarithmic axes, along with the convergence region. The curvature where the two heavy-fluid amplitudes are of the same magnitude $| Φ_{1} | = | Φ_{2} |$ is defined as the convergence curvature $ζ_{c n} = - 5 / 24$ ⁠. Therefore, the convergence-limit bubble solution is expressed as

\begin{matrix} ζ_{c n} = - \frac{5}{24} k, v_{c n} = \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} \frac{14}{3} \sqrt{\frac{10 A}{3 (45 + 53 A)}}, \\ Γ_{c n} = \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} 9 \sqrt{\frac{10 A}{3 (45 + 53 A)}} . \end{matrix}

(41)

FIG. 9.

View large Download slide

Fourier amplitudes of the heavier fluid: $| Φ_{1} | / \sqrt{g / k}$ (solid line) and $| Φ_{2} | / \sqrt{g / k}$ (dashed line).

To compare the convergence-limit bubble to the critical bubble, the ratio of the curvatures is $ζ_{c n} / ζ_{c r} = 5 / 9$ ⁠. For any Atwood number A, the convergence bubble is less curved than the critical bubble by its definition, thus implying its velocity is larger whereas its shear is smaller. The solutions for the velocity and curvature of the convergence-limit bubble are illustrated in Fig. 8.

3. Taylor bubble

The Taylor bubble has the curvature $ζ_{T} = - 1 / 8 k$ ⁠. The solution for the Taylor bubble can be derived from Davies and Taylor⁴⁷ upon modifying the wavevector $k = 2 π / λ$ to $k = 2 β_{1} / λ$ ⁠, where β₁ is the first zero of the Bessel function J₁. See Ref. 50 for details. The Taylor bubble solutions are written as

\begin{matrix} ζ_{T} = - \frac{1}{8} k, v_{T} = \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} 2 \sqrt{\frac{2 A}{3 + 5 A}}, \\ Γ_{T} = \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} 3 \sqrt{\frac{2 A}{3 + 5 A}} . \end{matrix}

(42)

Only at N = 1 for the Taylor bubble, the second Fourier amplitude for heavy fluid is zero $Φ_{2} = 0$ ⁠. To compare the solutions to those for the convergence-limit bubble, the curvature is less curved $ζ_{T} / ζ_{c n} = 3 / 5$ ⁠, and thus, the Taylor bubble velocity is larger whereas its shear is smaller (e.g., see Fig. 8).

4. Layzer-drag bubble

The Layzer-drag bubble is defined as the bubble with the velocity $v = {(k τ)}^{- 1} {(t / τ)}^{a / 2} \sqrt{2 A / (1 + A)}$ ⁠. This special solution is called Layzer-drag bubble since Layzer's first-order approximation for bubble velocity at A = 1 applies this velocity rescaling.⁴⁸ The experiments and simulations are often compared well to the Layzer-drag solutions. These special solutions are expressed as

\begin{matrix} ζ_{L} = ζ_{L} (k, A), v_{L} = \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{\frac{2 A}{1 + A}}, \\ Γ_{L} = \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} 12 {(9 - 64 {(\frac{{\hat{ζ}}_{L}}{k})}^{2})}^{- 1} \sqrt{\frac{2 A}{1 + A}} . \end{matrix}

(43)

The curvature of the Layzer-drag bubble is cumbersome and is given as $ζ_{L} (k, A)$ ⁠. The dimensionless curvature $ζ_{L} / k$ is plotted as a function of A in Fig. 8(a). However, a simplified expression is only available at the end of the limits of the Atwood number A. At the limit of $A \to 1^{-}$ ⁠, where the two fluids have highly contrasting densities, the solutions for the Layzer-drag bubble are

ζ_{L} \approx - \frac{1}{8} k, v_{L} \approx \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}}, Γ_{L} \approx \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} \frac{3}{2} .

(44)

This result is the same solution as for the Taylor bubble. At the limit of $A \to 0^{+}$ ⁠, where the two fluids have very similar densities,

\begin{matrix} ζ_{L} \approx - \frac{1}{8} k \sqrt{9 - 4 \sqrt{3}}, v_{L} \approx \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{2 A}, \\ Γ_{L} \approx \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{6 A} . \end{matrix}

(45)

The Layzer-drag bubble is less curved for the contrasting-density fluid pair $A \to 1^{-}$ than for the similar-density fluid pairs $A \to 0^{+}$ ⁠. For a finite density ratio $A \in (0, 1)$ ⁠, the Layzer-drag bubble is related to the critical bubble and Taylor bubble as $ζ_{L} / ζ_{c r} \in (\sqrt{9 - 4 \sqrt{3}} / 3, 1 / 3)$ and $ζ_{L} / ζ_{T} \in (\sqrt{9 - 4 \sqrt{3}}, 1)$ ⁠, respectively. Similar to the previous analysis, the velocity for Layzer-drag bubble solutions is large whereas the shear is smaller than that for the Taylor bubble. The solution of Eq. (43) is presented in Fig. 8.

5. Atwood bubble

As introduced earlier, the Atwood bubble is the bubble with the fastest velocity within the family of nonlinear solutions. In order to derive the curvature for the Atwood bubble ζ_A, the characteristic equation for p4 mm is sought as

4096 A ξ_{A}^{4} + 4096 ξ_{A}^{3} + 1152 A ξ_{A}^{2} - 27 A = 0,

(46)

which is the same equation as Eq. (46) but in the infinite domain. Recall that $ξ_{A} = - ζ_{A} / k$ and $ζ_{A} < 0$ for the bubble. We find four roots by solving for Eq. (46), and only the roots that are positive in $A \in [0, 1]$ are relevant. Substitution of these roots into Eq. (36) leads to the velocity and shear expressions for the Atwood number. However, they have a complex dependence on the Atwood number A. Thus, they are illustrated in Fig. 8 without the exact analytical expression.

The fundamental scales in Eq. (34) are based on the assumption that the vertical length scale characterizing the nonlinear regime solutions for bubble becomes relevant to horizontal length scales. This emergence of another length scale can be demonstrated mathematically by deriving an invariant property of the Atwood bubble. To proceed, we solve Eq. (46) for the Atwood number A and substitute the result for the velocity in Eq. (29). Then, the invariant property of the Atwood bubble for p4 mm is written as

\frac{1}{512} v_{A}^{2} {(k τ)}^{2} ξ_{A}^{- 3} {(\frac{t}{τ})}^{- a} = {\frac{S_{k}^{4}}{S_{g}} |}_{k Z \to \infty} = 1.

(47)

This invariant property in Eq. (47) is remarkable since it implies the multiscale dynamics of nonlinear RT dynamics. More specifically, the left-hand side of Eq. (47) is composed of Atwood-bubble velocity v_A, the horizontal length scale from the initial wavelength $k = 2 π / λ$ ⁠, and the vertical length scale from the amplitude ζ_A, whereas the right-hand side is a constant for a given finite domain size kZ. Thus, relating the horizontal and vertical length scale, Eq. (47) demonstrates the multiscale feature of nonlinear RTI unlike linear RTI, which is governed by a single horizontal length scale.

Similar to the solutions for the Layzer-drag bubble, solutions at the two end limits of A can be found in a simplified form. First, for fluids with highly contrasting densities $A \to 1^{-}$ ⁠, the solutions are

\begin{matrix} ζ_{A} \approx - \frac{{8 - (1 - A)}}{64} k, v_{A} \approx \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} \frac{{16 - 3 (1 - A)}}{16}, \\ Γ_{A} \approx \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} \frac{3 {32 - 7 (1 - A)}}{64} . \end{matrix}

(48)

These results are the same solutions as for the Taylor bubble and the Layzer-drag bubble at the same limit. This convergence of the three different solutions at A = 1 is shown in Fig. 8. At the limit of $A \to 0^{+}$ ⁠, where the two fluids have very similar densities,

\begin{matrix} ζ_{A} \approx - \frac{3}{16} k A^{\frac{1}{3}}, v_{A} \approx \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(\frac{3}{2})}^{\frac{3}{2}} \sqrt{A}, \\ Γ_{A} \approx \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{6 A} . \end{matrix}

(49)

For the finite fluid density ratio $A \in (0, 1)$ ⁠, the curvatures of the Atwood bubble are related to previous solutions as $ζ_{A} / ζ_{c r} \in (0, 1 / 3), ζ_{A}, ζ_{L} \in (0, 1)$ and $ζ_{A} / ζ_{T} \in (0, 1)$ ⁠. The properties of the Atwood-bubble properties are presented in Fig. 8. As can be seen in the figure, the velocity is larger and the shear is smaller for the Atwood bubble than for the Layzer-drag and Taylor bubble.

6. Flat bubble

The flat bubble is a stagnant bubble with zero curvature such that v_f = 0 and $ζ_{f} = 0$ ⁠. The flat bubble solution is expressed as

ζ_{f} = 0, v_{f} = 0, Γ_{f} = 0.

(50)

We note that while the velocity is nonzero for both the critical bubble and the flat bubble, the curvature, and thus the shear, is zero only for the flat bubble.

7. Effect of the domain size

As discussed in Sec. III E, the effect of the finite domain size kZ on special solutions can be demonstrated by scaling Eqs. (33)–(35). The Layzer-drag and Atwood bubble solutions are illustrated with two different domain sizes kZ = 1.5 and kZ = 5 in Fig. 10. Due to the multiscale behavior of the solutions in the nonlinear regime, the curvature, velocity, and shear are scaled separately. The curvature ζ and velocity v decrease for smaller domain size, whereas the shear Γ increases. Figure 10(d) shows enhanced shear for the same velocity at the tip of the bubble. In summary, the bubble moves slower with a less curved front and larger shear in the finite-sized domain compared to one in an infinite domain.

FIG. 10.

View large Download slide

The effect of the finite domain size kZ on the Layzer-drag and Atwood bubble solutions. (a) Curvature, (b) velocity, and (c) shear are illustrated as functions of the Atwood number A. A one-to-one relation between the shear and the velocity is presented in (d). The various lines correspond to the solution, the Layzer-drag solution (gray), the Atwood solution (black), kZ = 1.5 (dashed lines), and kZ = 5 (solid lines).

H. Special solutions: Spikes

Similar to the nonlinear regime solutions for bubbles, special solutions for spikes are established in this section—the critical spike, the convergence-limit spike, the Taylor spike, the Layzer-drag spike, the Atwood spike, and the flat spike.

1. Critical spike

The critical spike defined above has zero velocity (⁠ ${\hat{v}}_{c r} = 0$ ⁠) at the critical curvature $\hat{ζ} = {\hat{ζ}}_{c r}$ ⁠. The solutions are

{\hat{ζ}}_{c r} = \frac{3}{8} k, {\hat{v}}_{c r} = 0, {\hat{Γ}}_{c r} = - \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{\frac{6 A}{1 - A}} .

(51)

We note that the critical shear is ${\hat{Γ}}_{c r} = {\hat{Γ}}_{max}$ ⁠.

2. Convergence-limit spike

The convergence-limit spike is categorized for the spike whose Fourier amplitudes in Eq. (37) follow $| Φ_{1} | > | Φ_{2} |$ at $\hat{ζ} \in ({\hat{ζ}}_{r}, {\hat{ζ}}_{c r})$ and $| {\tilde{Φ}}_{1} | > | {\tilde{Φ}}_{2} |$ at $\hat{ζ} \in ({\hat{ζ}}_{r}, {\hat{ζ}}_{c n})$ ⁠. As was shown in Fig. 9 for bubble amplitudes, the Fourier amplitudes for the lighter fluid ${\tilde{Φ}}_{1}$ and ${\tilde{Φ}}_{2}$ for spike solutions manifest the same trend as the spike curvature $\hat{ζ}$ varies. Thus, the Fourier amplitudes of light fluid are balanced $| {\tilde{Φ}}_{1} | = | {\tilde{Φ}}_{2} |$ for spike at convergence-limit curvature $\hat{ζ} = {\hat{ζ}}_{c n}$ ⁠, whereas the Fourier amplitudes of the heavy fluid are balanced for bubbles at $ζ = ζ_{c n}$ ⁠. The solutions for the convergence-limit spike at $\hat{ζ} = {\hat{ζ}}_{c n}$ are written as

\begin{matrix} {\hat{ζ}}_{c n} = \frac{5}{24} k, {\hat{v}}_{c n} = - \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} \frac{14}{3} \sqrt{\frac{10 A}{3 (45 - 53 A)}}, \\ {\hat{Γ}}_{c n} = - \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} 9 \sqrt{\frac{10 A}{3 (45 - 53 A)}} . \end{matrix}

(52)

The curvature for the convergence-limit spike has the same magnitude as that of the bubble, while the velocity and shear approach the negative infinity as $A \to 45 / 53 \approx 0.85$ due to the singularity. Also, we note that the velocity and shear become zero as $A \to 0^{+}$ ⁠. The convergence-limit curvature is related to the critical curvature as ${\hat{ζ}}_{c n} / {\hat{ζ}}_{c r} = 5 / 9$ ⁠, which implies that the convergence-limit spike is less curved than the critical spike.

3. Taylor spike

The Taylor spike is a spike with a curvature ${\hat{ζ}}_{T} = 1 / 8 k$ ⁠. This definition stems from the same analogy for the Taylor bubble. The solutions for the Taylor spike at $\hat{ζ} = {\hat{ζ}}_{T}$ are

\begin{matrix} {\hat{ζ}}_{T} = \frac{1}{8} k, {\hat{v}}_{T} = - \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} 2 \sqrt{\frac{2 A}{3 - 5 A}}, \\ {\hat{Γ}}_{T} = - \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} 3 \sqrt{\frac{2 A}{3 - 5 A}} . \end{matrix}

(53)

We note that only at N = 1 for the Taylor spike, the second Fourier amplitude for light fluid is zero, ${\tilde{Φ}}_{2} = 0$ ⁠. The curvature for the Taylor spike has the same magnitude as for the Taylor bubble curvature. The velocity ${\hat{v}}_{T}$ and shear ${\hat{Γ}}_{T}$ have finite values for $0 < A < 3 / 5 \approx 0.6$ and approach zero as $A \to 0^{+}$ ⁠. As the Atwood number is increased to the other limit $A \to 3 / 5 \approx 0.6, {\hat{v}}_{T}$ and Γ_T diverge to negative infinity (Fig. 11). Comparisons of the Taylor-spike curvature with the previous special solutions for spikes are ${\hat{ζ}}_{T} / {\hat{ζ}}_{c r} = 1 / 3$ and ${\hat{ζ}}_{T} / {\hat{ζ}}_{c n} = 3 / 5$ ⁠.

FIG. 11.

View large Download slide

Special solutions for spikes among the one-parameter family of asymptotic solutions for p4 mm. (a) Curvature, (b) velocity and (c) shear are illustrated as functions of Atwood number A. A one-to-one relation between the shear and the velocity is presented in (d). The various lines correspond to solution, the critical spike (dash-dotted lines), the convergence-limit spike (dotted lines), the Taylor spike (dashed lines), the Layzer-drag spike (gray solid lines), and the Atwood spike (black solid lines).

4. Layzer-drag spike

The Layzer-drag spike is a spike with velocity $\hat{v} = {(k τ)}^{- 1} {(t / τ)}^{a / 2} \sqrt{2 A / (1 - A)}$ and has the same origin as the Layzer-drag bubble. The solutions for the Layzer-drag spike are expressed as

\begin{matrix} {\hat{ζ}}_{L} = {\hat{ζ}}_{L} (k, A), {\hat{v}}_{L} = - \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{\frac{2 A}{1 - A}}, \\ {\hat{Γ}}_{L} = - \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} 12 {(9 - 64 {(\frac{{\hat{ζ}}_{L}}{k})}^{2})}^{- 1} \sqrt{\frac{2 A}{1 - A}} . \end{matrix}

(54)

The curvature of the Layzer-drag spike is cumbersome and is given as ${\hat{ζ}}_{L} (k, A)$ ⁠. ${\hat{ζ}}_{L} = {\hat{ζ}}_{L} (k, A)$ is also shown in Fig. 11(a). The velocity ${\hat{v}}_{L}$ and shear ${\hat{Γ}}_{L}$ of the solutions for the Layzer-drag spike approach zero as $A \to 0^{+}$ ⁠. ${\hat{v}}_{L}$ and ${\hat{Γ}}_{L}$ remain finite as negative values within the range $0 < A < 1$ and become singular, diverging to the negative infinity as $A \to 1^{-}$ ⁠. The Layzer-drag solutions are further simplified at the two limits of $A \in {0, 1}$ ⁠. At the limit of $A \to 1^{-}$ ⁠, where the two fluids have very contrasting densities, the solutions for the Layzer-drag spike are

\begin{matrix} {\hat{ζ}}_{L} \approx \frac{1}{8} k, {\hat{v}}_{L} \approx - \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{\frac{2}{1 - A}}, \\ {\hat{Γ}}_{L} \approx - \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} \frac{1}{4} {(\frac{3}{8} - \frac{{\hat{ζ}}_{L}}{k})}^{- 1} \sqrt{\frac{2}{1 - A}} . \end{matrix}

(55)

Both ${\hat{v}}_{L}$ and ${\hat{Γ}}_{L}$ at the limit $A \to 1^{-}$ become singular, with the latter being more singular than the former. At the limit of $A \to 0^{+}$ ⁠, where the two fluids have very similar densities,

\begin{matrix} {\hat{ζ}}_{L} \approx \frac{1}{8} k \sqrt{9 - 4 \sqrt{3}}, {\hat{v}}_{L} \approx - \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{2 A}, \\ {\hat{Γ}}_{L} \approx - \frac{1}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} \sqrt{6 A} . \end{matrix}

(56)

The Layzer-drag spike is more curved for the contrasting-density fluid pair $A \to 1^{-}$ than for the similar-density fluid pair $A \to 0^{+}$ ⁠. For a finite density ratio $A \in (0, 1)$ ⁠, the Layzer-drag spike is related to the critical spike and Taylor spike as ${\hat{ζ}}_{L} / {\hat{ζ}}_{c r} \in (\sqrt{9 - 4 \sqrt{3}} / 3, 1)$ and ${\hat{ζ}}_{L} / {\hat{ζ}}_{T} \in (\sqrt{9 - 4 \sqrt{3}}, 3)$ ⁠, respectively. This curvature ratio implies that the Layzer-drag velocity ${\hat{v}}_{L}$ is smaller and the Layzer-drag shear ${\hat{Γ}}_{L}$ is larger than those for the Taylor spike. The solutions in Eq. (54) are presented in Fig. 11. We note that not all solutions for the Layzer-drag spike are convergent. The ratio of the curvature for the Layzer-drag spike to that of the convergence-limit spike becomes unity ${\hat{ζ}}_{L} / {\hat{ζ}}_{c n} = 1$ at $A = 235 / 451 \approx 0.52$ ⁠. Thus, the solutions for the Layzer-drag spike are not convergent beyond A >0.52 [see, e.g., Fig. 11(a)].

5. Atwood spikes

The Atwood solutions are the fastest solutions among the family of solutions in Eq. (36). Thus, the curvature, velocity, and shear follow ${\hat{ζ}}_{A} = {\hat{ζ}}_{r}, {\hat{v}}_{A} = {\hat{v}}_{max}$ ⁠, and ${\hat{Γ}}_{A} = {\hat{Γ}}_{max}$ ⁠, respectively. The special solutions for the Atwood spike are written as

\begin{matrix} {\hat{ζ}}_{A} = κ (A) k, {\hat{v}}_{A} = - \frac{ϕ (A)}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(\frac{k}{\hat{ζ} - {\hat{ζ}}_{A}})}^{\frac{1}{2}}, \\ {\hat{Γ}}_{A} = - \frac{γ (A)}{τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(\frac{k}{\hat{ζ} - {\hat{ζ}}_{A}})}^{\frac{1}{2}}, \end{matrix}

(57)

where

\begin{matrix} κ (A) = \frac{3 (1 - \sqrt{1 - A^{2}})}{8 A}, ϕ (A) = \frac{9 {(1 - A^{2})}^{\frac{1}{4}} {(1 - \sqrt{1 - A^{2}})}^{\frac{3}{2}}}{4 A^{2}}, \\ γ (A) = \frac{3 {(1 - \sqrt{1 - A^{2}})}^{\frac{1}{2}}}{2 {(1 - A^{2})}^{\frac{1}{4}}} . \end{matrix}

(58)

For the Atwood spike, the curvature ${\hat{ζ}}_{A}$ is regular whereas the velocity ${\hat{v}}_{A}$ and the shear ${\hat{Γ}}_{A}$ are singular as shown in Fig. 11. The solutions in Eq. (58) are simplified at the two end limits of the Atwood number A. At the limit of $A \to 1^{-}$ ⁠, where the two fluids have very contrasting densities,

\begin{matrix} κ (A) \approx \frac{3}{8} (1 - \sqrt{2 (1 - A)}), ϕ (A) \approx \frac{9}{4} {2 (1 - A)}^{\frac{1}{4}}, \\ γ (A) \approx \frac{3}{2} {2 (1 - A)}^{- \frac{1}{4}} . \end{matrix}

(59)

At the limit of $A \to 0^{+}$ ⁠, where the two fluids have very similar densities,

κ (A) \approx \frac{3}{16} A, ϕ (A) \approx \frac{9}{8 \sqrt{2}} A, γ (A) \approx \frac{3}{2 \sqrt{2}} A .

(60)

The Atwood-spike is less curved than the Layzer-drag spike in $A \in [0, a)$ ⁠. At A = 1, the curvature ${\hat{ζ}}_{A}$ becomes the same curvature as the Layzer-drag curvature ${\hat{ζ}}_{L}$ and the critical curvature ${\hat{ζ}}_{c r}$ ⁠. The Atwood spike curvature increases as the Atwood number A is increased. Its curvature ${\hat{ζ}}_{A}$ intersects with the curvature for the Taylor spike ${\hat{ζ}}_{A} = {\hat{ζ}}_{T}$ at $A = 3 / 5 \approx 0.6$ and with the curvature for convergence-limit spike ${\hat{ζ}}_{A} = {\hat{ζ}}_{c n}$ at $A = 45 / 53 \approx 0.85$ ⁠. Then, ${\hat{ζ}}_{A} = {\hat{ζ}}_{c n}$ implies that the solutions for the Atwood spike are no longer convergent for A >0.85 [see, e.g., Fig. 11(a)].

Especially noteworthy is that while the curvatures of both the Atwood spike and the Atwood bubble are finite, the curvature is more curved for the Atwood spike than for the Atwood bubble for fluid pairs with highly contrasting densities whereas the opposite is true for fluid pairs with similar densities.

6. Flat spike

A flat spike is a stagnant spike (⁠ ${\hat{v}}_{f} = 0$ ⁠) with zero curvature (⁠ ${\hat{ζ}}_{f} = 0$ ⁠). The solutions for the flat spike are expressed as

{\hat{ζ}}_{f} = 0, {\hat{v}}_{f} = 0, {\hat{Γ}}_{f} = 0.

(61)

We note that while the velocity is zero for both the critical spike and the flat spike, the curvature, and thus the shear, is zero only for the flat spike.

7. Effect of the domain size

Similar to the case of bubbles' special solutions, the effect of finite domain size kZ can be shown by multiplying Eqs. (33)–(35) to the solutions in the infinite-sized domain. The Layzer-drag and Atwood spike solutions are presented with two different domain sizes kZ = 1.5 and kZ = 5 in Fig. 12. The curvature $\hat{ζ}$ and velocity $\hat{v}$ decrease for smaller domain size, whereas the shear $\hat{Γ}$ increases. Thus, the spike moves slower with a less curved front and has larger shear in the finite domain compared to the infinite domain. Unlike bubbles, spikes have a singularity at the curvature $\hat{ζ} = {\hat{ζ}}_{A}$ ⁠. The velocity and shear for the Atwood spike diverge, so they do not appear in Figs. 12(b) and 12(c). In Fig. 12(d), however, the shear is increased for the same velocity at the spike tip with a reduced finite domain size for both the Layzer-drag and Atwood spike solutions.

FIG. 12.

View large Download slide

The effect of the finite domain size kZ on the Layzer-drag and Atwood spike solutions. (a) Curvature, (b) velocity, and (c) shear are illustrated as functions of Atwood number A. A one-to-one relation between the shear and the velocity is presented in (d). The various lines correspond to the solution, the Layzer-drag solution (gray), the Atwood solution (black), kZ = 1.5 (dashed lines), and kZ = 5 (solid lines).

I. Transition to mixing

Along the evolution of RTI, we demonstrated the properties of large coherent structures of RTI characterized by the dominant scales. In the linear regime, the horizontal length scale, the fundamental wavelength λ, serves as a single dominant length scale. In contrast, the two length scales, the amplitude ζ and the wavelength λ, become the dominant length scales characterizing nonlinear regime solutions, thus manifesting multiscale behavior. After the scale-dependent linear and nonlinear regimes, it is known that RTI transits into a scale-invariant regime where the vertical scale grows in a self-similar manner with time.

Classical understanding of the transition of the RT dynamics to self-similar mixing is based on a merger mechanism ascribed to the growth of horizontal scales.⁴⁹ However, prior research on the momentum model confirms the possibility of the transition of RT dynamics to self-similar mixing given that the vertical length scale becomes the dominant scale in the system.¹² Also, several experiments of RTI with high Reynolds number substantiate this occurrence of the transition via the amplitude dominance mechanism.¹ Since the group theory approach in this study demonstrates the self-similar RT dynamics within the framework of the momentum model, the transition from scale-dependent dynamics to self-similar mixing can be associated with the multiscale characteristics of nonlinear dynamics of RTI. Here, we discuss how the transition from the nonlinear regime to the mixing regime occurs in terms of a scaling analysis for bubbles and spikes.

1. Bubbles

In a laboratory framework, the interface is described as $Z (r, t) = z_{0} (r, t) + z^{*}$ ⁠, where z₀ is the location of the free surface and $z^{*}$ is the interface distortion in the reference frame of the tip in Eq. (13), while r and t denote a vector of arbitrary position on the interface and time, respectively. Then, the velocity is $v = \dot{z_{0}}$ ⁠. The Atwood-bubble solutions are the fastest solutions within the family of asymptotic solutions in the nonlinear regime. From Eq. (36), the velocity $v_{a} = v_{max}$ at the curvature $ζ = ζ_{A}$ of the Atwood bubble satisfies

| ζ_{A} | \sim k \sim λ^{- 1}, | v_{A} | \sim \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} = {(\frac{g}{k})}^{\frac{1}{2}} .

(62)

Then,

| v | \sim \sqrt{G λ t^{a}}, | z_{0} | \sim t \sqrt{G λ t^{a}} .

(63)

Following Eq. (47), the amplitude z₀ of the bubble grows to a vertical length scale comparable to the horizontal length scale λ. For $| z_{0} | > λ$ ⁠, the vertical length scale becomes the dominant length scale characterizing the RTI. From the studies of the momentum model,¹² the rate of loss of the specific momentum (i.e., the drag force per unit mass) scales as $\sim v^{2} / | z_{0} |$ ⁠. Applying the transition of the dominant length scale from λ to $| z_{0} |$ shows the reduction of the drag force from $\sim v^{2} / λ$ to $\sim v^{2} / | z_{0} |$ ⁠. This reduced drag force leads to the acceleration of the bubble and thus provokes the transition of RTI dynamics from scale-dependent dynamics to self-similar mixing. By substituting λ in Eq. (63) to $| z_{0} |$ ⁠, the scalings of the velocity $| v |$ and amplitude $| z_{0} |$ imply

| v | \sim G t^{a + 1}, | z_{0} | \sim G t^{a + 2} .

(64)

2. Spikes

Analogous to bubbles, the transition to the self-similar mixing of spikes is addressed with the Atwood spike. The solutions for the Atwood spikes follow

(65)

| {\hat{ζ}}_{A} | \sim k \sim λ^{- 1},

(65a)

| {\hat{v}}_{A} | \sim \frac{1}{k τ} {(\frac{t}{τ})}^{\frac{a}{2}} {(\frac{k}{\hat{ζ} - {\hat{ζ}}_{A}})}^{\frac{1}{2}} = {(\frac{g}{k})}^{\frac{1}{2}} {(\frac{k}{\hat{ζ} - {\hat{ζ}}_{A}})}^{\frac{1}{2}} .

(65b)

Then,

| \hat{v} | \sim \sqrt{G λ t^{a}}, | z_{0} | \sim t \sqrt{G λ t^{a}} .

(66)

The amplitude z₀ of the spikes grows to a vertical length scale comparable to the horizontal length scale λ. For $| z_{0} | > λ$ ⁠, the vertical length scale becomes the dominant length scale characterizing the RTI. From the studies of the momentum model,¹² the rate of loss of specific momentum (i.e., the drag force per unit mass) scales as $\sim {\hat{v}}^{2} / | z_{0} |$ ⁠. Applying the transition of the dominant length scale from λ to $| z_{0} |$ shows the reduction of the drag force from $\sim {\hat{v}}^{2} / λ$ to $\sim {\hat{v}}^{2} / | z_{0} |$ ⁠. This reduced drag force leads to the acceleration of the spike and thus provokes the transition of RTI dynamics from scale-dependent to self-similar mixing. By substituting λ in Eq. (66) to $| z_{0} |$ ⁠, the scalings of the velocity $| \hat{v} |$ and amplitude $| z_{0} |$ imply

| \hat{v} | \sim G t^{a + 1}, | z_{0} | \sim G t^{a + 2} .

(67)

Since the spikes are singular, the growth of the amplitude z₀ can be faster than the analytical solution in Eq. (65). The fast growth in amplitude indicates the dominance of the vertical length scale as RTI evolves. Thus, for $| z_{0} | ≫ λ$ ⁠, the drag force is reduced, thus leading to a self-similar mixing in Eq. (67). For fluids with very contrasting densities, $A \to 1^{-}$ ⁠, the Atwood spike is highly singular and can only be self-similar.

J. Interfacial dynamics

The interface for 3D RTI with group p4 mm (square symmetry) is presented in Fig. 13(a). The first Fourier harmonic is retained to provide the dominant mode shape. Also, the velocity field is achieved from the fluid potentials in Eq. (12). The gradient of the fluid potential yields the velocity field as shown in Figs. 13(b) and 13(c) for a bubble and a spike, respectively. In Figs. 13(b) and 13(c), the velocity near the interface is large in magnitude and decreases as we move away from the interface. This observation indicates that the RTI is an interfacial phenomenon.

FIG. 13.

View large Download slide

Schematic view of the RTI for p4 mm subspace group. (a) Interface contour of the first harmonic. Instantaneous velocity fields at the tips of the large coherent structures: (b) Atwood bubble velocity with the curvature $- ζ / k = 0.125$ and Atwood number A = 0.98, and (c) Atwood spike velocity with the curvature $\hat{ζ} / k = 0.153$ and Atwood number A = 0.7. Velocity fields are evaluated as the sum of the first and the second Fourier harmonics.

We note that the working fluid in the present study is idealized so that viscosity, compressibility, surface tension, and mass flux are not considered. Nonetheless, a qualitative comparison of the flow field corroborates the validity of the current approach. More specifically, our method for time-varying RTI with group velocity successfully captures not only the asymptotic behavior of the large coherent motions but also the bulk fluid motion in the region away from the interface. This capability to predict the flow field confirms that RTI results in an interfacial flow.

Moreover, we observe the velocity difference between the two fluids in the tangential direction to the interface. This difference in the tangential velocity gradient across the interface is quantified as the shear Γ in this study. According to Fig. 13(b), the shear increases as we move away from the tip of the large coherent structure along the interface. Note that our theory is valid near the tips; however, the current approach provides a better understanding of the generation of secondary instabilities in RTI.

Previous studies agree that the secondary instability is triggered due to shear across the interface during RTI evolution.^1,22,30,32 This secondary instability is in small-scale vortical structures and invokes shear-driven Kelvin–Helmholtz instability (KHI). The shear function Γ measures the intensity of the tangential velocity difference across the interface. Thus, Γ estimates the growth rate of KHI for bubbles and spikes as $ω_{KHI} \sim Γ$ and ${\hat{ω}}_{KHI} \sim \hat{Γ}$ ⁠, respectively. Note that the nonlinear solutions for spikes are singular. This singularity $\hat{Γ} \to - \infty$ as $\hat{ζ} \to {\hat{ζ}}_{A}$ implies the possibility of larger shear compared to the analytically predicted dynamics. On the other hand, nonlinear solutions for bubbles are regular. Since the growth rate of KHI is proportional to the shear $ω_{KHI} \sim Γ$ ⁠, small-scale vortical structures are expected to be more pronounced for spikes than for bubbles.

IV. DISCUSSION AND CONCLUSION

Our study investigates the dynamics of scale-dependent RTI in 3D periodic flow with square symmetry under time-varying acceleration in a finite-sized domain. The acceleration is acting normal to the periodic plane and follows a power law in time. We consider the acceleration exponents $a > - 2$ for RT-type dynamics. Group theory is applied for defining the irreducible representations. A boundary value problem is linked to the governing equations and is reduced to a dynamical system, which is solved with the expansions of velocity potentials. The solutions for the coherent structures are sought for both the linear and nonlinear regimes. A thorough investigation to reveal the characteristics of the solutions is conducted.

Five core findings emerge in this study. First, the framework established by applying group theory provides a systematic approach for both the linear and nonlinear regimes. Our approach respects first principles, and the resulting dynamical systems are solved with an appropriate truncation of the expansions at a desired order. The family of asymptotic solutions in both the linear and nonlinear regimes are derived. This universal form of the solutions is the second significant result of the current study. The identified form of the solutions covers a broad range of parameters, including acceleration strengths and exponents, Atwood numbers, and the size of the finite domain. Moreover, special solutions in the nonlinear regimes are established from a one-parameter family solution. Third, finite-domain-size effects on the solutions in both the linear and nonlinear regimes are demonstrated. The finite domain size reduces the penetration velocity and front curvature, whereas it increases the shear. Due to the multiscale character of the dynamics, characteristic scales in the nonlinear regime are separately scaled with the domain size in contrast to the linear regime. Fourth, the interplay of the length scales is revealed. The solutions' single-scale and multiscale characters are shown in the linear and nonlinear regimes, respectively. Furthermore, the emergence of the vertical length scale in the nonlinear regime sheds light on the transition to the mixing regime. Finally, the velocity field predicted by the current approach manifests the interfacial character of RTI. This interfaciality corroborates the validity and completeness of the present method.

The results encompass the cases with a constant acceleration in an infinite domain. Thus, the connection between the present and the prior studies is strongly established. An excellent qualitative agreement is observed: stability of the solutions at the tips of large coherent structures, development of small-scale vortical structures due to high shear, interfacial behavior rather than volumetric, Atwood number dependency on the nonlinear solution, and eventual transition to the mixing regime.

Interpretation of the findings suggests remarkable impacts to the phenomena at astrophysical and at molecular scales. The mixing of the two liquids is only pronounced at the vicinity of the interface, and the interface of coherent structures is affected by the time-dependent acceleration. Also, the dynamics are influenced by the deterministic initial condition. Thus, applications, such as ICF, can be better controlled by applying variable acceleration, and similarly, a particular type of supernova where RTI is the dominant mechanism can be explored in a more systematic way. The current research strengthens the capability to diagnose and predict the RTI and RT mixing events.

In the present work, we focused on the RTI subject to acceleration with power-law time dependence and with the exponent of the power-law greater than −2. We address in the future the study of RM dynamics with the exponent of the acceleration's power law smaller than −2, and the RT-to-RM transition when the exponent of the acceleration power-law is equal to −2. We also address in the future the study of linear and nonlinear RT and RM dynamics with the scale-dependent acceleration, such as exponential, smeared delta or Gaussian.

In summary, the highlighted results present a comprehensive reference and benchmark for experiments and simulations. Derived solutions at each stage of the RTI evolution and the effect of the various parameters as well as the interplay between the dominant length scales serve as excellent diagnostic measures for future studies. We revealed the connection between the interface dynamics and its morphology. We also find the features of universality in the dynamics in the linear and nonlinear regimes. Our theory can be extended to more complicated problems by permitting a finite mass flow rate at the interface; an inspection of the different acceleration parameter regimes preserving the same framework can explore the dynamics of RM-type instability. Therefore, our work provides the groundwork that improves our knowledge of RTI dynamics and lays a strong foundation for the benchmark of 3D time-varying RTI in finite-sized domains.

AUTHORS' CONTRIBUTIONS

H.H. contributed to investigation, analysis, and preparation of the manuscript. S.I.A. contributed to conceptualization, methodology investigations, resources, and supervision. W.H.R.C. and S.S.J. contributed to discussion.

ACKNOWLEDGMENTS

H. Hwang was supported by the Predictive Science Academic Alliance Program 3 (PSAAP3) and P. Franklin and Caroline M. Johnson Fellowship at the Stanford University. W. H. R. Chan was funded by a National Science Scholarship from the Agency for Science, Technology, and Research in Singapore and P. Franklin and Caroline M. Johnson Fellowship at the Stanford University. S. S. Jain was supported by P. Franklin and Caroline M. Johnson Fellowship from the Stanford University. S. I. Abarzhi was funded by the University of Western Australia and the National Science Foundation, USA.

DATA AVAILABILITY

The data that support the findings of this study are available within the article and from the corresponding author upon reasonable request.

References

1.

E. E.

Meshkov

and

S. I.

Abarzhi

, “

Group theory and jelly's experiment of Rayleigh–Taylor instability and Rayleigh–Taylor interfacial mixing

,”

Fluid Dyn. Res.

51

,

065502

(

2019

).

https://doi.org/10.1088/1873-7005/ab3e83

Google Scholar

Crossref

2.

L.

Rayleigh

, “

Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density

,”

Proc. London Math. Soc.

14

,

170

–

177

(

1882

).

Google Scholar

3.

G. I.

Taylor

, “

The instability of liquid surfaces when accelerated in a direction perpendicular to their planes

,”

Proc. R. Soc. London, Ser. A

201

,

192

–

196

(

1950

).

Google Scholar

4.

R. D.

Richtmyer

, “

Taylor instability in shock acceleration of compressible fluids

,”

Commun. Pure Appl. Math.

13

,

297

–

319

(

1960

).

https://doi.org/10.1002/cpa.3160130207

Google Scholar

Crossref

5.

E. E.

Meshkov

, “

Instability of the interface of two gases accelerated by a shock wave

,”

Fluid Dyn.

4

,

101

–

104

(

1972

).

https://doi.org/10.1007/BF01015969

Google Scholar

Crossref

6.

S. P.

Regan

,

R.

Epstein

,

B. A.

Hammel

,

L. J.

Suter

,

J.

Ralph

,

H.

Scott

,

M. A.

Barrios

,

D. K.

Bradley

,

D. A.

Callahan

,

C.

Cerjan

,

G. W.

Collins

,

S. N.

Dixit

,

T.

Doeppner

,

M. J.

Edwards

,

D. R.

Farley

,

S.

Glenn

,

S. H.

Glenzer

,

I. E.

Golovkin

,

S. W.

Haan

,

A.

Hamza

,

D. G.

Hicks

,

N.

Izumi

,

J. D.

Kilkenny

,

J. L.

Kline

,

G. A.

Kyrala

,

O. L.

Landen

,

T.

Ma

,

J. J.

MacFarlane

,

R. C.

Mancini

,

R. L.

McCrory

,

N. B.

Meezan

,

D. D.

Meyerhofer

,

A.

Nikroo

,

K. J.

Peterson

,

T. C.

Sangster

,

P.

Springer

, and

R. P. J.

Town

, “

Hot-spot mix in ignition-scale implosions on the NIF

,”

Phys. Plasmas

19

,

056307

(

2012

).

https://doi.org/10.1063/1.3694057

Google Scholar

Crossref

7.

S. I.

Abarzhi

,

A. K.

Bhowmick

,

A.

Naveh

,

A.

Pandian

,

N. C.

Swisher

,

R. F.

Stellingwerf

, and

W. D.

Arnett

, “

Supernova, nuclear synthesis, fluid instabilities and mixing

,”

Proc. Natl. Acad. Sci. U. S. A.

116

,

18184

(

2019

).

https://doi.org/10.1073/pnas.1714502115

Google Scholar

Crossref

PubMed

8.

G. S.

Novak

,

J. P.

Ostriker

, and

L.

Ciotti

, “

Feedback from central black holes in elliptical galaxies: Two-dimensional models compared to one-dimensional models

,”

Astrophys. J.

737

,

26

(

2011

).

https://doi.org/10.1088/0004-637X/737/1/26

Google Scholar

Crossref

9.

R. P. J.

Town

and

A. R.

Bell

, “

Three-dimensional simulations of the implosion of inertial confinement fusion targets

,”

Phys. Rev. Lett.

67

,

1863

(

1991

).

https://doi.org/10.1103/PhysRevLett.67.1863

Google Scholar

Crossref

PubMed

10.

S. I.

Abarzhi

,

K.

Nishihara

, and

R.

Rosner

, “

A multi-scale character of the large-scale coherent dynamics in the Rayleigh–Taylor instability

,”

Phys. Rev. E

73

,

036310

(

2006

).

https://doi.org/10.1103/PhysRevE.73.036310

Google Scholar

Crossref

11.

S. I.

Abarzhi

, “

Coherent structures and pattern formation in Rayleigh–Taylor turbulent mixing

,”

Phys. Scr.

78

,

015401

(

2008

).

https://doi.org/10.1088/0031-8949/78/01/015401

Google Scholar

Crossref

12.

S. I.

Abarzhi

, “

Review of theoretical modeling approaches of Rayleigh–Taylor instabilities and turbulent mixing

,”

Philos. Trans. R. Soc. A

368

,

1809

–

1828

(

2010

).

https://doi.org/10.1098/rsta.2010.0020

Google Scholar

Crossref

13.

A. W.

Cook

,

W.

Cabot

, and

P. L.

Miller

, “

The mixing transition in Rayleigh–Taylor instability

,”

J. Fluid Mech.

511

,

333

–

362

(

2004

).

https://doi.org/10.1017/S0022112004009681

Google Scholar

Crossref

14.

D. H.

Olson

and

J. W.

Jacobs

, “

Experimental study of Rayleigh–Taylor instability with a complex initial perturbation

,”

Phys. Fluids

21

,

034103

(

2009

).

https://doi.org/10.1063/1.3085811

Google Scholar

Crossref

15.

S. I.

Abarzhi

and

W. A.

Goddard

, “

Interfaces and mixing: Nonequilibrium transport across the scales

,”

Proc. Natl. Acad. Sci. U. S. A.

116

,

18171

(

2019

).

https://doi.org/10.1073/pnas.1818855116

Google Scholar

Crossref

PubMed

16.

S. I.

Abarzhi

,

S.

Gauthier

, and

K. R.

Sreenivasan

, “

Turbulent mixing and beyond: Non-equilibrium processes from atomistic to astrophysical scales

,”

Philos. Trans. R. Soc. A

371

,

20120435

(

2013

).

https://doi.org/10.1098/rsta.2012.0435

Google Scholar

Crossref

17.

H.

Hwang

and

S. I.

Abarzhi

,

Two-Dimensional Scale-Dependent Rayleigh–Taylor Dynamics with Variable Acceleration in a Finite-Size Domain

(

Center for Turbulence Research Annual Research Briefs

,

Stanford University

,

2020

), pp.

371

–

382

.

Google Scholar

18.

H.

Hwang

,

W. H. R.

Chan

,

S. S.

Jain

, and

S. I.

Abarzhi

,

Linear Two-Dimensional Rayleigh–Taylor Dynamics with Variable Acceleration: Comparisons Between Theory and Simulations

(

Center for Turbulence Research

,

Stanford University

,

2020

), pp.

409

–

420

.

Google Scholar

19.

N. C.

Swisher

,

C.

Kuranz

,

W. D.

Arnettt

,

O.

Hurricane

,

H.

Robey

,

B. A.

Remington

, and

S. I.

Abarzhi

, “

Rayleigh–Taylor mixing in supernova experiments

,”

Phys. Plasmas

22

,

102707

(

2015

).

https://doi.org/10.1063/1.4931927

Google Scholar

Crossref

20.

A.

Naveh

,

D. L.

Hill

,

M. T.

Matthews

, and

S. I.

Abarzhi

, “

Early- and late-time evolution of Rayleigh–Taylor instability in a finite-sized domain by means of group theory analysis

,”

Fluid Dyn. Res.

52

,

025504

(

2020

).

https://doi.org/10.1088/1873-7005/ab693d

Google Scholar

Crossref

21.

H. F.

Robey

,

Y.

Zhou

,

A. C.

Buckingham

,

P.

Keiter

,

B. A.

Remington

, and

R. P.

Drake

, “

The time scale for the transition to turbulence in a high Reynolds number, accelerated flow

,”

Phys. Plasmas

10

,

614

–

622

(

2003

).

https://doi.org/10.1063/1.1534584

Google Scholar

Crossref

22.

E. E.

Meshkov

, “

Some peculiar features of hydrodynamic instability development

,”

Philos. Trans. R. Soc. A

371

,

20120288

(

2013

).

https://doi.org/10.1098/rsta.2012.0288

Google Scholar

Crossref

23.

A. W.

Cook

and

P. E.

Dimotakis

, “

Transition stages of Rayleigh–Taylor instability between miscible fluids

,”

J. Fluid Mech.

443

,

69

–

99

(

2001

).

https://doi.org/10.1017/S0022112001005377

Google Scholar

Crossref

24.

D. L.

Youngs

, “

The density ratio dependence of self-similar Rayleigh–Taylor mixing

,”

Philos. Trans. R. Soc. A

371

,

20120173

(

2013

).

https://doi.org/10.1098/rsta.2012.0173

Google Scholar

Crossref

25.

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

).

Google Scholar

26.

A.

Hamzehloo

,

P.

Bartholomew

, and

S.

Laizet

, “

Direct numerical simulations of incompressible Rayleigh–Taylor instabilities at low and medium Atwood numbers

,”

Phys. Fluids

33

,

054114

(

2021

).

https://doi.org/10.1063/5.0049867

Google Scholar

Crossref

27.

L.

Yang

,

C.

Shu

,

Z.

Chen

,

G.

Hou

, and

Y.

Wang

, “

An improved multiphase lattice Boltzmann flux solver for the simulation of incompressible flow with large density ratio and complex interface

,”

Phys. Fluids

33

,

033306

(

2021

).

https://doi.org/10.1063/5.0038617

Google Scholar

Crossref

28.

T.

Lyubimova

,

A.

Vorobev

, and

S.

Prokopev

, “

Rayleigh–Taylor instability of a miscible interface in a confined domain

,”

Phys. Fluids

31

,

014104

(

2019

).

https://doi.org/10.1063/1.5064547

Google Scholar

Crossref

29.

Y.

Li

,

R.

Samtaney

, and

V.

Wheatley

, “

The Richtmyer–Meshkov instability of a double-layer interface in convergent geometry with magnetohydrodynamics

,”

Matter Radiat. Extremes

3

,

207

–

218

(

2018

).

https://doi.org/10.1016/j.mre.2018.01.003

Google Scholar

Crossref

30.

K.

Kadau

,

J. L.

Barber

,

T. C.

Germann

,

B. L.

Holian

, and

B. J.

Alder

, “

Atomistic methods in fluid simulation

,”

Philos. Trans. R. Soc. A

368

,

1547

–

1560

(

2010

).

https://doi.org/10.1098/rsta.2009.0218

Google Scholar

Crossref

31.

J.

Ding

,

P.

Sun

,

S.

Huang

, and

X.

Luo

, “

Single- and dual-mode Rayleigh–Taylor instability at microscopic scale

,”

Phys. Fluids

33

,

042102

(

2021

).

https://doi.org/10.1063/5.0042505

Google Scholar

Crossref

32.

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

).

https://doi.org/10.1063/1.4733396

Google Scholar

Crossref

33.

J.

Glimm

,

D. H.

Sharp

,

T.

Kaman

, and

H.

Lim

, “

New directions for Rayleigh–Taylor mixing

,”

Philos. Trans. R. Soc. A

371

,

20120183

(

2013

).

https://doi.org/10.1098/rsta.2012.0183

Google Scholar

Crossref

34.

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

).

Google Scholar

35.

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

).

Google Scholar

36.

G.

Terrones

and

T.

Heberling

, “

Rayleigh–Taylor instability at spherical interfaces between viscous fluids: The fluid/fluid interface

,”

Phys. Fluids

32

,

094105

(

2020

).

https://doi.org/10.1063/5.0018601

Google Scholar

Crossref

37.

G. G. H.

Wolf

, “

Dynamic stabilization of the Rayleigh–Taylor instability of miscible liquids and the related ‘frozen waves

,’”

Phys. Fluids

30

,

021701

(

2018

).

https://doi.org/10.1063/1.5017846

Google Scholar

Crossref

38.

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

).

https://doi.org/10.1063/5.0031179

Google Scholar

Crossref

39.

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

).

https://doi.org/10.1063/1.5138978

Google Scholar

Crossref

40.

S. I.

Abarzhi

and

K. C.

Williams

, “

Scale-dependent Rayleigh–Taylor dynamics with variable acceleration by group theory approach

,”

Phys. Plasmas

27

,

072107

(

2020

).

https://doi.org/10.1063/5.0012035

Google Scholar

Crossref

41.

J. W.

Jacobs

and

I.

Catton

, “

Three-dimensional Rayleigh–Taylor instability part 2. Experiment

,”

J. Fluid Mech.

187

,

353

–

371

(

1988

).

https://doi.org/10.1017/S0022112088000461

Google Scholar

Crossref

42.

D. L.

Hill

,

A. K.

Bhowmick

,

D. V.

Ilyin

, and

S. I.

Abarzhi

, “

Group theory analysis of early-time dynamics of Rayleigh–Taylor instability with time varying acceleration

,”

Phys. Rev. Fluids

4

,

063905

(

2019

).

https://doi.org/10.1103/PhysRevFluids.4.063905

Google Scholar

Crossref

43.

L. D.

Landau

and

E. M.

Lifshitz

,

Course of Theoretical Physics, Vol.6: Fluid Mechanics

, 2nd ed. (

Pergamon Press

,

New York

,

1987

), pp.

336

–

343

.

Google Scholar

44.

P. R.

Garabedian

, “

On steady-state bubbles generated by Taylor instability

,”

Proc. R. Soc. A

241

,

423

(

1957

).

Google Scholar

45.

S. I.

Abarzhi

and

N. A.

Inogamov

, “

Stationary solutions in the Rayleigh–Taylor instability for spatially periodic flow

,”

J. Exp. Theor. Phys.

80

,

132

–

143

(

1995

).

Google Scholar

46.

S. I.

Abarzhi

, “

Stationary solution of the Rayleigh–Taylor instability for spatially periodic flows: Questions of uniqueness, dimensionality, and universality

,”

J. Exp. Theor. Phys.

83

,

1012

–

1026

(

1996

).

Google Scholar

47.

R. M.

Davies

and

G. I.

Taylor

, “

The mechanics of large bubbles rising through extended liquids and through liquids in tubes

,”

Proc. R. Soc. A

200

,

375

–

390

(

1950

).

Google Scholar

48.

D.

Layzer

, “

On the instability of superposed fluids in a gravitational field

,”

Astrophys. J.

122

,

1

(

1955

).

https://doi.org/10.1086/146048

Google Scholar

Crossref

49.

U.

Alon

,

J.

Hecht

,

D.

Offer

, and

D.

Shvarts

, “

Power-laws and similarity of Rayleigh–Taylor and Richtmyer–Meshkov mixing fronts at all density ratios

,”

Phys. Rev. Lett.

74

,

534

(

1995

).

https://doi.org/10.1103/PhysRevLett.74.534

Google Scholar

Crossref

PubMed

50.

S. I.

Abarzhi

, “

Length scale for bubble problem in Rayleigh–Taylor instability

,”

Phys. Fluids

11

,

940

–

942

(

1999

).

https://doi.org/10.1063/1.869964

Google Scholar

Crossref

2021

Author(s)

Scale-dependent Rayleigh–Taylor dynamics with variable acceleration in a finite-sized domain for three-dimensional flows

I. INTRODUCTION

II. METHODS

A. Boundary value problem setup

B. Time-dependent acceleration

C. Group theory

D. Fourier series expansions

III. RESULTS

A. Dynamical system

B. The shear function

C. Early time dynamics

D. Late-time dynamics

1. One-parameter family of solutions

2. Asymptotic solutions

E. Universal nonlinear scaling

F. Singularity of the solutions for RT spikes

G. Special solutions: Bubbles

1. Critical bubble

2. Convergence-limit bubble

3. Taylor bubble

4. Layzer-drag bubble

5. Atwood bubble

6. Flat bubble

7. Effect of the domain size

H. Special solutions: Spikes

1. Critical spike

2. Convergence-limit spike

3. Taylor spike

4. Layzer-drag spike

5. Atwood spikes

6. Flat spike

7. Effect of the domain size

I. Transition to mixing

1. Bubbles

2. Spikes

J. Interfacial dynamics

IV. DISCUSSION AND CONCLUSION

AUTHORS' CONTRIBUTIONS

ACKNOWLEDGMENTS

DATA AVAILABILITY

References

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only