On Rayleigh–Taylor Dynamics

Abstract

In this work, we theoretically and numerically investigate Rayleigh–Taylor dynamics with constant acceleration. On the side of theory, we employ the group theory approach to directly link the governing equations to the momentum model, and to precisely derive the buoyancy and drag parameters for the bubble and spike in the linear, nonlinear, and mixing regimes. On the side of simulations, we analyze numerical data on Rayleigh–Taylor mixing by applying independent self-similar processes associated with the growth of the bubble amplitude and with the bubble merger. Based on the obtained results, we reveal the constituents governing Rayleigh–Taylor dynamics in the linear, nonlinear, and mixing regimes. We outline the implications of our considerations for experiments in plasmas, including inertial confinement fusion.

1. Introduction

Rayleigh–Taylor instability (RTI) plays an important role in a broad range of plasma processes in nature and technology under conditions of high and low energy densities [1,2,3,4,5,6]. Examples include the abundance of chemical elements in supernova remnants, the influence of irregularities in the Earth’s ionosphere on regional climate change, the formation of hot spots in inertial confinement fusion, and the efficiency of plasma thrusters [7,8,9,10,11,12]. While realistic plasma processes are necessarily electro- and magneto-dynamic, the control of Rayleigh–Taylor (RT) phenomena often require a better understanding of the fluid dynamic aspects [4]. In this work, we study, both analytically and numerically, RTI driven by constant acceleration in neutral plasma (fluid). We employ the group theory approach and the analysis of numerical data of Front Tracking simulations [4]. We reveal independent driving constituents that the RT flow possesses in the scale-dependent and self-similar regimes.

RT unstable interfaces and interfacial mixing play an important role in plasma processes at astrophysical and atomic scales [4]. In core-collapse abarzhi-2019, RT mixing of materials of the progenitor star enable conditions for synthesis of heavy mass chemical elements [7,13,14,15,16,17]. In the Sun, RT mixing is responsible for the accelerated ejection of matter in downdrafts propagating deeply inside of the convection zone [18]. In the inertial confinement fusion, the unstable interface between the hotter plasma, the colder plasma, and the ensued interfacial mixing influence the formation and shape of the hot spot [9,10,16,19,20,21,22,23,24]. In the Earth’s ionosphere, RT unstable plasma structures are directly relevant to the climate change on regional scales [25]. In plasma processes in technology, RT unstable interfaces are inherent in light–matter interaction, nanofabrication, and in the plasma discharges formed in and interacting with liquids [26,27,28,29].

In even vastly distinct physical conditions, RT flows exhibit similar features in their evolution [4,30,31,32]. RTI develops when fluids of different densities are accelerated against their density gradients [1,2]. Small perturbations at the fluid interface grow quickly. The flow proceeds to the nonlinear regime. The interface is transformed to a composition of small-scale shear-driven vortical structures and a large-scale coherent structure, with the heavy (light) fluid penetrating the light (heavy) fluid in bubbles (spikes) [3,4]. Eventually, the flow transitions to a state of intense interfacial mixing, where the amplitude of the mixing layer increases quadratically with time [5]. RT mixing is self-similar, anisotropic, heterogeneous, and sensitive to deterministic conditions [4,33].

RTI and RT mixing are challenging to study in experiment, in theory, and in simulation [4]. RT experiments impose tight requirements on the flow implementation, diagnostics, and control [6,34]. RT simulations need to track unstable interfaces, capture small-scale processes, and ensure large spans of spatial and temporal scales [31,35]. In RT theory, we have to solve a boundary value problem at the unstable freely evolving interface, identify asymptotic solutions, and capture symmetries of the singular and ill-posed RT dynamics [4].

Significant success has recently been achieved in theory, simulation, and experiments of RT dynamics [6]. In particular, the group theory approach has identified the asymptotic solutions and the invariant forms of RT dynamics in the scale-dependent and self-similar regimes, and has explained the experimental observation that RT mixing may keep order even at high Reynolds numbers [6,34,36]. The Front Tracking simulations and the associated merger model of the numerical data have provided experiments with accurate practical mechanisms for the evolution of RT mixing and found that numerical values of the mixing zone growth agreed with experiments [5].

In RT flow, two key constituents govern the specific (per unit mass) dynamics of a fluid parcel: the rate of momentum gain (a buoyant force) and the rate of momentum loss (an effective dissipation force) [4,5,37]. The buoyant force accelerates the fluid parcel, whereas the effective dissipation decelerates it and drags it from free fall [5,33]. We need to better understand whether these buoyancy and drag constituents can be treated as independent processes in the scale-dependent linear regime, scale-dependent nonlinear regime, and self-similar mixing regime. We have to examine whether the group theory approach and the analysis of data of the Front Tracking simulations are consistent with one another [5,33].

In this work, we investigate, both theoretically and numerically, RT dynamics with constant acceleration. In theoretical analysis, we employ the group theory approach to directly link the governing equations to the momentum model and to precisely derive the buoyancy and drag parameters for the bubble and spike in the linear, nonlinear, and mixing regimes [38]. In numerical analysis, we analyze 2D experimental and 3D simulation studies of RT mixing. We construct 2D and 3D merger trees illustrating the merger process. We find that predictions made by the merger model are in agreement with physical experiments and are better than previously reported. Based on our analytical and numerical results, we reveal that the buoyant force and the effective drag force are independent constituents governing RT dynamics in each of the regimes (linear, nonlinear, and mixing). We discuss the implications of our considerations for experiments in plasmas, including inertial confinement fusion [9,34].

Fluid instabilities and interfacial mixing are a subject of active research in contemporary plasma physics [4,5,10,12,13,21,22,24,27,28,39,40]. The novelty of our work is in the rigorous theoretical solution of RT dynamics and in the detailed analysis of the numerical data of Front Tracking numerical simulations. The group theory approach precisely derives from the governing equations the independent driving constituents of RT dynamics. The interpretation of the analyzed numerical data of Front Tracking simulations agrees with experiments. Our results provide extensive benchmarks for plasma processes to which RT instability and RT mixing are relevant.

Our paper has the following structure. After Introduction, we proceed to the mathematical problem of RT dynamics in Section 2. Section 3 presents the group theory approach. Section 4 provides analytical results for the buoyancy and drag parameters in linear and nonlinear RTI. Section 5 illustrates in detail the properties of special analytical solutions in the nonlinear regime. Section 6 investigates analytically the regime of self-similar mixing. Section 7 presents the numerical approach. Section 8 studies the merger process of two-dimensional bubbles and presents the two-dimensional merger trees. Section 9 focuses on the merger process of three-dimensional bubbles and presents three-dimensional merger tree figures. Section 10 investigates simulation growth rates. Section 11 reviews the merger model. Section 12 discusses the application of the merger model to experiments. Section 13 discusses the outcomes of our work. Section 14 presents the conclusions.

2. Mathematical Problem of RT Dynamics

We consider immiscible, inviscid fluids of differing densities separated by a well-defined interface. The heavier fluid sits above the lighter fluid and the entire system is subject to a time-dependent downwards acceleration field. The dynamics of such fluids are governed by conservation of mass, momentum, and energy:

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}\left(\rho v_i\right)=0,$$

$$\frac{\partial }{\partial t}\left(\rho v_j\right)+\frac{\partial }{\partial x_i}\left(\rho v_i{v}_j\right)+\frac{\partial P}{\partial x_j}=0,$$

$$\frac{\partial E}{\partial t}+\frac{\partial }{\partial x_i}\left((E+P)v_i\right)=0,$$

(1)

where $i=1,2,3$, $(x_1,x_2,x_3)=(x,y,z)$ are the spatial coordinates, t is time, ($\rho ,\mathbf{v},P,E)$ are the fields of density $\rho$, velocity $\mathbf{v}$, pressure P, and energy density $E=\rho (e+\frac{1}{2}{\mathbf{v}}^2)$, and e is the specific internal energy [4,13,32].

It is necessary that momentum is conserved at the interface and that there can be no mass flow across it. Hence, the boundary conditions at the interface are

$$\left[\rho \left(\frac{\dot{\theta }}{|\nabla \theta |}+\mathbf{v}·\mathbf{n}\right)\right]=0,\left[\mathbf{v}·\mathbf{n}\right]=0,\left[P\right]=0,$$

$$\left[\mathbf{v}·\mathit{\tau }\right]=\mathrm{arbitrary},\left[W\right]=\mathrm{arbitrary},$$

(2)

where $[\cdots ]$ denotes the jump of functions across the interface; $W=e+P/\rho$ is the specific enthalpy; $\mathbf{n}$ and $\mathbf{\tau }$ are the normal and tangential unit vectors of the interface with $\mathbf{n}=\mathbf{\nabla }\theta /|\mathbf{\nabla }\theta |$ and $\mathbf{n}·\mathbf{\tau }=0$; $\theta =\theta (x,y,z,t)$ is a local scalar function, with $\theta =0$ at the interface and $\theta >0$ $(\theta <0)$ in the bulk of the heavy (light) fluid, indicated hereafter by subscript $h\left(l\right)$. The initial conditions consist of initial perturbations of the flow fields and the interface. The flow is periodic in plane $(x,y)$ normal to the z-direction, as is set by the initial conditions [4,32,38].

Instability is driven by acceleration. Acceleration can be due to a body force and be directed from the heavy to the light fluid. It can also be the acceleration of the (nearly planar) interface (i.e., the inertial acceleration) and be directed from the light to the heavy fluid. In our work, acceleration $\mathbf{g}$ is due to a body force and is directed from the heavy to the light fluid along the z-axis $\mathbf{g}=(0,0,-g)$. This acceleration modifies the pressure field: $P\to P-\rho gz$. We assume that there are no mass sources and hence the velocity in the upper (lower) fluid moves to zero as $z\to \infty$ ($z\to -\infty$). The Atwood number $A=({\rho }_h-{\rho }_l)/({\rho }_h+{\rho }_l)$ parameterizes the ratio of the fluid densities [4,13,32].

