2. Different behaviours between spikes and bubbles

It is well known that spikes and bubbles can exhibit qualitatively different behaviours. Spikes are more unstable than bubbles and consequently it is more difficult to predict the behaviour of spikes. Even among fingers of the same type, different fluid density ratios can lead to quantitatively different behaviours. The larger the density ratio, the more unstable is the finger. In addition, the spike in a system with an infinite density ratio shows completely different behaviour at late times from those in systems with finite density ratios. For RM instability, the former attains a constant velocity asymptotically while the latter decay with time; for RT instability, the former asymptotically grows with a constant acceleration while the latter grow with approximately constant velocities at the quasi-steady stage. These phenomena are well illustrated in figures 2–4 for RM instability, and in figure 5 for RT instability. In the captions of these figures, $a_{0}$ is the initial amplitude of the interface perturbation, $v_{0}$ is the initial velocity, $k$ is the wavenumber, ${\it\lambda}$ is the wavelength and $g$ is the gravitational acceleration. The data for figures 2 and 5 are taken from Sohn (Reference Sohn2004), in which the numerical simulations were performed with the vortex method (E & Hou Reference E and Hou1990). The data for figure 3 are taken from Dimonte & Ramaprabhu (Reference Dimonte and Ramaprabhu2010), and the data for figure 4 are taken from Alon et al. (Reference Alon, Hecht, Ofer and Shvarts1995). To describe density ratios, the Atwood number $A=({\it\rho}_{1}-{\it\rho}_{2})/({\it\rho}_{1}+{\it\rho}_{2})$ is commonly used, where ${\it\rho}_{i}~(i=1,2)$ denotes the density of fluid $i$ . Conventionally only a positive Atwood number is used for both spikes and bubbles of the same system, but in this paper we specifically use the positive Atwood number for bubbles, and its negative counterpart for spikes with the same density ratio. This allows us to express the results for both spikes and bubbles by the same function in the later analysis.

Figure 2.

Spike (solid curves) and bubble (dashed curves) velocities at different Atwood numbers $A$ from numerical simulations for RM instability (Sohn Reference Sohn2004). The initial conditions are $a_{0}=0.5$  cm and $v_{0}=0.8$  cm ms $^{-1}$ , with $A=0,\pm 0.3,\pm 0.7,\pm 1$ .

Figure 3.

Spike (solid curves) and bubble (dashed curves) velocities at different Atwood numbers $A$ from numerical simulations for RM instability (Dimonte & Ramaprabhu Reference Dimonte and Ramaprabhu2010). The initial condition is $a_{0}k=0.125$ , with $A=\pm 0.25,\pm 0.5,\pm 0.75,\pm 0.88,\pm 0.94,\pm 0.98$ .

Figure 4.

Spike (solid curves) and bubble (dashed curves) velocities at different Atwood numbers $A$ from numerical simulations for RM instability (Alon et al. Reference Alon, Hecht, Ofer and Shvarts1995). The initial conditions are $v_{0}=1$  cm ms $^{-1}$ and ${\it\lambda}=1$  cm, with $A=\pm 0.2,\pm 0.5,\pm 0.9,-1$ .

Figure 5.

Spike (solid curves) and bubble (dashed curves) velocities at different Atwood numbers $A$ from numerical simulations for RT instability. Blue curves: vortex model results from Sohn (Reference Sohn2004), with $a_{0}k=0.5$ , $v_{0}=0$ and $A=\pm 0.05,\pm 0.3,\pm 0.7,\pm 1$ . Black curve (labelled as c.m.): conformal mapping result from Menikoff & Zemach (Reference Menikoff and Zemach1983), with the same initial conditions and $A=1$ .

In this study, we investigate both spikes and bubbles in systems with all density ratios up through the quasi-steady stage, and report a very surprising result: by appropriately scaling the physical quantities, the main behaviours of growth rates of all fingers collapse onto a single curve. We call this curve ‘universal’ because it applies (1) among bubbles of different density ratios, (2) among spikes of different density ratios and (3) between bubbles and spikes of different density ratios.

The Layzer-type method has been extensively adapted in studying fingering instabilities. This method was initially introduced by Layzer (Reference Layzer1955) to study the motion of a vacuum bubble in RT instability. The model was later extended by Hecht et al. (Reference Hecht, Alon and Shvarts1994) to study bubbles in RM instability and by Mikaelian (Reference Mikaelian1998) to study bubbles in cylindrical geometry. Zhang (Reference Zhang1998) extended the model to the spike. All these studies are for systems with infinite density ratios. Goncharov (Reference Goncharov2002) generalized the method to finite density ratios. However, Mikaelian (Reference Mikaelian2008) pointed out the problems in Goncharov’s model: ‘The failure for spikes becomes perhaps obvious for this model $\ldots$  it is well known that spikes are ‘sharper,’ i.e. have larger $|{\it\eta}_{2}|$ than bubbles, especially at large $A$ . (Goncharov’s) equation (4) is independent of $A$ and clearly predicts the wrong curvature for Goncharov’s spikes.’ Later Abarzhi et al. (Reference Abarzhi, Glimm and Lin2003) and Sohn (Reference Sohn2003) obtained similar growth rates with different approaches, however only for bubbles. Mikaelian (Reference Mikaelian2003) proposed an analytic model for the evolution of the bubble amplitude from the linear to the nonlinear regime. From these works, we observe a list of characteristics that the finger curvature possesses:

1. (1) for infinite density ratio, the bubble curvature asymptotically tends to one-sixth of the wavenumber (Hecht et al. Reference Hecht, Alon and Shvarts1994; Mikaelian Reference Mikaelian1998; Zhang Reference Zhang1998; Goncharov Reference Goncharov2002; Abarzhi et al. Reference Abarzhi, Glimm and Lin2003; Sohn Reference Sohn2003);

2. (2) for infinite density ratio, the spike curvature asymptotically tends to infinity (Zhang Reference Zhang1998);

3. (3) the bubble curvature in the quasi-steady stage is insensitive to the density ratio (Sohn Reference Sohn2004).

Based on these properties we will deduce the finger growth rate for a system with any density ratio and show the universality of all fingers at late times, but before the possible re-acceleration stage.

3. Theoretical formulation

We consider incompressible, inviscid and irrotational fluids with arbitrary density ratio in two dimensions, whose governing equations (Layzer Reference Layzer1955) are

Here $z={\it\eta}(t,x)$ is the vertical position of the interface at time $t$ , $g$ is the gravitational acceleration, ${\it\phi}_{i}$ is the velocity potential, ${\it\rho}_{i}$ is the fluid density, subscript $i~(i=1,2)$ denotes the upper fluid and lower fluid, respectively, and $f(t)$ is an arbitrary function that depends only on time $t$ .

In our configuration, a signed Atwood number is used to describe the density ratio. To elaborate, when studying bubbles, we consider the heavy fluid is on top of the light fluid with the acceleration pointing downwards. Therefore, the Atwood number $A=({\it\rho}_{1}-{\it\rho}_{2})/({\it\rho}_{1}+{\it\rho}_{2})>0$ and $g>0$ for bubbles. When studying spikes, we consider the light fluid is on top of the heavy fluid with the acceleration pointing upwards. Therefore, $A<0$ and $g<0$ for spikes. This setting allows us to consider only the case of the finger penetrating from the lower fluid into the upper fluid for both finger types (bubbles and spikes).

The shape of a finger tip (bubble or spike) is approximated as a parabola, ${\it\eta}(t,x)=z_{0}(t)+{\it\xi}(t)kx^{2}$ . Here $z_{0}(t)$ and ${\it\xi}(t)$ are the vertical position and the curvature of the finger tip, respectively. In our configuration, the curvature is always non-positive, i.e. ${\it\xi}\leqslant 0$ .