The mathematical problem of RTI requires one to solve the system of nonlinear partial differential equations in four-dimensional space-time, the boundary value problems at the unstable nonlinear interface and at the outside boundaries, and also the ill-posed initial value problem with account for non-locality and singularities [4,13,41,42]. In our work, we solve this problem analytically by employing group theory approach [4,13] and numerically by analyzing data from Front Tracking simulations.

3. Group Theory Approach

3.1. Large-Scale Coherent Structures

The large-scale coherent structure is an array of bubbles and spikes periodic in the plane normal to the acceleration direction. The initial conditions determine the length and time scales of the problem. The length scale is the wavelength $\lambda$ of the array, or the associated wave-vector $k^{-1}$, with $\lambda \sim k^{-1}$. As a specific example, we consider the flow which has square symmetry, in which case the length scale $\lambda$ and the wave-vector k are related as $k=2\pi /\lambda$.

The group theory approach [4,13] accurately describes both the linear and nonlinear dynamics of the unstable interface. In particular, for the nonlinear dynamics, there are families of asymptotic solutions with the number of relevant parameters defined by flow symmetry; the dynamics are multi-scale and characterized by the contributions of two macroscopic length scales, these being the spatial period and the amplitude of the interfacial coherent structure. Furthermore, in any regime, the dynamics are essentially interfacial, with intense motion of the fluids near the interface and effectively no motion away from the interface; see works [4,13].

Even for ideal incompressible fluids, with $\dot{e}=0$ and $\nabla e=0$, system Equations (1) and (2) are extremely challenging [4,13,32]. Even with this amount of complexity, the dynamics nevertheless have remarkable features of universality and order and can thus be approached using group theory [4]. Here, we implement the group theory approach and compare the principal resulting equation with the momentum model, the equations of which have the same symmetries and scaling transformations as conservation laws [4].

3.2. Analytical Methodology

RT dynamics are known to possess certain features of universality and order, and can be analyzed from first principles by applying group theory [4]. For scale-dependent dynamics, group theory can employ space groups to derive from governing equations a dynamical system and find its solutions. For scale-dependent and self-similar dynamics, group theory can be realized in a momentum model with the same symmetry transformations as governing equations [13].

For spatially periodic flows, scale-dependent RT dynamics are invariant with respect to space group G, whose generators are translations in the plane, rotations and reflections, and which also has anisotropy in the acceleration direction and inversion in the normal plane, such as the groups of square $p4mm$ or hexagon $p6mm$ in three-dimensional flows.

By deriving the structures of flow fields in the bulk from the governing equations, considering dynamics with potential velocity fields ${\mathbf{v}}_{h\left(l\right)}=\nabla {\Phi }_{h\left(l\right)}$, and applying irreducible representations of the relevant group, we can expand the flow fields as Fourier series and perform a spatial expansion in the vicinity of a regular point of the interface—the tip of a bubble or a spike with coordinates $(0,0,z_0\left(t\right))$ and velocity $(0,0,v(t\left)\right)$ with $v\left(t\right)={\dot{z}}_0\left(t\right)$ (the dot indicating a partial time derivative). This reduces the governing equations to a dynamical system for moments (each of which is a correlation function represented by an infinite sums of Fourier amplitudes) and surface variables. For square group $p4mm$, to the first order, $N=1$, the interface is $z-z_0\left(t\right)=\zeta \left(t\right)(x^2+y^2)$ and the dynamical system is

$${\rho }_h\left(\dot{\zeta }-4M_1\zeta -\frac{M_2}{2}\right)=0,{\rho }_l\left(\dot{\zeta }-4{\tilde{M}}_1\zeta +\frac{{\tilde{M}}_2}{2}\right)=0,$$

(3)

$$\begin{gathered} (1+A)\left(\frac{{\dot{M}}_1}{2}+\zeta {\dot{M}}_0-\frac{M_1^2}{2}-\zeta g\right) \\ =(1-A)\left(\frac{\dot{\widetilde{M_1}}}{2}-\zeta \dot{{\tilde{M}}_0}-\frac{{\tilde{M}}_1^2}{2}-\zeta g\right), \end{gathered}$$

(4)

$$M_0=-{\tilde{M}}_0=-v\left(t\right).$$

(5)

Here, $M\left(\tilde{M}\right)$ are the moments of the heavy (light) fluid, and $v\ge 0(\le 0),\zeta \le 0(\ge 0)$ are the velocity and curvature of the interface at the tip of the bubble (spike). The system length-scale is $k^{-1}$, where the wavevector is $k\sim {\lambda }^{-1}$ and the wavelength is $\lambda$. The time scale is $\tau =1/\sqrt{kg}$. Group theory is further applied to solve the closure problem and find linear and nonlinear solutions for scale-dependent RT dynamics. The initial curvature and velocity are $\zeta \left(t_0\right)$ and $v\left(t_0\right)$, respectively, and $t_0\gg \tau$ is the initial time.

3.3. The Momentum Model

The momentum model [4,13,33,38] has the same symmetries and scaling transformations as the conservation laws, and balances, per unit mass, the rates of gain and loss of specific momentum. In doing so, it identifies invariant, scaling, and spectral properties of the dynamics.

The dynamics of a parcel of fluid are governed by the balance per unit mass of the rate of momentum gain $\mu$ and the rate of momentum loss $\tilde{\mu }$ as

$$\frac{dh}{dt}=v,\frac{dv}{dt}=\tilde{\mu }-\mu ,$$

(6)

where h is the displacement along acceleration direction $\mathbf{g}$, v is the corresponding velocity, $\tilde{\mu }$ and $\mu$ are the magnitudes of the rates of gain and loss of specific momentum in that direction [4,13]. Rates of gain $\tilde{\mu }$ and loss $\mu$ of specific momentum are associated with rates of gain $\widetilde{ϵ}$ and loss $ϵ$ of specific energy as $\tilde{\mu }=\widetilde{ϵ}/v$ and $\mu =ϵ/v$. The rates of energy gain and dissipation are $\widetilde{ϵ}=Bgv$ and $ϵ=C{v}^3/L$, with B and C being the buoyancy and drag parameters and L being the length scale for energy dissipation. This leads to $\tilde{\mu }=Bg$ and $\mu =C{v}^2/L$, and hence

$$\frac{dh}{dt}=v,\frac{dv}{dt}=Bg-C\frac{v^2}{L},$$

(7)

with $h=h\left(t\right)>0$ and $v=v\left(t\right)>0$ for $t>t_0>0$. Length scale L for energy dissipation can be the horizontal scale (wavelength) $L\sim \lambda \sim k^{-1}$, the vertical scale (amplitude) $L\sim h$, or a combination of these. Cases $L\sim \lambda \sim k^{-1}$ and $L\sim h$ correspond to scale-dependent and scale-invariant dynamics, respectively.

The dynamical system in terms of moments and surface variables enables the study of the interplay of harmonics and asymptotic solutions, whereas the momentum model reveals explicit solutions and the asymptotic limits of RT dynamics.

3.4. Effect of Symmetry and Dimension

The group theory approach accommodates the analysis of two-dimensional and three-dimensional RT dynamics. In this work, we present analytical solutions for a three-dimensional RT flow. Two-dimensional analytical solutions can be obtained from three-dimensional solutions with proper transformations [13,38].

4. Analytical Quantities for Buoyancy and Drag

To derive expressions for buoyancy and drag parameters in RT dynamics, we directly link the dynamical system and the momentum model.

4.1. Scale-Dependent Linear Dynamics

In the dynamical system, we associate the curvature with the amplitude, $k{z}_0=-4\zeta /k$, and retain only first-order harmonics in expressions for the moments. In the momentum model, we relate the vertical scale with amplitude $h=z_0$, and the length-scale with wavevector, $L=k^{-1}$ [38]. This transforms the governing equations to

$$\dot{h}=v,\dot{v}=B_{l}g-C_{l}k{v}^2,\zeta =-\frac{k^2}{4}{z}_0,v=-M_0={\tilde{M}}_0,$$

(8)

where the buoyancy parameter is $B_l=Ak{z}_0=-4A\zeta /k$ with $B_l>0$ for bubbles and $B_l<0$ for spikes, and the drag parameter is constant, $C_l=A/2>0$; subscript l stands for linearity. By integrating these equations, we obtain solutions for RT dynamics of bubbles (spikes) moving up (down) and being concave down (up) with $z_0>0$, $v>0$, $\zeta <0$ ($z_0<0$, $v<0$, $\zeta >0$):

$$z_0\left(t\right)=C_1{e}^{\sqrt{Ak}t}+C_2{e}^{-\sqrt{Ak}t},$$

$$v\left(t\right)=\sqrt{Ak}{C}_1{e}^{\sqrt{Ak}t}-\sqrt{Ak}{C}_2{e}^{-\sqrt{Ak}t},\zeta \left(t\right)=-\frac{k^2}{4}{z}_0\left(t\right),$$

where $C_1$ and $C_2$ are determined by initial conditions $z\left(t_0\right)$ and $v\left(t_0\right)$. These results hold for $k{z}_0\ll 1$ and $v\ll \sqrt{B_{l}g/C_{l}k}$.

The formation of bubbles and spikes is defined by the initial conditions, with bubbles developing for $\zeta \left(t_0\right)<0$ and spikes developing for $\zeta \left(t_0\right)>0$. For both bubbles and spikes, the magnitude of the curvature increases with time, and the magnitude of the velocity also increases.

In the linear regime, buoyancy grows with amplitude, and the constant drag is set by the initial conditions. The linear RT dynamics of bubbles and spikes are single-scale and defined by the horizontal scale—the wavelength of the coherent structure.

4.2. Scale-Dependent Nonlinear Dynamics

As time progresses, the bubble/spike amplitude increases, its curvature approaches a constant value, and RT dynamics become nonlinear.

In a dynamical system, we impose a constant curvature condition and retain higher-order harmonics in the moments. In the momentum model, we relate the vertical scale with the amplitude, $h=z_0$, and the length-scale with the wavevector, $L=k^{-1}$. This transforms the governing equations to

$${\dot{z}}_0=v,\dot{v}=B_{n}g-C_{n}k{v}^2,\zeta =k\xi ,v=-M_0={\tilde{M}}_0.$$

(9)

For square symmetry, buoyancy parameter $B_n$ and drag parameter $C_n$ are, respectively,

$$B_n(A,\zeta )=\frac{-2A\zeta (9k^2-64{\zeta }^2)}{3k^3+10A{k}^2\zeta -128A{\zeta }^3},$$

$$C_n(A,\zeta )=\frac{k^3(9A{k}^2-48k\zeta +64A{\zeta }^2)}{(9k^2-64{\zeta }^2)(3k^3+10A{k}^2\zeta -128A{\zeta }^3)},$$

(10)

in domain $[-{\zeta }_{\mathrm{cr}},{\zeta }_{\mathrm{cr}}]$ with ${\zeta }_{\mathrm{cr}}=\frac{3}{8}k$. Sub-script n stands for nonlinearity.

Figure 1 and Figure 2 are the plots of buoyancy and drag parameters, respectively, in nonlinear RT dynamics.

For RT bubbles, the values of buoyancy and drag parameters are $B_n\ge 0$, $C_n\ge 0$ for $\zeta \in (-{\zeta }_{\mathrm{cr}},0)$ with $B_n\in [0,{\widehat{B}}^{*}]$, ${\widehat{B}}^{*}={\widehat{B}}^{*}\left(A\right)$, and $C_n\in [0,\infty )$. For RT spikes, the values are $B_n\le 0$ with $[{\check{B}}^{*},0]$, ${\check{B}}^{*}={\check{B}}^{*}\left(A\right)$, for $\zeta \in (0,{\zeta }_{\mathrm{cr}})$, and $C_n\le 0$ with $C_n\in (-\infty ,0)]$, for $\zeta \in ({\zeta }_{\mathrm{A}},{\zeta }_{\mathrm{cr}})$. Here, ${\zeta }_{\mathrm{A}}={\zeta }_{\mathrm{A}}\left(A\right)$, ${\zeta }_{\mathrm{A}}\in (0,{\zeta }_{\mathrm{cr}})$ is the curvature of the fastest spike; we call it the Atwood spike. A spike has unbounded velocity when the drag parameter is zero. This occurs for $\zeta ={\zeta }_{\mathrm{A}}$, where

$$\frac{{\zeta }_{\mathrm{A}}}{{\zeta }_{\mathrm{cr}}}=\frac{1-\sqrt{1-A^2}}{A}.$$

(11)

Hence, we find that for nonlinear RT dynamics, there is a family of buoyancy and drag values parameterized by the interface morphology, that is, the curvature of the bubble/spike; for given $(A,\xi )$, values $(B_n,C_n)$ are constant. The non-uniqueness is due to the presence of shear at the interface. This can be seen by defining the shear function $\mathsf{\Gamma }$ as the spatial derivative of the jump of tangential velocity at the interface. It is $\mathsf{\Gamma }={\tilde{M}}_1-M_1=6k^{3}v/(9k^2-64{\zeta }^2)$ near the tip of the bubble/spike.

Asymptotically, the solution of the second of Equations (9) is

$$v=\pm \sqrt{\frac{B_{n}g}{C_{n}k}}.$$

For square symmetry $p4mm$, this velocity is

$$v=\pm \sqrt{\frac{g}{k}}\frac{9k^2-64{\zeta }^2}{k^2}\sqrt{\frac{-2Ak\zeta }{64A{\zeta }^2-48k\zeta +9A{k}^2}},$$

(12)

and the corresponding shear function is

$$\mathsf{\Gamma }=\pm 6\sqrt{kg}\sqrt{\frac{-2Ak\zeta }{64A{\zeta }^2-48k\zeta +9A{k}^2}}.$$

(13)

In Equations (12) and (13), the positive sign applies for bubbles and the negative sign applies for spikes. Figure 3 and Figure 4 show the bubble tip velocity and shear function, respectively, as functions of the bubble curvature. Figure 5 and Figure 6 show the spike tip velocity and shear function, respectively, as functions of the spike curvature.

Differentiating Equation (12) with respect to $\zeta$ and setting the result to zero shows that curvature ${\zeta }_{\mathrm{vmax}}$ which maximizes velocity satisfies equation

$${\zeta }_{\mathrm{vmax}}^4-\frac{k}{A}{\zeta }_{\mathrm{vmax}}^3+\frac{9k^2}{32}{\zeta }_{\mathrm{vmax}}^2-k^4{\left(\frac{3}{16}\right)}^3=0.$$

(14)

By solving this equation and substituting the solution into Equation (12), we find that the corresponding maximum velocity is

$$v_{\max }=16\sqrt{2}{\left|\frac{{\zeta }_{\mathrm{vmax}}}{k}\right|}^{\frac{3}{2}}.$$

(15)

This is a universal relation between the curvature and the maximum velocity. We refer to this bubble as the ‘Atwood’ bubble to emphasize its dependence on the Atwood number.

RT dynamics are set by the interplay of buoyancy and drag for a given interface morphology and shear. The dynamics of RT bubbles are regular and influenced by a few competing factors: bubbles with larger curvature have larger buoyancy and move faster than flattened bubbles; yet, bubbles with larger curvature have larger shear and larger drag reducing their velocities; for curved bubbles, the shear alone can maintain the pressure at the interface, leading to zero buoyancy and infinite drag.

The dynamics of RT spikes are singular. While the magnitude of the spike’s buoyancy is qualitatively similar to that of the bubble, the spike’s drag vanishes $C_n\to 0$ for $\zeta \to {\zeta }_{\mathrm{A}}$. This leads to a singularity indicating that the RT Atwood spike has velocity and shear growing quickly with $\left|v\right|/\sqrt{g/k}\to \infty$ and $\left|\mathsf{\Gamma }\right|/\sqrt{kg}\to \infty$.

The nonlinear RT dynamics of bubbles/spikes are multi-scale and defined by the two macroscopic length scales—the vertical and horizontal scales—the amplitude and wavelength of the structure, respectively. Mathematically, the (as we call it) multi-scalability is associated with the dependence of the buoyancy and the drag parameters on the interface morphology and shear. Physically, it can be understood by viewing RT coherent structures as standing waves with growing amplitudes. This reveals the mechanisms of the transition of RT dynamics from the scale-dependent to scale-invariant regimes. Particularly, the traditional merger and amplitude dominance mechanisms can both lead to the acceleration of RT bubbles/spikes due to drag force reduction, and can transfer scale-dependent dynamics $L\sim k^{-1}\sim \lambda$ to scale-invariant mixing $L\sim \left|h\right|$.

5. Special Solutions

We note the presence of special solutions in the nonlinear family. Some of these special solutions include the Taylor bubble/spike with curvature magnitude $k/8$, similar to works [3,38]; the Layzer-drag bubble/spike with velocity dependence $\sqrt{2A/(1\pm A)}$, similar to works [37,43]; the fastest Atwood bubble with $({\widehat{\zeta }}_{\mathrm{A}},{\widehat{v}}_{\mathrm{A}},{\widehat{\mathsf{\Gamma }}}_{\mathrm{A}})$ and the fastest Atwood spike with $({\check{\zeta }}_{\mathrm{A}},{\check{v}}_{\mathrm{A}},{\check{\mathsf{\Gamma }}}_{\mathrm{A}})$; see works [13,38].

5.1. Special Solutions for Nonlinear Bubbles

We introduce nondimensional variables

$$\widehat{\zeta }=-\frac{\zeta }{k},\widehat{v}=\frac{v}{\sqrt{g/k}},\widehat{\mathsf{\Gamma }}=\frac{\mathsf{\Gamma }}{\sqrt{kg}},\widehat{B}=B,\widehat{C}=C.$$

5.1.1. The Taylor Bubble

We refer to this bubble as a ‘Taylor bubble’ since its curvature is the same as in work [3] of Davis and Taylor except for a difference in the wavevector value [4]. For the Taylor bubble, the curvature, velocity, shear function, and buoyancy and drag parameters are

$${\widehat{\zeta }}_{\mathrm{T}}=\frac{1}{8},{\widehat{v}}_{\mathrm{T}}=\sqrt{\frac{8A}{3+5A}},{\widehat{\mathsf{\Gamma }}}_{\mathrm{T}}=\frac{3}{4}{\widehat{v}}_{\mathrm{T}},$$

$${\widehat{B}}_{\mathrm{T}}=\frac{2A}{3-A},{\widehat{C}}_{\mathrm{T}}=\frac{3+5A}{4(3-A)}.$$

(16)

5.1.2. The Layzer-Drag Bubble

The family of bubbles has a solution with velocity $\widehat{v}=\sqrt{2A/(1+A)}$. We call this solution the ‘Layzer-drag bubble’ since the drag model applies this velocity rescaling with the Atwood number to the single-mode Layzer first-order approximation at $A=1$ [37,43]. Experiments and simulations tend to compare well with this rescaling [37,43,44]. For the Layzer-drag bubble, the curvature, velocity, and shear function are

$${\widehat{\zeta }}_{\mathrm{L}}={\widehat{\zeta }}_{\mathrm{L}}\left(A\right),{\widehat{v}}_{\mathrm{L}}=\sqrt{\frac{2A}{1+A}},{\widehat{\mathsf{\Gamma }}}_{\mathrm{L}}=\frac{6{\widehat{v}}_{\mathrm{L}}}{9-64{\widehat{\zeta }}_{\mathrm{L}}^2},$$

(17)

and the buoyancy and drag parameters are as defined in Equation (10). The dependence of the Layzer-drag bubble curvature on the Atwood number is cumbersome. For fluids with very different densities, $A\to 1$, in the solution for the Layzer-drag bubble,

$${\widehat{\zeta }}_{\mathrm{L}}\to \frac{1}{8}\left(1+\sqrt{\frac{1-A}{6}}\right),$$

$${\widehat{v}}_{\mathrm{L}}\to 1-\frac{1-A}{4},{\widehat{\mathsf{\Gamma }}}_{\mathrm{L}}\to \frac{3}{4}+\frac{\sqrt{6(1-A)}}{32},$$

$${\widehat{B}}_{\mathrm{L}}\to 1+\sqrt{\frac{1-A}{6}},{\widehat{C}}_{\mathrm{L}}\to 1+\sqrt{\frac{1-A}{6}}.$$