In theory, the velocity potential in phase 1 can be expressed in a series such as ${\it\phi}_{1}=\sum _{j=1}^{\infty }a_{1}^{(j)}(t)\cos (jkx)\text{e}^{-jkz}=a_{1}^{(1)}(t)\cos (kx)\text{e}^{-kz}+\sum _{j=2}^{\infty }a_{1}^{(j)}(t)\cos (jkx)\text{e}^{-jkz}$ . This usual standard perturbation expansion with modes $k,2k,3k,\ldots$ is good for early-time behaviour, but not suitable for late-time approximation, since it is a divergent series for large times. Since our study focuses on the late-time behaviour of single-mode RT and RM instabilities, an approach that captures the dominant behaviour of fingers at late times is needed. The works in Layzer (Reference Layzer1955), Hecht et al. (Reference Hecht, Alon and Shvarts1994), Mikaelian (Reference Mikaelian1998) and Zhang (Reference Zhang1998) showed that, for $A=1$ , modelling the system with only the primary mode $k$ is sufficient to obtain a good approximation. However, as pointed out by Mikaelian (Reference Mikaelian2008), for a system with $|A|\neq 1$ , still modelling the system with only one mode had led to incorrect predictions for the asymptotic curvature. Thus, it is necessary to model the collective behaviour of all modes beyond the primary mode $k$ . In our model, the first term $a_{1}^{(1)}(t)\cos (kx)\text{e}^{-kz}$ , abbreviated as $a_{1}(t)\cos (kx)\text{e}^{-kz}$ , represents the dominant behaviour of the finger, and we approximate the collective behaviour of all remaining terms $\sum _{j=2}^{\infty }a_{1}^{(j)}(t)\cos (jkx)\text{e}^{-jkz}$ by $b_{1}(t)\cos (c(t)kx)\text{e}^{-c(t)kz}$ . Note that in our approximation the coefficient $c(t)$ depends on the Atwood number $A$ , which is different from other approximations. Most approximations assume that the finger can be well expressed by the first few terms of the series $\sum _{j=1}^{n}a^{(j)}(t)\cos (jkx)\text{e}^{-jkz}$ , and each frequency is proportional to $k$ and independent of $A$ . Since we use an effective mode to represent the collective behaviour of all remaining modes $j=2,3,\ldots ,$ and the coefficients of these modes depend on $A$ , it follows that $c(t)$ in our approximation must also depend on $A$ . We comment that $c(t)k$ is not an approximation for only the first few terms in the expansion, and it is not for the whole interface. It is a local approximation for the collective behaviour of all remaining modes near the finger tip. Although the same functional form $c(t,A,g)$ is used for both bubbles and spikes, in our signed Atwood number notation, namely, $A>0,~g>0$ for bubbles and $A<0,~g<0$ for spikes, we have $c_{bb}=c(t,|A|,|g|)$ and $c_{sp}=c(t,-|A|,-|g|)$ . Therefore, the value of $c_{bb}$ is different from that of $c_{sp}$ in a system with a given density ratio. The velocity potential in phase 2 is approximated in a similar way. Therefore, we write the velocity potentials in the following form:

We comment that $a_{i}(t)$ , $b_{i}(t)$ , $c(t)$ , $z_{0}(t)$ and ${\it\xi}(t)$ are also functions of $A$ and $g$ . For conciseness, we do not display them explicitly.

By substituting expression (3.4) into (3.2) and expanding the resulting equations through the order of $x^{2}$ , we obtain:

One can solve for $a_{i}$ and $b_{i}$ explicitly from (3.5) and (3.6). After substituting the resulting expressions into (3.3), one can obtain an ordinary differential equation for the term proportional to $x^{2}$ :

Here $v=\text{d}z_{0}/\text{d}t$ is the velocity of the finger tip, and