For fluids with very similar densities, $A\to 0$, in the solution for the Layzer-drag bubble,

$${\widehat{\zeta }}_{\mathrm{L}}\to \frac{\sqrt{9-4\sqrt{3}}}{8},{\widehat{v}}_{\mathrm{L}}\to \sqrt{2A},{\widehat{\mathsf{\Gamma }}}_{\mathrm{L}}\to \sqrt{\frac{3A}{2}},$$

$${\widehat{B}}_{\mathrm{L}}\to \frac{8\sqrt{3}{\widehat{\zeta }}_{\mathrm{L}}}{3}A,{\widehat{C}}_{\mathrm{L}}\to \frac{4\sqrt{3}}{3}{\widehat{\zeta }}_{\mathrm{L}}.$$

The Layzer-drag bubble is more curved for fluids with similar densities, $A\to 0$, when compared to fluids with very different densities, $A\to 1$. For fluids with very different densities, $A\to 1$, the values of the curvature, velocity, shear, and group-momentum coefficients of the Layzer-drag bubble are the same as the values of the corresponding quantities of the Taylor bubble. For $0<A<1$, the Layzer-drag bubble is more curved and has smaller velocity and larger shear when compared to the Taylor bubble.

5.1.3. The Atwood Bubble

We refer to the fastest member of the family of solutions as the ‘Atwood bubble’ to emphasize its dependence on the Atwood number [38]. The curvature, velocity, shear function, and buoyancy and drag parameters of the Atwood bubble are

$${\widehat{\zeta }}_{\mathrm{A}}={\widehat{\zeta }}_{\mathrm{vmax}},{\widehat{v}}_{\mathrm{A}}=16\sqrt{2{\widehat{\zeta }}_{\mathrm{A}}^3},{\widehat{\mathsf{\Gamma }}}_{\mathrm{A}}=\frac{6{\widehat{v}}_{\mathrm{A}}}{9-64{\widehat{\zeta }}_{\mathrm{A}}^2},$$

$${\widehat{B}}_{\mathrm{A}}=\frac{8192{\widehat{\zeta }}_{\mathrm{A}}^4}{9-320{\widehat{\zeta }}_{\mathrm{A}}^2-8192{\widehat{\zeta }}_{\mathrm{A}}^4},{\widehat{C}}_{\mathrm{A}}=\frac{16{\widehat{\zeta }}_{\mathrm{A}}}{9-320{\widehat{\zeta }}_{\mathrm{A}}^2-8192{\widehat{\zeta }}_{\mathrm{A}}^4},$$

(18)

where ${\zeta }_{\mathrm{vmax}}=-k{\widehat{\zeta }}_{\mathrm{vmax}}$ is given by Equation (14). We note that

$${\widehat{B}}_{\mathrm{A}}={\widehat{C}}_{\mathrm{A}}{\widehat{v}}_{\mathrm{A}}^2.$$

For fluids with very different densities, $A\to 1$, in the solution for the Atwood bubble,

$${\widehat{\zeta }}_{\mathrm{A}}\to \frac{1}{8}\left(1-\frac{1-A}{8}\right)$$

$${\widehat{v}}_{\mathrm{A}}\to 1-\frac{3}{16}(1-A),{\widehat{\mathsf{\Gamma }}}_{\mathrm{A}}\to \frac{3}{4}-\frac{21}{128}(1-A),$$

$${\widehat{B}}_{\mathrm{A}}\to 1-\frac{13}{8}(1-A),{\widehat{C}}_{\mathrm{A}}\to 1-\frac{5}{4}(1-A).$$

For fluids with very similar densities, $A\to 0$, in the solution for the Atwood bubble,

$${\widehat{\zeta }}_{\mathrm{A}}\to \frac{3}{16}{A}^{\frac{1}{3}},{\widehat{v}}_{\mathrm{A}}\to \sqrt{\frac{27A}{8}},{\widehat{\mathsf{\Gamma }}}_{\mathrm{A}}\to \sqrt{\frac{3A}{2}},$$

$${\widehat{B}}_{\mathrm{A}}\to \frac{9}{8}{A}^{\frac{4}{3}},{\widehat{C}}_{\mathrm{A}}\to \frac{1}{3}{A}^{\frac{1}{3}}.$$

For fluids with very different densities, $A\to 1$, the values of the curvature, velocity, shear function, and group-momentum coefficients of the Atwood bubble are very close to the values of the corresponding quantities of the Taylor and the Layzer-drag bubbles. For fluids with a finite density ratio, the curvatures relate as ${\widehat{\zeta }}_{\mathrm{A}}/{\widehat{\zeta }}_{\mathrm{cr}}\in (0,1/8)$, ${\widehat{\zeta }}_{\mathrm{L}}/{\widehat{\zeta }}_{\mathrm{cr}}\in (0,1)$ and ${\widehat{\zeta }}_{\mathrm{A}}/{\widehat{\zeta }}_{\mathrm{T}}\in (0,1)$. For $0<A<1$, the Atwood bubble is less curved and has larger velocity and smaller shear when compared to the Layzer-drag and the Taylor bubble.

Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11 are plots of the nondimensional variables as functions of the Atwood number A for these special bubbles.

5.2. Special Solutions for Nonlinear Spikes

We introduce nondimensional variables

$$\check{\zeta }=\frac{\zeta }{k},\check{v}=-\frac{v}{\sqrt{g/k}},\check{\mathsf{\Gamma }}=-\frac{\mathsf{\Gamma }}{\sqrt{kg}},\check{B}=-B,\check{C}=-C.$$

5.2.1. The Taylor Spike

We refer to this spike as a ‘Taylor spike’ since its curvature magnitude is the same as in work [3] of Davis and Taylor except for a difference in the wavevector value [4]. For the Taylor bubble, the curvature, velocity, shear function, and buoyancy and drag parameters are

$${\check{\zeta }}_{\mathrm{T}}=\frac{1}{8},{\check{v}}_{\mathrm{T}}=\sqrt{\frac{8A}{3-5A}},{\check{\mathsf{\Gamma }}}_{\mathrm{T}}=\frac{3}{4}{\check{v}}_{\mathrm{T}},$$

$${\check{B}}_{\mathrm{T}}=\frac{2A}{3+A},{\check{C}}_{\mathrm{T}}=\frac{3-5A}{4(3+A)}.$$

(19)

For the Taylor spike, the velocity and the shear function have non-positive values, which are zero for $A\to 0^{+}$, finite for $0<A<3/5$, and which become singular approaching negative infinity for $A\to 3/5$.

5.2.2. The Layzer-Drag Spike

The family of spikes has amongst them a solution with velocity $v=-\sqrt{2A/(1-A)}$. We call this solution the ‘Layzer-drag spike’ since the drag model applies this velocity rescaling with the Atwood number to the single-mode Layzer first-order approximation at $A=1$ [37,43]. Experiments and simulations tend to compare well with this rescaling [37,43,44]. For the Layzer-drag bubble, the curvature, velocity, and shear function are

$${\check{\zeta }}_{\mathrm{L}}={\check{\zeta }}_{\mathrm{L}}\left(A\right),{\check{v}}_{\mathrm{L}}=\sqrt{\frac{2A}{1-A}},{\check{\mathsf{\Gamma }}}_{\mathrm{L}}=\frac{6{\check{v}}_{\mathrm{L}}}{9-64{\check{\zeta }}_{\mathrm{L}}^2},$$

(20)

and the buoyancy and drag parameters are as defined in Equation (10). The dependence of the Layzer-drag spike curvature on the Atwood number is cumbersome. For fluids with very different densities, $A\to 1$, in the solution for the Layzer-drag spike,

$${\check{\zeta }}_{\mathrm{L}}\to \frac{3}{8}\left(1-\sqrt{2(1-A)}\right),$$

$${\check{v}}_{\mathrm{L}}\to \sqrt{\frac{2}{1-A}},{\check{\mathsf{\Gamma }}}_{\mathrm{L}}\to \frac{1}{6}\sqrt{\frac{2}{1-A}},$$

$${\check{B}}_{\mathrm{L}}\to \frac{9}{11}\left(1-\frac{4}{11}\sqrt{2(1-A)}\right),{\check{C}}_{\mathrm{L}}\to \frac{9}{22}(1-A).$$

For fluids with very similar densities, $A\to 0$, in the solution for the Layzer-drag spike,

$${\check{\zeta }}_{\mathrm{L}}\to \frac{\sqrt{9-4\sqrt{3}}}{8},{\check{v}}_{\mathrm{L}}\to \sqrt{2A},{\check{\mathsf{\Gamma }}}_{\mathrm{L}}\to \sqrt{\frac{3A}{2}},$$

$${\check{B}}_{\mathrm{L}}\to \frac{8\sqrt{3}{\check{\zeta }}_{\mathrm{L}}}{3}A,{\check{C}}_{\mathrm{L}}\to \frac{4\sqrt{3}{\check{\zeta }}_{\mathrm{L}}}{3}.$$

The Layzer-drag spike is more curved for fluids with very different densities $A\to 1$ than for fluids with very similar densities $A\to 0$. For a two-fluid system with $A\in (0,1)$, ${\check{\zeta }}_{\mathrm{L}}/{\check{\zeta }}_{\mathrm{T}}\in (\sqrt{9-4\sqrt{3}},3)$, and the Layzer-drag spike is more curved than the Taylor spike. The Layzer-drag spike has smaller absolute velocity and smaller absolute shear than the Taylor spike.

For fluids with very similar densities $A\to 0$, the magnitude of the curvature of the Layzer-drag spike is the same as that for the Layzer-drag bubble, with ${\check{\zeta }}_{\mathrm{L}}={\widehat{\zeta }}_{\mathrm{L}}$. The magnitudes of the velocity and shear of the Layzer-drag spike differ from those of the Layzer-drag bubble: the velocity and shear magnitudes of the spike approach infinity for $A\to 1$, whereas the velocity and the shear magnitudes of the bubble are finite for any A.

5.2.3. The Atwood Spike