The results from numerical simulations (Tryggvason Reference Tryggvason1988; Sohn Reference Sohn2004) and from theoretical analysis (Alon et al. Reference Alon, Hecht, Ofer and Shvarts1995; Goncharov Reference Goncharov2002) showed that, in the quasi-steady state, spikes have approximately constant velocities in RT instability with finite density ratios. Only in an infinite density ratio system, i.e. $A=-1$ , does the spike approach free fall asymptotically. Therefore, we will consider the fingers with $A\in (-1,1]$ in our analysis. If a theoretical model is consistent, one would expect that its solution should have the property that, as $A\rightarrow -1$ , it approaches the behaviour of spikes in an infinite density ratio system. As we will see later, our solution indeed has this important consistency property.

We now examine the velocity of the finger in the quasi-steady stage. It has been shown that both the velocity and the curvature of the finger tip are insensitive to time in this stage (Glimm et al. Reference Glimm, Li and Lin2002; Sohn Reference Sohn2004), namely, $\text{d}v/\text{d}t=0$ , $\text{d}{\it\xi}/\text{d}t=0$ and $\text{d}c/\text{d}t=0$ . From this property, we obtain an equation from (3.7) for RT instability in the quasi-steady stage:

where the constants $v_{qs}$ and ${\it\xi}_{qs}$ stand for the quasi-steady velocity and the quasi-steady curvature, respectively. The right-hand side of (3.14) must be non-negative for all values of $A$ . In this equation, $A$ and $g$ share the same sign, and ${\it\xi}_{qs}$ is always non-positive. Therefore, to maintain a non-negative $v_{qs}^{2}$ it is required that

which is equivalent to

where $q=(c_{qs}+1)/{\it\xi}_{qs}$ . In the case of spikes, i.e. when $A<0$ , the solution of (3.16) is

In particular, (3.17) must hold in the neighbourhood of $A=-1$ , which gives

namely, $c_{qs}/{\it\xi}_{qs}=-6-1/{\it\xi}_{qs}+o(1)=-6+o(1)$ , since $1/{\it\xi}_{qs}=o(1)$ by property (2) at the end of § 2. Therefore, we obtain $c_{qs}=-6\,{\it\xi}_{qs}(1+o(1))$ in the neighbourhood of $A=-1$ . To the leading-order approximation, we will approximate $c_{qs}$ by $-6{\it\xi}_{qs}$ . We now determine the behaviour of ${\it\xi}_{qs}$ near $A=-1$ . After approximating $c_{qs}$ by $-6{\it\xi}_{qs}$ , a condition for ${\it\xi}_{qs}$ can be deduced from (3.15):

where ${\it\zeta}={\it\xi}_{qs}\,\sqrt{1+A}$ . In the neighbourhood of $A=-1$ , the left-hand side of (3.19) will be determined by the leading term:

The equality in (3.19) holds for the case in which the expression inside the square bracket in the denominator of (3.14) tends to zero. This corresponds to the largest growth rate for the spike near $A=-1$ . Therefore, to obtain the most unstable growth rate, we require $72{\it\zeta}^{2}=1$ , i.e.

Equation (3.21) suggests that ${\it\xi}_{qs}(A)$ should have a Laurent series in terms of $\sqrt{1+A}$ , and we will keep only the three leading terms, namely,

Property (1) of § 2 requires, when $A=1$ , that ${\it\xi}_{qs}(1)=-1/6$ , which means

To capture property (3), we take $(\text{d}{\it\xi}_{qs}/\text{d}A)|_{A=1}=0$ , which provides another equation:

Equations (3.23) and (3.24) give $d_{0}=0$ and $d_{1}=-1/(12\sqrt{2})$ . Thus, we obtain the quasi-steady curvature

from (3.22), and the quasi-steady velocity

from (3.14) for RT instability. We comment that the result given by (3.26) is valid for both spikes ( $A<0,~g<0$ ) and bubbles ( $A>0,~g>0$ ).

As given by (3.26) $v_{qs}$ is the velocity at the quasi-steady stage. The leading-order contribution at early times is just the initial condition $v_{0}$ . It is desirable to have an asymptotic matched solution for both early and quasi-steady stages. To achieve this, we approximate ${\it\xi}(t)$ by ${\it\xi}_{qs}$ and $c(t)$ by $-6{\it\xi}_{qs}$ . This is motivated by the analytical results for $A=1$ (Mikaelian Reference Mikaelian1998; Zhang Reference Zhang1998). Equations (4) and (7) in Zhang (Reference Zhang1998) show that ${\it\xi}(t)$ approaches its asymptotic value much faster than $v(t)$ does: ${\it\xi}(t)$ approaches its asymptotic value at the rate of $\text{e}^{-3kz_{0}}$ ; $v(t)$ approaches its asymptotic value at the rate of $\text{e}^{-[3kz_{0}-\text{ln}(kz_{0})]/2}$ , which is much slower than $\text{e}^{-3kz_{0}}$ . Therefore, it is reasonable to approximate ${\it\xi}(t)$ as a constant. Since, in our analysis, the asymptotic value of $c(t)$ is related to that of ${\it\xi}(t)$ , it is reasonable to approximate $c(t)$ as a constant as well. Under these approximations, (3.7) becomes

where

which is a function of the Atwood number $A$ . Therefore, a matched solution for $v$ can be obtained by integrating (3.27) from 0 to $t$ .

4. Universality of fingers at all density ratios

For RM instability, we take the limit $g\rightarrow 0$ in (3.27) and obtain

and the matched solution is

where $v_{0}$ denotes the initial velocity. When $A=1$ , the expression given by (4.2) recovers the results shown by Mikaelian (Reference Mikaelian1998) and Buttler et al. (Reference Buttler, Oró, Preston, Mikaelian, Cherne, Hixson, Mariam, Morris, Stone, Terrones and Tupa2012).

If we introduce the scaled dimensionless variables

then (4.2) can be written as

Figure 6.

Comparison between universality curve given by (4.4) and scaled data from figure 2 for RM instability at different Atwood numbers. The solid curves are the scaled spike velocities, and the dashed curves are the scaled bubble velocities. Red curve: universality model prediction. Blue curves: vortex model results from Sohn (Reference Sohn2004) with $a_{0}=0.5$  cm  $v_{0}=0.8$  cm ms $^{-1}$ and $A=0,\pm 0.3,\pm 0.7,1$ . Black curve and green curve: $a_{0}=0$ and $v_{0}=0.5$  cm ms $^{-1}$ , with $A=1$ . The black curve (labelled as c.m.) is the conformal mapping result from Menikoff & Zemach (Reference Menikoff and Zemach1983), and the green curve is the vortex model result from Sohn (Reference Sohn2004).

It is completely unexpected that, despite the different behaviours among fingers with various Atwood numbers, and even the qualitative distinctions between spikes and bubbles, the growth rates of all finger tips, up through the quasi-steady stage, can be approximately described by a single equation in terms of the appropriately scaled dimensionless variables given by (4.3). This equation is independent of finger types and of the Atwood number. This means that, in terms of scaled variables, the dominant behaviours of the growth rates of all fingers (both spikes and bubbles) at all Atwood numbers follow a universal curve for systems governed by (3.1)–(3.3). Furthermore, the result of the bubble growth rate at any specific Atwood number can be used to predict (through scaling) not only the bubble growth rates at other Atwood numbers, but also the spike growth rates at all Atwood numbers. Similarly, one can also use the result of a spike to do such predictions. This implies that, to study the growth rates for both finger types and for all Atwood numbers, one only needs to study the behaviour of one finger at any Atwood number. This is particularly important and useful for numerical simulations and experiments, since it is highly time-consuming and also possibly expensive to perform studies for multiple Atwood numbers up through the nonlinear stage.