The fastest spike is where the drag parameter is zero. Its curvature is that given in Equation (11). We refer to this spike as the ‘Atwood spike’ to emphasize its dependence on the Atwood number [13,38]. It is the fastest growing, and hence the physically significant spike. For the Atwood spike, the curvature is regular, $\check{\zeta }\to {\check{\zeta }}_{\mathrm{A}}$, whereas the velocity and the shear function are singular, with ${\check{v}}_{\mathrm{A}}/\sqrt{g/k}\to \infty$, for any Atwood number.

For fluids with very different densities, $A\to 1$, the solution for the Atwood spike is

$${\check{\zeta }}_{\mathrm{A}}=\frac{3}{8}\left(1-\sqrt{2(1-A)}\right),B_{\mathrm{A}}=\frac{9}{11}\left(1-\frac{4}{11}\sqrt{2(1-A)}\right).$$

For fluids with very similar densities, $A\to 0$, the solution for the Atwood spike is

$${\check{\zeta }}_{\mathrm{A}}=\frac{3}{16}A,B_{\mathrm{A}}=\frac{9}{8}{A}^2.$$

The Atwood spike is less curved than the Layzer-drag spike for $A\in [0,1)$ and has the same curvature as the Layzer-drag spike ${\check{\zeta }}_{\mathrm{A}}/{\check{\zeta }}_{\mathrm{L}}=1$ at $A=1$. The Atwood spike attains the same curvature as the Taylor spike ${\check{\zeta }}_{\mathrm{A}}/{\check{\zeta }}_{\mathrm{T}}=1$ at $A=3/5$.

The properties of the curvature, velocity, and shear of the Atwood spike differ dramatically from those of the Atwood bubble: for any density ratio A, the magnitudes of the velocity and shear of the Atwood spike are singular, whereas the magnitudes of the velocity and shear of the Atwood bubble are regular and finite. While the magnitudes of the curvature of the Atwood spike and the Atwood bubble are both finite, the Atwood spike is more curved (less curved) when compared to the Atwood bubble for fluids with very different (very similar) densities $A\to 1(A\to 0)$.

Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 are plots of nondimensional variables as functions of the Atwood number A for these special spikes.

6. Self-Similar Mixing

In the mixing regime, the dynamical system is generalized as

$$\dot{v}=B_{m}g-C_m\frac{v^2}{h},\dot{h}=v.$$

(21)

In the momentum model, we relate the length-scale for energy dissipation with the amplitude as $L=\left|h\right|=|z_0|$, and we set $B=B_m$, $C=C_m$, where subscript m stands for mixing. The buoyancy and drag are free values due to the many scales contributing, with $B_m,C_m>0(<0)$ for bubbles (spikes) with $h>0,v>0$ ($h<0,v<0)$. Herein, we consider $B_m>0,C_m>0,h>0,v>0$. Our consideration can also be applied for $B_m<0,C_m<0,h<0,v<0$; see works [4,33,36].

We solve the model equation in domain $h,v>0$, $t>{\tilde{t}}_0$ for constant $B_m,C_m>0$ and with initial conditions ${\tilde{h}}_0=h\left({\tilde{t}}_0\right)$, ${\tilde{v}}_0=v\left({\tilde{t}}_0\right)$:

$$\ddot{h}+C_m\frac{{\dot{h}}^2}{h}-B_{m}g=0.$$

The solution is sought in the physical domain with $h,v,t>0$ and with initial conditions $t=t_0$, $h_0=h\left(t_0\right)$, and $v_0=\dot{h}\left(t_0\right)$ [4,36].

The associated homogeneous equation is $\ddot{h}=-C{\dot{h}}^2/h$. It can be solved by the method of separation of variables leading to the general solution, with integration constants $H_0$, $V_0$, as

$$\frac{h_d^{C+1}-H_0^{C+1}}{H_0^C}=(C+1)V_0(t-t_0),$$

where subscript d stands for dissipative dynamics. For $t/t_0\gg 1$, the solution is

$$h_d=B_d{t}^{\frac{1}{1+C}},B_d=H_0{\left((C+1)\frac{V_0}{H_0}\right)}^{\frac{1}{C+1}}.$$

To find a particular solution for the non-homogeneous equation $\ddot{h}+C{\dot{h}}^2/h=g$, we invoke Lie group analysis and employ one-parameter scaling transformation $\widehat{t}=e^{\psi }t$, $\widehat{h}=e^{\gamma \psi }h$, where $\psi$ and $\gamma$ are the transformation’s parameter and constant. This leads to $e^{(\gamma -2)\psi }(\ddot{h}+C{\dot{h}}^2/h)=g$ and hence $\gamma =2$ and the corresponding invariant $h/t^2$. The particular solution is $h_a=B_{a}g{t}^2/(2+4C)$, where subscript a stands for accelerative dynamics. The general solution is a function of the dissipative and accelerative solutions, $h=h(h_d,h_a)$.

Whilst solutions $h_d$ and $h_a$ should be coupled because of the nonlinear nature of the differential equation, it has, in fact, the remarkable property that in the respective asymptotic regimes, $(t-t_0)/t)\ll 1$ and $t/t_0\gg 1$, solutions $h_a$ and $h_d$ are effectively decoupled. This decoupling is due to the distinct symmetries of the solutions, $h_a$ and $h_d$ [4,36].

For $L\sim \lambda$, the gains and losses are fully balanced, $\mu =\tilde{\mu }$ and $ϵ=\widetilde{ϵ}$. For $L\sim h$, the gains and losses are unbalanced, $\mu \ne \tilde{\mu }$ and $ϵ\ne \widetilde{ϵ}$ [4,36]. According to data, the imbalance is small, $(\tilde{\mu }-\mu )/\tilde{\mu }\ll 1$ [5]. This slight imbalance results in the flow acceleration and in the development of self-similar RT mixing: parcel of the heavy fluid moves in the acceleration direction, parcel of the light fluid moves against the acceleration direction, the position of the center of mass of the fluid system changes, energy is gained, almost all of this gain effectively dissipates, and its remnant leads to flow acceleration [4,36]. The transition from nonlinear RTI to RT mixing may occur via the growth of the horizontal scale, as the merger model predicts [5].

7. Numerical Approach

Assuming that acceleration is uniform in time, we recall that the penetration of the bubble interface into the heavier fluid occurs on an acceleration-determined time scale and follows a $t^2$ scaling law. The penetration distance of the light fluid into the heavy fluid is known to be given by

$$h_b={\alpha }_{b}Ag{t}^2,$$

(22)

where $A=({\rho }_2-{\rho }_1)/({\rho }_2+{\rho }_1)$ is the Atwood number, a relative density difference that modifies the uniform acceleration, g. Here, ${\rho }_1$ (${\rho }_2$) denotes the density of the light (heavy) fluid. The constant ${\alpha }_b$ is a growth rate parameter.

This remainder of this paper analyzes RT mixing data in depth, with archived data that could be used in future studies of bubble merger. Here, we analyze 2D experimental and 3D simulation studies of RT mixing; 2D and 3D merger trees are constructed as a result of this process. In these figures, the bubble merger process occurring over multiple time steps can be observed. Detailed analysis of all measured properties of each bubble are recorded in a spreadsheet (not included here due to its size), and will aid in any further analyses; see Section 8 and Section 9.

The mean bubble radius, $\overline{r}$, and the separation height between adjacent bubbles undergoing a merger, $h_m$, grow on a quadratic time scale [45], in common with ${\alpha }_b$. These growth rates are related by the model, with ${\alpha }_r$ and ${\alpha }_{h_m}$ being their respective growth rate parameters. We determine these three growth rates for the numerically simulated experiment. The extraction and analysis of these data, including relevant plots, are presented in Section 10.

The bubble merger model correctly predicts the measured merger rates relative to experimental data. The model has one adjustable parameter, which is a noise level (fluctuation or variance) in the bubble radii. A brief overview of the merger model and all relevant equations is given in Section 11.

On a quantitative level, we revisit predictions made by the merger model and reported in [45]. We find that for the Smeeton and Youngs (SY) experiments [46], the merger model better predicts growth rate ${\alpha }_b$ than originally reported. Applications of this model to numerical simulation data are also presented; see Section 12. We comment that ${\alpha }_b$ is sensitive to the amplitude of long wavelength random (noise) perturbations to the initial data, as was observed in the experiments of Meusche, Andrews and Schilling [47,48,49] and the LEM experiments of Dimonte [50] whose long-wavelength noise in the experiment apparatus made its use for comparison to simulations difficult. The bubble merger model studied here excludes such initial perturbations, which appear to be missing in the SY experimental data [46] and are missing in the simulation data studied here.

We recall that the fluids are accelerated; this is customarily achieved by a rapid downward acceleration of a gravitationally stable two-fluid system, with the light fluid above the heavy. These turbulent events accelerating the bubble penetration occur at the base of the bubble, or even above the base, in the form of a vortex ring wrapping around the stem above the downward moving bubble. According to this picture, we see the origin of the randomness as lying in the fluids at the base of the bubble, or even above it, and coming from turbulent aspects of the mixing flow, a well-known source of randomness. In fact, in our data, we observe that the merger occurs in stages, starting at the base and progressing downwards to the tip. Not explored here, but likely, is the fact that the merger begins earlier, at the level of vortex rings driving the bubble.

We mention an alternate bubble merger model [37] which predicts one observed parameter (${\alpha }_b$) on the basis of one input parameter: a noise level in the initial data that is assumed to originate in the initial perturbations of the initial interface. This model mispredicts the bubble radius and velocity for non-interacting bubbles. The misprediction is significant as its velocity contains a term proportional to its radius. The bubble merger model proposed by Cheng et al. [45] yields a range of predictions for the bubble height-to-width ratio that is consistent with the values observed in experiment and simulation. However, as noted in [51], the model [52,53] predicts a value for this ratio that falls outside the observed range.

8. 2D Bubble Merger

Five experiments of SY are analyzed (numbered 99, 103, 104, 105, and 114). Plates A-F in each experiment reflect six subsequent times in the evolution of the mixing front. The elapsed time between plates varies by experiment and are shown in

Table 1. Elapsed time between plates in SY experiments [46]. Experiments differed in the densities of the two fluids involved, the tank pressure, the interfacial tension, and/or the mean acceleration. Plates A–F in each experiment reflect six subsequent times in the evolution of the mixing front.

	Time Elapsed (ms)
Experiment No.	A→B	B→C	C→D	D→E	E→F
99	28.5	10.6	15.9	16.0	10.6
103	16.1	10.7	16.1	12.4	11.5
104	16.3	10.8	11.5	11.0	11.0
105	20.6	10.3	10.2	10.3	20.4
114	10.3	10.3	10.3	13.1	10.3

. These experiments are presented as photographs taken from face on cameras, so the resulting data are essentially two-dimensional. The acceleration is tangential to the y-axis in these experiments. The data show a few of the bubbles surviving until the final plate, while many smaller bubbles, present initially, are left behind and disappear from the mixing front and from the plots as time evolves. As the fronts are moving downward in time, the later bubbles are below the earlier ones. We show the larger of the observed bubbles at each time frame. The edges of the same bubble, observed at successive times, are joined to provide a picture of the continuous evolution of each major bubble.

Analysis is omitted for Plates A and B across all experiments due to the low resolution and small amplitudes of bubbles. Plate C is also omitted for Experiments 99 and 103 for identical reasons. Plate photographs are scaled using tank dimensions from SY and center-line values from [54].

To minimize the possibility of interpretation errors in assigning heights and widths to each bubble, both values are independently measured for all bubbles in the ensemble by two authors. Leading bubbles are defined as bubbles with the furthest penetration and with sufficiently large radii. Two authors independently identify the leading bubbles of each plate.

We construct a merger path for a qualitative analysis of the bubble merger using bubble positions and widths. Construction of merger paths begins with the terminal leading bubbles of each experiment and continues working backwards in time to identify its predecessor(s).

Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21 present bubble merger trees, which contain the leading bubbles of each plate along with their adjacent neighbors. For all bubbles in the ensemble, including those not involved in a merger, we represent their widths by a horizontal line. Coloring and line styles are used to identify a bubble’s plate of origin. One can observe the merger process, where a given bubble grows systematically by expanding into space vacated by smaller bubbles. We quantify this observation in

Table 2. We quantify the bubble merger process in Rayleigh–Taylor dynamics by counting the number of bubbles lost between successive time steps, denoted as Plates A-F, in SY experiments [46]. Analysis was omitted for Plates A and B across all experiments due to the low resolution and small amplitudes of bubbles. (a) 2D experiments. The time evolution of bubble loss across various physical experiments of SY [46]. Details regarding the elapsed time between plates are given in Table 1. Plate C in Experiment 99 and Experiment 103 could not be analyzed due to insufficient bubble amplitudes and resolution. (b) 3D simulation. The time evolution of bubble loss in 3D simulation data modeling SY Experiment 112 [46]. The times elapsed between Plates C, D, E, and F in this simulation are all 10 ms.

(a)
	Number of Bubbles Lost
Experiment No.	C→D	D→E	E→F
99	N/A	3	1
103	N/A	2	3
104	6	7	2
105	11	5	2
114	12	8	0
(b)
	Number of Bubbles Lost
Experiment No.	C→D	D→E	E→F
112 (sim.)	5	2	2

a.

9. 3D Bubble Merger

We construct the three-dimensional analogue of the merger trees found in Section 8. The 3D nature of the data introduces additional difficulties in its analysis and presentation. We examine simulation data modeling the miscible and compressible SY [46], Experiment 112. This experiment involves sodium iodide solution (${\rho }_h=1.89$ g/cm${}^3$) and a solution consisting of water, hydrochloric acid, and phenolphthalein (${\rho }_l=1.00$ g/cm${}^3$), with an approximate Atwood number of 0.3. The mean acceleration is $39g_0$, where $g_0=9.81$ m/s${}^2$. The domain in this simulation has dimensions 3.75 cm × 0.625 cm × 15 cm, with periodic boundaries in the x-y plane, and it has no edge effects. The spatial mesh has spacing $\Delta x=\Delta y=\Delta z=5/192$, and the numbers of grid cells are $N_x=144$, $N_y=24$, and $N_z=567$. The initial conditions contain short-wavelength perturbations and the equation of state for this simulation is a gamma law gas. We thank Tulin Kaman for providing the simulation data used in this study. For details on a related simulation study of Experiment 112, see [55].

Simulations are analyzed at four time steps (40, 50, 60, and 70 ms), each corresponding to one of the Plates C–F from Experiment 112 of SY [46]. Data from these simulations consist of a meshed interface modeling the mixing front and discrete state variables, namely density. Unlike the 2D case of the preceding section, we note greater ambiguity in identifying bubbles of light fluid in the mixing front. We define a bubble as a structure of reasonable shape that grows from a distinct tip. An example is shown in Figure 22. We observe that in this example, the bubble structure contains a filament of heavy fluid within its enclosure, which we ignore in subsequent analyses.

At a fixed time step, we identify distinct bubble tips and record each bubble’s height at three distinct levels. The lowest height is the bubble tip, above this is a larger cross-section at the middle of the downward moving bubble, and above that is a rim that forms the upper limit of the bubble within the fluid. A bubble’s terminal rim is defined as the largest of the bubble’s heights at which its cross-section has a consistent and maximum shape. The data analysis consists of examination of cross-sections at different height (z) levels through the simulation data.

Following data collection and bubble identification, we construct a merger path for each terminal bubble. That is, for each terminal bubble, we determine what bubble(s) must have preceded it in earlier time steps. The resulting five merger paths are presented in

Table 3. Three-dimensional merger paths for successive time steps (denoted as Plates C–F) in numerical RT mixing data simulating Experiment 112 of SY [46]. The domain of this simulation is periodic in the x-y plane and has dimensions 3.75 cm × 0.625 cm × 15.0 cm. The acceleration is tangential to the z-axis. Starting with Plate C, we identify each bubble by recording the $(x,y,z)$ coordinates of its approximate bubble center and follow its position as time evolves until Plate F. This defines a merger path (MP). Here, p is used to denote splitting by the periodic boundary. MP 1–4 display the expected dynamics of a single bubble, expanding in time into space vacated by smaller bubbles. MP 5 originates from two distinct tips but terminates at a common base, and so the merger occurs at heights beyond the tips. See Figure 23, Figure 24, Figure 25, Figure 26 and Figure 27 for visualization.

MP	C	D	E	F
1.	$(3.75p,0.50,7.05)$	$(3.75p,0.625p,6.43)$	$(3.75p,0.625p,5.85)$	$(3.75p,0.625p,5.44)$
2.	$(1.40,0.625p,7.26)$	$(1.40,0.625p,6.69)$	$(1.25,0.625p,5.96)$	$(1.25,0.625p,5.31)$
3.	$(1.75,0.18,6.82)$	$(1.70,0.20,6.38)$	$(1.75,0.20,5.98)$	$(1.90,0.20,5.46)$
4.	$(2.30,0.35,7.26)$	$(2.40,0.40,6.52)$	$(2.45,0.55p,5.67)$	$(2.10,0.625p,4.71)$
	$(3.50,0.30,7.16)$
5.	$(3.40,0.625p,7.21)$	$(3.75p,0.47p,7.00)$	$(3.75p,0.53p,6.64)$	$(3.75p,0.625p,6.22)$
	$(3.25,0.15,7.21)$

; here, p is used to denote a bubble split periodically in that respective coordinate. We note in particular Bubble 5 of Table 3, which originates from two distinct tips but terminates at a common base, and so the merger must have occurred at a height beyond the tip.

We now outline the process of constructing the 3D merger trees. For a given bubble, we first extract three 2D subsets from the 3D simulation data, one for each of the heights of interest previously discussed. Each slice is digitized, converting all density values to binary values, with each grid block taking a value of zero if it contains density values closer to that of the lighter fluid, and one for values closer to the heavier fluid. A marching squares algorithm, a two-dimensional analogue of the marching cubes algorithm, is then used to identify the vertices of each bubble’s contour. We filter out the non-leading bubbles that do not have the potential to participate as a leading bubble at any time step, using bubble centers to identify the leading bubbles and then take their convex hull. The remaining non-leading bubbles are analyzed separately. Bubbles split by the periodic boundary are reconstructed using appropriate coordinate shifts. Contours are then plotted to create 3D merger trees. Each terminal bubble, and thus each merger tree, is plotted separately. Each bubble is shown at three heights, from the tip to the middle to its base. To suggest time evolution, the bases at the different times are joined by a few nearly vertical line segments. Time evolution is coded by color and line style, consistent with the 2D figures shown in Section 8, with red corresponding to Plate C, green for Plate D, blue for Plate E, and black for Plate F. See Figure 23, Figure 24, Figure 25, Figure 26 and Figure 27, with two alternate viewing angles in most cases.

We recall the picture behind the bubble merger model, which states that the larger bubbles expand by occupying space vacated by smaller bubbles that were swept aside. This process can be observed in Figure 28, Figure 29 and Figure 30, where all bubbles in the ensemble, including smaller (non-leading) bubbles, are shown. As in Figure 23, Figure 24, Figure 25, Figure 26 and Figure 27, the time evolution of each bubble is shown using nearly vertical lines connecting three of its cross-sections: the bubble tip, an intermediate height, and the bubble rim. The leading bubbles can be seen expanding with time, and non-leading bubbles are swept aside and ultimately vanish from the mixing front. We quantify these phenomena, both for the 2D and 3D cases, in Table 2b.

Another example of the bubble merger process can be found in Figure 31. This example follows a fixed region of the mixing front as time evolves. The images in these figures are taken directly from the numerical simulation’s meshed interface data and visualized using the VisIt visualization software. In particular, we point out two smaller, less advanced bubbles, which are indicated by a solid red and dashed blue arrows, respectively. The first bubble survives until Plate D, whereas the second survives until Plate E. By Plate F, we observe that both bubbles have vanished, while the nearby leading bubbles have grown in size and have penetrated further into the heavy fluid.

10. Data Collection and Simulation Growth Rates