To verify our prediction on universality, we apply (4.3) to scale all the results from numerical simulations shown in figures 2–4, and compare the scaled results with our prediction given by (4.4). The comparison results are presented in figures 6–8. All results from three sets of data show that the scaled data are in good agreement with our universality curve. The scaled results also show that, in terms of scaled time, longer simulations are needed for spikes with large Atwood numbers to reach the quasi-steady stage. This may not be obvious in the original unscaled data.

The scaling technique can also be applied to RT instability. By expressing (3.27) in terms of the scaled dimensionless variables

we obtain a universal curve for the scaled velocity in RT instability:

where $r=[1-u_{RT}(0)]/[1+u_{RT}(0)]$ .

Figure 7.

Comparison between universality curve given by (4.4) and scaled data from figure 3 for RM instability at different Atwood numbers. The solid curves are the scaled spike velocities, and the dashed curves are the scaled bubble velocities. Red curve: universality model prediction. Blue curves: numerical results (FLASH) from Dimonte & Ramaprabhu (Reference Dimonte and Ramaprabhu2010), with $a_{0}k=0.125$ and $A=\pm 0.25,\pm 0.5,\pm 0.75,\pm 0.88,\pm 0.94,\pm 0.98$ .

Figure 8.

Comparison between universality curve given by (4.4) and scaled data from figure 4 for RM instability at different Atwood numbers. The solid curves are the scaled spike velocities, and the dashed curves are the scaled bubble velocities. Red curve: universality model prediction. Blue curves: numerical results (LEEOR2D) from Alon et al. (Reference Alon, Hecht, Ofer and Shvarts1995), with $v_{0}=1$  cm ms $^{-1}$ , ${\it\lambda}=1$  cm and $A=\pm 0.2,\pm 0.5,\pm 0.9$ .

Figure 9.

Comparison between universality curve given by (4.6) and scaled data from figure 5 for RT instability at different Atwood numbers. The solid curves are the scaled spike velocities, and the dashed curves are the scaled bubble velocities. Red curve: universality model prediction. Blue curves: vortex model results from Sohn (Reference Sohn2004), with $a_{0}k=0.5$ , $v_{0}=0$ and $A=\pm 0.05,\pm 0.3,\pm 0.7,1$ . Black curve (labelled as c.m.): conformal mapping result from Menikoff & Zemach (Reference Menikoff and Zemach1983), with the same initial conditions and $A=1$ .

Figure 9, the scaled version of figure 5, also shows the similarity and the universality between the growth rate of spikes and that of bubbles at different Atwood numbers for RT instability. Despite the different developments in the early stage, the scaled physical quantities for all fingers are gradually approaching the same universality curve. The agreement between the theoretical prediction and the scaled numerical data for RT instability shown in figure 9 is not as good as those for RM instability shown in figures 6–8. This is due to several reasons. (1) The fingers in RT instability are driven by gravity, and therefore are more unstable than those in RM instability. Thus, it is much harder to carry out the simulations for RT instability to the asymptotic regime. They need a larger computational domain and longer computing time. Owing to these difficulties, in terms of scaled time, the numerical simulations for RT instability were actually terminated earlier than those for RM instability. (2) Since in terms of scaled time the simulations for RT instability are shorter than those for RM instability, they contain more transient effects, which are not included in our model. (3) In the numerical simulations presented by Sohn (Reference Sohn2004), a numerical desingularization parameter is introduced to avoid singularities in computation, which reduces the growth rates. Therefore, the numerical data are consistently below our theoretical prediction. For the same reason, the numerical results shown in figure 6 are also slightly lower than our theoretical prediction. In particular, figure 9 shows that the simulation for the spike at $A=-0.7$ ends too soon and has not approached the quasi-steady stage yet. Therefore, we call for further numerical studies of RT instability in late times, especially for spikes with large Atwood numbers.