We further analyze numerical simulation modeling Experiment 112 discussed in Section 9. In particular, we wish to determine measured values for three growth rate parameters: ${\alpha }_b$, ${\alpha }_r$, and ${\alpha }_{h_m}$. Of most importance is the bubble growth rate, ${\alpha }_b$, which describes the asymptotic growth behavior of the mixing layer resulting from RT instability [51]. We determine these parameters first by measurement, and later predict these growth rates using the merger model (Section 11 and Section 12). Values for the measured parameters presented in this section can be found in Table 8.

Required quantities are measured using interface data, as opposed to state variables, due to their superior resolution. An example of these interface data is pictured in Figure 31 using the VisIt visualization software.

10.1. Bubble Growth Rate, ${\alpha }_b$

We recall from Equation (22) that the penetration distance of the bubble front follows a $t^2$ scaling law. The leading bubble of each time step is defined as the bubble which penetrates furthest into the heavy fluid. The penetration of the leading bubble is taken as the overall penetration distance of the bubble ensemble, and is what we define as $h_b$. Penetration distances for the numerical simulation are provided by Tulin Kaman and are presented in

Table 4. Penetration distance of the leading bubble, $h_b$, of each time step for numerical simulation data modeling Experiment 112 of SY [46].

Time (ms)	${\mathit{h}}_{\mathit{b}}$ (mm)
0	0.009
2.87	0.109
11.5	0.398
25.8	0.543
45.9	0.617
71.8	0.806
103	1.51
141	3.17
184	5.64
232	8.82
287	12.4
347	16.3
413	21.7
485	26.3
563	30.8

.

The values for $h_b$ are then regressed against scaled acceleration distances $X=Ag{t}^2$ to determine ${\alpha }_b$; see Figure 32. We observe ${\alpha }_b=0.056$ in the numerical experiment. The nonzero intercept in Figure 32 suggests a delay in mixing, which is consistent with the original SY experiment.

SY report ${\alpha }_b=0.052$ for Experiment 112 [46]. We compare our measured value to SY, but do not expect and do not achieve exact agreement due to the 3.5 mm initial thickness present in the SY experiments, which is in contrast to the sharp interface of the numerical simulation. The role of this initial thickness parameter was studied in [55], and the discrepancy of the simulation ${\alpha }_b$ to the SY value is consistent with study [55].

10.2. Radius Growth Rate, ${\alpha }_r$

It was shown in [45] that the mean bubble radius, $\overline{r}$, satisfies a self-similar scaling law similar to Equation (22):

$$\overline{r}\approx {\alpha }_{r}Ag{t}^2.$$

(23)

We now determine the radius growth rate, ${\alpha }_r$, for numerical simulation data. We note that the transition from 2D to 3D introduces additional difficulties; in particular, the base of the bubble is not perfectly circular, but rather (approximately) elliptical. Thus, to have a well-defined notion of a bubble’s radius, we define the radius at the base to be half the length of its major axis. We collect radius values for all bubbles in the ensemble, including smaller non-leading bubbles. The mean radius of each plate is presented in

Table 5. Mean radius, $\overline{r}$, of each plate for numerical simulation data modeling Experiment 112. of SY [46].

Plate	Time (ms)	$\overline{\mathit{r}}$ (mm)
C	40	1.308
D	50	1.764
E	60	2.125
F	70	2.870

.

By regressing the values of Table 5 against the scaled acceleration distances, we obtain ${\alpha }_r=0.004$; see Figure 33. The results from Table 5 are consistent with the qualitative picture behind the bubble merger, where bubbles are expanding as time evolves.

10.3. Separation Height Growth Rate, ${\alpha }_{h_m}$

Another parameter within the merger model is the maximum separation height between adjacent bubbles undergoing a merger, denoted as $h_m$. According to [45], this parameter is also governed by a $t^2$ scaling law:

$$h_m={\alpha }_{h_m}Ag{t}^2.$$

(24)

This can be seen in Figure 31, where the height between adjacent bubbles increases as time evolves.

We wish to determine the growth rate parameter, ${\alpha }_{h_m}$, for numerical simulation data. To this end, we record the maximum separation distance of adjacent bubbles for each plate. The results are presented in

Table 6. Maximum separation height between adjacent bubbles, $h_m$, of each plate for numerical simulation data modeling Experiment 112. of SY [46].

Plate	Time (ms)	${\mathit{h}}_{\mathit{m}}$ (mm)
C	40	4.729
D	50	7.355
E	60	12.326
F	70	16.663

. As before, the growth rate, ${\alpha }_{h_m}$, is determined through regression; see Figure 34. We measure ${\alpha }_{h_m}=0.0323$.

11. Review of Merger Model

Here, we present a brief overview of the model and model equations; we refer the reader to [45] for details. The merger model predicts three measured parameters:

1.

the mixing rate, ${\alpha }_b$,

2.

the bubble height separation at the time of merger, $h_m$,

3.

the mean bubble radius, $\overline{r}$,

and

takes in a noise level as an adjustable input.

This noise level, as presented in the original merger model, is taken to be the variance of the fluctuations amongst bubble radii. More explicitly, for fixed time, the randomness within the bubble ensemble can be measured by first applying map $r\mapsto (r/\overline{r}-1$), where r is the bubble radius, and computing the variance of the transformed distribution. Area preservation and geometrical arguments concerning bubble merger lead to defining parameter

$$k\equiv \left(r_{\mathrm{new}}-\frac{r_1+r_2}{2}\right)/\left(\frac{r_1+r_2}{2}\right),$$

(25)

where

$$\begin{array}{cc}\hfill r_1& =\overline{r}(1+\sigma ),\hfill \\ \hfill r_2& =\overline{r}(1-\sigma ),\hfill \\ \hfill r_{\mathrm{new}}& =\sqrt{r_1^2+r_2^2},\hfill \end{array}$$

and $\sigma$ is the standard deviation of fluctuations amongst bubble radii. After some elementary manipulations, Equation (25) can be rewritten to be independent of $\overline{r}$:

$$k=\sqrt{2(1+{\sigma }^2)}-1.$$

(26)

We recall that $h_m$ and $\overline{r}$ satisfy self-similar scaling laws:

$$\overline{r}\left(t\right)\approx {\alpha }_{r}Ag{t}^2,{\alpha }_r=\frac{k}{2}{〈\frac{1}{t_m^{\prime }}〉}_{*}^2,$$

(27)

and

$$h_m={\alpha }_{h_m}Ag{t}^2,{\alpha }_{h_m}=\frac{k^2}{4}{h}_m^{\prime }{〈\frac{1}{t_m^{\prime }}〉}_{*}^2=h_m^{\prime }{\alpha }_r,$$

(28)

where $h_m^{\prime }=h_m/\overline{r}$, and ${\langle 1/t_m^{\prime }\rangle }_{*}$ is the averaged merger rate evaluated at the renormalization group fixed point.

Using the above, the model predicts ${\alpha }_b$ as

$$\begin{array}{c}\hfill {\alpha }_b=\frac{k}{4}\left(c_b+\frac{1+k}{2}{h}_m^{\prime }{〈\frac{1}{t_m^{\prime }}〉}_{*}\right){〈\frac{1}{t_m^{\prime }}〉}_{*}=\frac{1}{2}{c}_b{\alpha }_r^{1/2}+\left(\frac{1}{2k}+\frac{1}{2}\right){\alpha }_{h_m},\end{array}$$

(29)

where $c_b$ is a constant obtained from experiment or from theoretical determination. In all subsequent analyses, we follow [45] and take $c_b=0.532$ for hexagonal bubbles.

Thus, given a measure of randomness, a prediction for ${\alpha }_b$ using Equation (29) reduces to determining $h_m^{\prime }$ and ${\langle 1/t_m^{\prime }\rangle }_{*}$. To this end, we consider two interacting adjacent bubbles with radii $r_1$ and $r_2$. We let $r_l=min(r_1,r_2)$ and $r_h=max(r_1,r_2)$. We note that ratio $r_h/r_l$ is then itself a measure of randomness. When using experimental data, this value can be computed directly from the standard deviation of bubble radii using

$$\frac{r_h}{r_l}=\frac{\overline{r}(1+\sigma )}{\overline{r}(1-\sigma )}=\frac{1+\sigma }{1-\sigma }.$$

(30)

The velocity of a single bubble, when used in conjunction with geometric considerations of interacting bubble tips on a larger length scale, provides a criterion for completion of merger:

$$\begin{array}{c}\hfill \frac{6}{{\chi }^{1/2}}=\frac{5}{{\lambda }^{1/2}}{\left(\frac{e^{\lambda x_m/2}-\lambda x_m/2-1}{e^{\lambda x_m/2}+1}\right)}^{1/2}+{\left(\frac{e^{x_m}-x_m-1}{e^{x_m}+1}\right)}^{1/2},\end{array}$$

(31)

where

$$\begin{array}{cc}\hfill R& \equiv r_h+r_l\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill x_m& \equiv h_m{\beta }_1/R\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \chi & \equiv 1+r_h/r_l\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \lambda & \equiv \frac{\chi }{1+\sqrt{\left(1+{(\chi -1)}^2\right)/2}},\hfill \end{array}$$

(35)

and ${\beta }_1$ is the first zero of the Bessel function, $J_1\left(r\right)$.

We observe that solving Equation (31) for $x_m$ and using Equation (33), along with the fact that $R/\overline{r}=(r_h+r_l)/[(r_h+r_l)/2]=2$, yields $h_m^{\prime }$, as desired.

The model determines ${\langle 1/t_m^{\prime }\rangle }_{*}$ by using the already found value for $x_m$:

$${〈\frac{1}{t_m^{\prime }}〉}_{*}=\frac{1}{x_m^2}{\int }_0^{x_m}\frac{x_m-x}{t_m^{\prime }\left(x\right)}dx,$$

(36)

where

$$\begin{array}{c}\hfill t_m^{\prime }\equiv 3\sqrt{\frac{2}{{\beta }_1}}{\int }_x^{x_m}{\left[\sqrt{\frac{e^x-x-1}{e^x+1}}+\frac{5}{\sqrt{\lambda }}\sqrt{\frac{e^{\lambda x/2}-\lambda x/2-1}{e^{\lambda x/2}+1}}\right]}^{-1}dx.\end{array}$$

(37)

Having determined the averaged merger rate at the RNG fixed point, we may use Equations (27)–(29) to predict theoretically the growth rates, ${\alpha }_r$, ${\alpha }_{h_m}$, and ${\alpha }_b$.

12. Application of Model to Experiment

We use the process outlined in Section 11 to predict model parameters. We begin first by applying the model to various radius ratios $r_h/r_l$ and reproduce Table III of [45] keeping all notation consistent; the results are shown in

Table 7. Merger model predictions with various radius ratios $r_h/r_l$ given as an input. This table is a reproduction of Table III of [45], but with a minor discrepancy in the numerical solution for Equation (31) corrected.

${\mathit{r}}_{\mathit{h}}/{\mathit{r}}_{\mathit{l}}$	k	$\mathit{\chi }$	${\mathit{x}}_{\mathit{m}}$	${\mathit{h}}_{\mathit{m}}^{\prime }$	${\langle 1/{\mathit{t}}_{\mathit{m}}^{\prime }\rangle }_{*}$	${\mathit{\alpha }}_{\mathit{r}}$	${\mathit{\alpha }}_{\mathit{b}}$	${\mathit{h}}_{\mathit{b}}/\left(2\overline{\mathit{r}}\right)$
1.0	0.414	2.0	3.488	1.82	0.432	0.00799	0.049	3.04
1.1	0.416	2.1	3.325	1.74	0.440	0.00839	0.049	2.93
1.2	0.420	2.2	3.181	1.66	0.448	0.00886	0.050	2.82
1.3	0.427	2.3	3.051	1.59	0.455	0.00945	0.051	2.70
1.4	0.434	2.4	2.934	1.53	0.462	0.01005	0.052	2.59
1.5	0.442	2.5	2.829	1.48	0.468	0.01071	0.053	2.49

. We then apply the model to Experiments 104, 105, and 114 of SY, and compare the results with those reported in [45]. We present these predictions in

Table 8. Merger model predictions for various SY [46] experiments, one of which is simulated numerically.

Experiment	$\mathit{\sigma }$	${\mathit{r}}_{\mathit{h}}/{\mathit{r}}_{\mathit{l}}$	k	${\langle 1/{\mathit{t}}_{\mathit{m}}^{\prime }\rangle }_{*}$	${\mathit{\alpha }}_{\mathit{r}}^{\mathbf{exp}}$	${\mathit{\alpha }}_{\mathit{r}}$	${\mathit{\alpha }}_{{\mathit{h}}_{\mathit{m}}}^{\mathbf{exp}}$	${\mathit{\alpha }}_{{\mathit{h}}_{\mathit{m}}}$	${\mathit{\alpha }}_{\mathit{b}}^{\mathbf{exp}}$	${\mathit{\alpha }}_{\mathit{b}}$
104	0.405	2.36	0.526	0.5065	0.0177	0.0087	0.0205	0.0278	0.065	0.068
105	0.39	2.28	0.518	0.5037	0.0170	0.0117	0.0200	0.021	0.064	0.072
114	0.31	1.90	0.481	0.4887	0.0138	0.0108	0.0179	0.0359	0.059	0.060
112 (sim.)	0.43	2.48	0.534	0.5104	0.0004	0.0188	0.0323	0.0211	0.056	0.067

. For consistency, we use the measured standard deviation reported in [45] as an input into the model.

Minor discrepancies in the numerical solution for nonlinear Equation (31) are corrected, as well as the modifications to the experimental predictions. ${\alpha }_b$ prediction improves in comparison to the experiment, while predictions for ${\alpha }_r$ and ${\alpha }_{h_m}$ become somewhat worse. These changes are made with the acceptance of the authors of [45].

Additionally, we apply the model to the 3D numerical simulation modeling Experiment 112 of SY that is discussed in Section 9 and Section 10. We measure the radius at the base of each bubble in the ensemble and measure the fluctuations within the resulting distribution as was discussed in Section 11, similar to the 2D analysis carried out in [45]. Due to the 3D nature of the simulation, the cross-sections at the base are not perfectly circular and so we note some ambiguity in interpreting its radius. We map these standard deviations to radius ratios using Equation (30) as input to the model. The results can be found in the last row of Table 8.

We recall from Section 10.1 that we cannot confirm that the sharp interface in the numerical simulations is relevant to SY’s Experiment 112 [46]. Furthermore, the unavailability of data from earlier time steps renders our application of the merger model to simulation incomplete. We posit that this exclusion has an influence on the predicted value for ${\alpha }_b$. For these reasons, we do not compare the merger model prediction to the ${\alpha }_b$ value measured in either experiment or simulation.

13. Outcomes

We studied RTI driven by constant acceleration in neutral plasmas (fluids and/or gases). We employed the group theory analytical approach and the analysis of numerical data in Front Tracking simulations. We revealed the independent driving constituents that RT flow possesses in each of the linear, nonlinear, and mixing regimes.

13.1. Outcomes from Analysis

In our analysis, we employed the group theory approach and directly linked the governing equations to the dynamical system and the momentum model [4,38]. We precisely derived the buoyancy and drag parameters for RT bubbles and RT spikes in the linear, nonlinear, and mixing regimes. The governing equations include the conservation laws in the bulk, the boundary value problem at the interface and the initial value problem; see Equations (1), (2) and (6). They are reduced to the dynamical system by using an irreducible representation of the group of RT flow. The momentum model has the same symmetries and scaling transformations as the governing equations; see Equations (3)–(7).

In each regime, we precisely derived from the governing equations the momentum model parameters for RT bubbles and RT spikes. We integrated the momentum model equations and investigated their asymptotic limits, finding complete consistency with asymptotic analytical solutions for the dynamical system; see Equations (8) and (9) and Figure 3, Figure 4, Figure 5 and Figure 6. In addition to the interplay of harmonics, the group theory approach captures the interfacial character of RT dynamics, with intense fluid motion at the interface, effectively no motion far from the interface, and with shear-driven vortical structure at the interface; see works [13,38]. Our approach can be applied to three-dimensional and two-dimensional RT flows [4].

We revealed that the buoyancy and drag parameters are independent characteristics of RT dynamics. The expressions for the buoyancy and drag parameters are distinct in the early-time linear, later-time nonlinear, and mixing regimes. They are also distinct for RT bubbles and RT spikes. In addition to the density ratio, they depend on the interface morphology and the interfacial shear. However, they are independent of the acceleration, thus expanding the validity of our results to variable, in space and time, accelerations; see Equations (8)–(10) and (21) and Figure 1, Figure 2 and Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16.

13.2. Outcomes from Simulations

In regard to simulations, we presented 2D and 3D merger tree diagrams showing the dynamics of the merger model. In these figures, we observed that smaller and less advanced bubbles are pushed downstream and away from the mixing front. Larger and leading bubbles expand into this vacated space. Hence, we concluded that the qualitative picture behind the bubble merger is correct. Figures showing the growth of a single bubble as a result of this process are also presented, including an incomplete merger with two distinct tips merging to a common base; see Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21 and Figure 23, Figure 24, Figure 25, Figure 26, Figure 27, Figure 28, Figure 29, Figure 30 and Figure 31.

We analyzed numerical simulation data modeling Experiment 112 of SY and presented a detailed data analysis; see Section 10. Due to the significant difference in the initial thickness of the mixing layer in the SY experiment (3.5 mm) and the numerical simulation (0 mm), we did not compare our parameter values to the SY data.

On the quantitative side, we revisited predictions made by the merger model [45] and corrected the numerical solution to a nonlinear equation pertaining to the model. The result is an improvement in ${\alpha }_b$ for all three experiments, and some worsening of the two other growth rate parameters ${\alpha }_r$ and ${\alpha }_{h_m}$. We applied this model to the numerical simulation data and presented similar findings; we predicted ${\alpha }_b=0.067$. We did not compare this prediction to the value measured in experiment or simulation due to the absence of early time data; see Table 7 and Table 8.

The success of this model indicates that the direct measurement of randomness related to the bubble front evolution, in terms of a measured variance and standard deviation, is feasible both for simulation data and data collected from physical (2D) experiments. The results of this analysis may motivate future studies of bubble merger to further investigate the growth of randomness in hopes of constructing a fully theoretical model.

13.3. Outcomes for Plasma Experiments

Our analysis and simulations suggests the following implications in our considerations for plasma experiments, including inertial confinement fusion [9,34,36]. (1) RT dynamics are driven by two independent constituents—the buoyancy and the drag. The values of the buoyancy and drag parameters depend on the regime and on the deterministic conditions, and are distinct for bubbles and spikes; see Equations (8), (10) and (21) and Figure 1, Figure 2 and Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16. (2) In order to accurately describe experimental and numerical simulation data for uncontrolled broadband initial conditions, one needs to consider two independent self-similar processes—the merger process, as introduced in Section 1, and the amplitude growth process; see Equations (8), (9) and (21) and Figure 1 and Figure 2. (3) RT dynamics can, in principle, be controlled, by means of the deterministic conditions (i.e., the initial and the flow experimental conditions), in agreement with experiments in fluids and plasmas at high Reynolds numbers [13,38].

Our results suggest new opportunities for control of plasma flows by means of these deterministic conditions [9,33,34,36]. For instance, in inertial confinement fusion (ICF), it may be worth scratching the target in order to impose proper deterministic conditions and to gain better control of fluid instabilities and the interfacial mixing. This may help existing methods focused on fine polishing of ICF targets in order to fully eliminate RTI [9]. Our results can also be applied in nanofabrication, by describing the dynamics and morphology of the RT unstable interface in the scale-dependent linear and nonlinear regimes, and their scaling properties, along with the sensitivity to deterministic conditions in the self-similar mixing regime [26].

14. Conclusions

We studied, analytically and numerically, RTI with constant acceleration in neutral plasma (fluid) by employing the group theory approach and the merger model. Based on the obtained analytical solutions and conclusions drawn from the data analysis of numerical and physical experiments, we revealed the independent constituents governing RT dynamics in the linear, nonlinear, and mixing regimes. Our results can be used to enhance methods of control of RT unstable plasma flows.