A unified theoretical model for spatiotemporal development of Rayleigh–Taylor and Richtmyer–Meshkov fingers

Abstract

An initially perturbed interface between two fluids of different densities is usually unstable when driven by an acceleration or a shock wave; it is known as a Rayleigh–Taylor instability or a Richtmyer–Meshkov instability. One of the most significant issues in these instabilities is the spatiotemporal development of fingers generated at the interface, which plays an important role in both scientific research (e.g. supernova explosion) and engineering applications (e.g. inertial confinement fusion). Accurate theoretical solution of these interfacial fingers remains as an unsolved and challenging problem since Taylor's seminal work more than seven decades ago. This paper reports a unified theory established for such phenomena by combining the classical potential-flow theory and a dual-source model to address the long-standing difficulty highlighted by the initial-value sensitivity and strong nonlinearity. It is the first time for a theory to accurately predict the long-time developments in both growth rate and shape curvature of interfacial fingers at all density ratios in two and three dimensions. Moreover, the new theory clearly reveals the nonlinear coupling mechanism for interfacial evolution, and especially explains the origin of overshot in the growth rate curve.

1. Introduction

The initially perturbed interface of two fluids is usually unstable when driven by an acceleration ${g}$ pointing from a heavy fluid to a light one or impinged by a shock wave; they are known as a Rayleigh–Taylor instability (RTI) (Rayleigh 1883; Taylor 1950) and a Richtmyer–Meshkov instability (RMI) (Richtmyer 1960; Meshkov 1969), respectively. With the development of such hydrodynamic instabilities, light and heavy fluids penetrate into each other, resulting in the formation of bubble- and spike-like finger structures (Zhou 2017a,b). These phenomena play important roles in both scientific research and engineering applications, such as supernova explosions (see, e.g. Adam 2000; Isobe et al. 2005; Cabot & Cook 2006; Uchiyama et al. 2007) and inertial confinement fusion (see, e.g. Clery 2015; Betti & Hurricane 2016; Casey et al. 2017; Ding et al. 2017; Matsuo et al. 2021; Sabet et al. 2021). Two recent survey studies provide a systematic review on RTI and RMI phenomena, which make the corresponding applications accessible to a broad and varied audience (Zhou et al. 2019, 2021). An important issue of these classical hydrodynamic instabilities is how these finger structures on an unstable interface grow in time and space. There have been sets of experimental (Jacobs & Catton 1988; Luo, Wang & Si 2013; Liu et al. 2018; Luo et al. 2018; Liang et al. 2019; Lherm et al. 2022) and numerical (Ristorcelli & Clark 2004; Ramaprabhu, Dimonte & Andrews 2005; Li et al. 2021; Yan et al. 2022) data available for their growth rate and shape evolution since Taylor's pioneering work in 1950, and yet there is still a lack of a unified and accurate theory to characterize them.

The classical work of Rayleigh (1883) on hydrodynamic instability, along with its modern development by Taylor (1950), has been expanded by several milestone-like theories based on Layzer's potential-flow model (LPM) (see, e.g. Layzer 1955; Hecht, Alon & Shvarts 1994; Alon et al. 1995; Qiang & Sohn 1996; Mikaelian 1998; Zhang 1998; Goncharov 2002; Abarzhi, Glimm & Lin 2003; Mikaelian 2003; Sohn 2003; Zhang & Guo 2016; Zhang, Deng & Guo 2018; Guo & Zhang 2020; Zhao et al. 2020). These studies show that the interface evolution of single-mode RTI usually undergoes a linear stage, a nonlinear stage, quasisteady stages with constant growth rate, and other possible stages. Specifically, Mikaelian (1998) and Zhang (1998) are among the first to give the evolutions of bubble and spike, respectively, in a fluid–vacuum configuration. It is Goncharov (2002) who generalizes this method to an arbitrary density ratio. Recently, Zhang & Guo (2016) and Guo & Zhang (2020) obtained a universal finger evolution for arbitrary density ratio and Zhao et al. (2020) further applied it to RTI in a cylindrical geometry.

In the framework of original LPM, the sinusoidal potential function is employed to describe the flow of the fingertip and its immediate neighbourhood, which only approximately satisfies the dynamic condition of constant surface pressure. Meanwhile, most of these theories usually assume that the shape of the fingertip does not change with time in prediction of the growth rate, which is in conflict with numerical simulations from Sohn (2004) and Ramaprabhu & Dimonte (2005) and causes unphysical predictions in growth rate (Krechetnikov 2009). According to Kull (1983), the correctness of such potential functions without harmonics and possible singularities is questionable or hard to justify theoretically. In addition, it is found that the gravitational acceleration obtained based on the sinusoidal potential model increases rapidly below the fingertip and the increase trend becomes stronger for larger initial amplitudes (Kull 1986). This indicates that the acceleration in the flow field given by the sinusoidal potential flow model is non-uniform, which is obviously different from the fact that the acceleration, as a conservative force, should be uniform in the whole flow field. Since the growth rate depends weakly on the acceleration non-uniformities, the original LPM proves rather successful in describing the growth rate compared with the result from numerical simulation. However, the original LPM fails to capture the curvature of the fingertip accurately, which is a more sensitive parameter and is found to be only approximately ${2/3}$ of the value obtained numerically at the quasisteady stage (Kull 1986). Different from others, Zufiria (1988) allows time dependence of the shape based on a single-point source model introduced by Kull (1986), which can lead to a better approximation for shape, but only for a vacuum–bubble configuration. Sohn & Zhang (2001) also make use of such an assumption in comparison with the original LPM in a fluid–vacuum configuration. Nevertheless, no further results are secured due to mathematical complexity. Therefore, there still lacks a well-unified theory that can simultaneously predict the growth rate and shape evolution from early to late stages and is suitable for all density ratios, which motivates the present work.

This paper is intended to develop a new theory to address the long-standing problem within the LPM framework. The present theory is established qualitatively by introducing a dual-source model to mimic the pattern of real flow, and is further formulated quantitatively with the constraint of the geometric similarity principle. The derived equations can simultaneously predict the growth rate and shape evolution. It should be stressed that the new theory focuses on the development of instability up to the quasisteady stage similar to most of previous studies. The results are in good agreement with numerical and experimental data in both RTI and RMI regimes at all density ratios in different dimensions. Moreover, the theoretical model can help reveal the underlying nonlinear coupling mechanism between growth rate and shape evolution, and explain the origin of the overshot phenomenon in variation of the perturbation growth rate.

The remaining paper is structured as follows. First, theoretical analysis and derivations are given in § 2. Then, the dual-source model is validated for both RTI and RMI flow regimes in § 3. Finally, conclusion and discussion are made in § 4. For completeness, the detailed formulations for the present new theory are provided as Appendix A.

2. Theoretical analysis and formulations

For the purpose of clarity in subsequent discussions, we begin with an introduction of the LPM. In the two-dimensional (2-D) case, we consider two incompressible inviscid fluids in a vertically infinite strip with its left and right boundaries satisfying a no-penetration condition. The initial interface is characterized by a single-mode sinusoidal perturbation. The LPM for RTI problem can be described by the following governing equations (Goncharov 2002):

(2.1)\begin{gather} {\nabla ^2}\varphi _{i} = 0, \end{gather}

(2.2)\begin{gather}\dot{\eta} - \varphi _{i,x} \eta_x + \varphi _{i,z} = 0, \end{gather}

(2.3)\begin{gather}\sum_{i = 1}^2 {{{( - 1)}^i}} {\rho _i} \left[{-}g\eta + \dot{\varphi}_i -\frac{1}{2}({\varphi _{i,x}^2}+{\varphi _{i,z}^2})\right]=f(t). \end{gather}

Here $x$ and $z$ are the coordinates perpendicular and parallel to the acceleration ${g}$, respectively. Subscripts ${i=1,2}$ denote the upstream and downstream fluids, respectively. Here ${\varphi _{i}}$ is the velocity potential function, and ${f(t)}$ only depends on time $t$. According to the LPM, the shape of the fingertip is approximated as a parabola of the form ${\eta ( x,t )=z_0(t) + \xi (t)x^2}$, where ${z_0(t)}$ and ${\xi (t)}$ are the vertical position and curvature of the finger, which are used to characterize the growth rate and shape evolution, respectively. It is worth stressing that the traditional LPM only considers behaviour of the single finger located in the middle of the channel, and assumes that the fingers near the two sides of the strip have no dynamic effect on the middle one, as shown in figure 1(a). Furthermore, the LPM only describes the flow in the neighbourhood of the interface fingertip instead of the whole interface because the vorticity generated in the flow field may invalidate the potential function hypothesis far away from the fingertip. In three-dimensional (3-D) geometry, the above process can be easily extended to a cylindrical coordinate system with the $z$ axis in the direction of density gradient and the $x$ axis interpreted as the radial coordinate by assuming cylindrical symmetry of the fingers.

Figure 1. (a) Schematics of the dual-source model. The finger is generated by the penetration of downstream fluid ${(\rho _{1})}$ into the upstream fluid ${(\rho _{2})}$. Two sources are separated by the finger with each being located at a distance ${a_{i}(t)}$ from the fingertip with an intensity ${Q_{i}(t)}$. Here ${\lambda _{i}}$ is the characteristic wavelength of the potential function ${\varphi _{i}}$. The red solid parabola is the middle finger approximated by LPM and the blue dashed parabolas are fingers neglected in LPM. The dashed lines with arrows are streamlines induced by the present model. (b) Streamlines of the flow field obtained from the present numerical simulation of RTI ($A=0.15$), which are overlaid on contours of densities.

For the LPM, the key to accurately solving the interfacial fingers is to construct a reasonable velocity potential function ${\varphi _{i}}$. As for the treatment in the original LPM, the simplest sinusoidal potential function ${\varphi _{i}=\sum _{n} a_{i}^{n}(t)\cos (nkx)\exp ({(-1)^{i}nkz})}$ is introduced to describe the flow immediately around the fingertip with the constraint of approximately constant surface pressure. This form of ${\varphi _{i}}$ proves to capture the evolution of finger growth rate well, but the predicted curvature can deviate strongly from the numerical result (especially at the quasisteady stage), which is ascribed to the violation of uniform acceleration. Although efforts have been taken to improve the prediction ability for uniform acceleration and finger curvature by introducing source(s) (Kull 1986; Zufiria 1988; Sohn & Zhang 2001), few advances are made beyond the fluid–vacuum configuration. Therefore, suitable potential functions ${\varphi _{i}}$ for all density ratios are expected to describe the development of the interface finger.

Previous experiments and numerical simulations show that the finger development is closely related to the initial conditions. Therefore, a good theory should consider and solve the issue of initial-value sensitivity. Furthermore, it is observed that the growth rate and shape evolution of fingers are coupled with each other, and their strong nonlinear interaction should also be taken into account. These two challenges may be the major obstacles to develop an accurate theory for classical hydrodynamic instabilities. As discussed above, however, no available theory satisfies these requirements to date. To this end, we develop a novel theory by including the evolution of the fingers’ shape. Through the physical image of RTI (see figure 1b), we are surprised that there are always dual-source singularities combined with a uniform stream on both sides of the finger. Both of them are generated at the initial time and only change the distance from the interface due to the temporal variation of source intensity. Thus, the streamlines in a flow field should always meet the condition of dual-source no matter how the initial perturbation changes. Moreover, the initial effect gradually decays with the further development of fingers, and the strong nonlinear coupling between growth rate and shape evolution begins to dominate. Although the finger undergoes a nonlinear development, the fingertip structure always maintains a self-similar pattern. With these understandings, our approach consists of two steps: (i) construct a dual-source field by simulating the initial flow pattern to satisfy the initial-value sensitivity; and (ii) maintain the dual-source characteristics to guarantee the self-similar evolution of the interface and further restrict the nonlinear interaction. Based on the assumption of LPM and the numerical flow field (figure 1b), the RTI flow field is simplified as that depicted in figure 1(a). It is assumed that the flow field near the finger is generated by two sources distributed on the upstream and downstream sides of the interface, respectively. Each source is located in a distance ${a_{i}(t)}$ from the fingertip with source intensity ${Q_{i}(t)}$. It should be stressed that the potential functions ${\varphi _{i}}$ must obey the following constraints: (i) being a solution of (2.1); (ii) satisfying the symmetric and far-field boundary conditions as sketched in figure 1(a); (iii) avoiding an imaginary component in the solution (Goncharov 2002). By taking these constraints into consideration and referring to the research works by Zufiria (1988) and Sohn & Zhang (2001), a new potential function ${\varphi _{i}}$ for the 2-D case can be obtained in the form

(2.4)\begin{equation} \varphi _{i}= {\rm Real}\left\{ {{Q_i}(t)\ln {\frac{{\sinh [(z + jx + {a_i(t)}){k_i}/2]}}{{\exp [(z + jx + {a_i(t)}){k_i}/2]}}} } \right\}, \end{equation}

where ${j^2=-1}$, the hyperbolic part (${\sinh }$) represents the source flow and the exponential part ($\exp$) represents the uniform flow. For comparison purposes, the streamlines for the pure source flow (induced by the hyperbolic part) and combined flow (induced by dual-source singularities in combination with a uniform stream) are shown in figures 2(a) and 2(b), respectively. Here, the parameters are set to $Q_{i}=1$, $k_{i}=1$ and $a_{i}=2(-1)^{i}$ without loss of generalization. From figure 2(b), one can clearly see that the patterns of combined flow are quite similar to those observed for the simulated RTI flow both upstream and downstream of the interface as shown in figure 1(b). It needs to be mentioned that the flow patterns in figure 2(b) are only a qualitative display of the potential function (2.4), and the realistic flow field should be further constrained by (2.2). Nevertheless, it is shown that the potential function (2.4) can successfully reflect typical characteristics of the target flow field. Since there is no concise complex potential function for 3-D flow, the series forms of this new potential function ${\varphi _{i}}$ are employed to facilitate the unified solution, which reads

(2.5)\begin{gather} \varphi _{i}= Q_{i}(t)\sum_{n=0}^{\infty}\frac{\cos\left[(2n+1)k_{i}x\right]\exp({({-}1)^i(2n+1)k_{i}(z+a_{i}(t))})}{2n+1} , \end{gather}

(2.6)\begin{gather}\varphi _{i}= Q_{i}(t)\sum_{n=0}^{\infty}\frac{{\rm J}_0\left[(2n+1)k_{i}x\right]\exp({({-}1)^i(2n+1)k_{i}(z+a_{i}(t))})}{2n+1} , \end{gather}

where (2.5) and (2.6) are the expressions for the 2-D and 3-D cases, respectively. Here ${k_{i}\equiv 2{\rm \pi} /\lambda _{i}}$ is the wavenumber, and ${{\rm J}_{0}(x)}$ is the zero-order Bessel function. For a specific ${\varphi _{i}}$, one can readily obtain the flow field in the vicinity of the interface as illustrated by the streamlines (dashed lines) in figure 1(a). As can be seen, the streamlines are consistent with the current numerical simulation shown in figure 1(b), which implies the rationality of the present dual-source figure for ${\varphi _{i}}$. It should be mentioned that the present model readily recovers previous potential flow models (e.g. Zhang 1998; Goncharov 2002; Zhang & Guo 2016) when the expansion order $n\leq 2$ and wavenumbers $k_1=k_2$. Meanwhile, the new velocity potential function ${\varphi _{i}}$ involves only two time-dependent unknown quantities ($Q_{i}(t),a_{i}(t)$) due to the dual-source constraint.

Figure 2. The streamlines for flow induced by (a) pure sources (i.e. the hyperbolic part in (2.4)) and (b) source-uniform flow combination (dual-source singularities combined with a uniform stream, i.e. both the hyperbolic and exponential parts in (2.4)).The parameters are taken as $Q_{i}=1$, $k_{i}=1$ and $a_{i}=2(-1)^{i}$ for reference.

As one can imagine from figure 1(a), the channel should not be completely occupied by the middle finger (approximated by a parabola) due to the existence of fingers near both channel walls (which are neglected in LPM). In the real physical picture, when the convex curve (the red solid parabola) and concave curve (the blue dashed parabola) intersects, there are always two points where the tangent line is parallel to the channel wall, and the non-penetrability condition of the interface is satisfied. The distance between these two points ${\lambda _{1}}$ represents the characteristic scale of the finger. Therefore, according to the parabolic approximation of LPM, the velocity potential function ${\varphi _{1}}$ should be characterized by ${\lambda _{1}}$ (corresponding to wavenumber ${k_{1}}$), while ${\varphi _{2}}$ should always reflect the scale of channel width ${\lambda _{2}}$ (corresponding to wavenumber ${k_{2}}$). For such a purpose, the potential function on each side of the finger does not have to share the domain. Note that the wavenumbers ${k_{1}}$ and ${k_{2}}$ in the potential functions are not employed to reflect the lateral scale of the whole interface, but to characterize the local fluid motion near the fingertip. In addition, the results from numerical simulations and experiments show that the evolution of fingers are self-similar. Qualitatively, it implies that the fluid between the source and fingertip approximately satisfies the geometric similarity, which leads to ${\lambda _{1}^{0}a_{1}^{0}=\lambda _{1}a_{1}}$, with the superscript ‘0’ representing the initial time. Based on this assumption, one ends up with the following relation for the ratio of wavenumbers:

(2.7)\begin{equation} c\left( {A,t} \right) \equiv \frac{{{\lambda _1}}}{{{\lambda _2}}} = \frac{{{k_2}}}{{{k_1}}} = \frac{{\lambda _1^0}}{{{\lambda _2}}}\frac{{a_1^0}}{{{a_1}}} = {c^0}\frac{{a_1^0}}{{{a_1}}}. \end{equation}

In light of the numerical results (Sohn 2004), the wavelength ratio is close to 2 ($\lambda _{2}/\lambda _{1}\sim 2$) when the density ratio is very close to 1, and the shapes of spike and bubble are almost the same. However, the wavelength ratio is close to ${1 (\lambda _{2}/\lambda _{1}\sim 1)}$ when one considers a fluid–vacuum configuration (Layzer 1955), which corresponds to the extreme case. For the latter, the channel is almost occupied by the bubble, and the spike almost disappears. Therefore, it is concluded that ${c^0}$ is a function of the density ratio and falls within the interval ${c^0\in [0.5,1]}$. In consideration of the extreme cases of ${A=1}$ (fluid–vacuum configuration) and ${A=0}$ (equal density) as well as the dimensional effects, it is suggested that ${c^0(A)=(A^3+1)/2}$ for the 2-D case and ${c^0(A)=(A^{1.5}+1)/2}$ for the 3-D case. Here, the Atwood number ${A\equiv ( \rho _{1} - \rho _{2} )/ ( \rho _{1} + \rho _{2} )\in [ -1,1 ]}$ is employed to realize a unified description (Zhang & Guo 2016). Specifically, one considers bubbles (light fluid penetrating into heavy fluid with a downward gravity) for ${A>0, g>0}$, and spikes (heavy fluid penetrating into light fluid with a upward gravity) for ${A<0, g<0}$.

From the above analysis, there are six unknowns in the dual-source model, i.e. $z_0$, $\xi$, $Q_{1}$, $Q_{2}$, $a_{1}$, $a_{2}$. Following the LPM, four constraints can be obtained by substituting expansions of $\varphi _{1}$ and $\varphi _{2}$ into (2.2) and expanding the resultant equations to the order of ${x^2}$. In order to close the current model, one needs to expand equation (2.3) through the orders of ${x^2}$ and ${x^4}$, respectively, to yield another two constraints. Finally, one ends up with the following ordinary differential equations:

(2.8)\begin{gather} \dot{\xi} = F\left( {{a_1},{k_1},\xi ,v} \right), \end{gather}

(2.9)\begin{gather}{S_1}g + {S_2}{v^2} + {S_3}\dot{v} + {S_4}v{\dot{a}_1} = 0, \end{gather}

(2.10)\begin{gather}{T_2}{v^2} + {T_3}\dot{v} + {T_4}v{\dot{a}_1} = 0. \end{gather}

Here, $v={\rm d}z_0/{\rm d}t$ is the growth rate, and ${F(a_{1},\xi,k_{1},v)}$ can be obtained from (2.2) in both the 2-D and 3-D cases. Here ${S_{p}(A,c,a_{1},\xi,k_{1})}$ and ${T_{p}(A,c,a_{1},\xi,k_{1})}$ $(p=1,2,3,4)$ are the coefficients of the ${x^2}$ and ${x^4}$ terms of (2.3). The detailed expressions of these coefficients are obtained with the aid of Mathematica (Wolfram Research Inc. 2017) and are provided in Appendix A for reference. From (2.8)–(2.10), it can be clearly seen that there is a strong nonlinear coupling effect between growth rate and shape evolution. Through simple manipulation, (2.8)–(2.10) can also be written as a general ternary ordinary differential equation (ODE) system of the form

(2.11)\begin{gather} \dot{\xi} = F\left( {{a_1},{k_1},\xi ,v} \right), \end{gather}

(2.12)\begin{gather}\dot{v} = \frac{S_1 T_4 g + (S_2 T_4 - S_4 T_2)v^2}{S_4 T_3 - S_3 T_4}, \end{gather}

(2.13)\begin{gather}\dot{a}_1 = \frac{-S_1 T_3 g + (S_3 T_2 - S_2 T_3)v^2}{S_4 T_3 - S_3 T_4}. \end{gather}

Therefore, one can accurately predict the development of fingers only by solving the above ternary ODE system. Again, it should be pointed out that, the settings of $A>0, g>0$ in (2.11)–(2.13) lead to the solution for bubbles. Accordingly, these parameters are set to $A<0, g<0$ if one wants to obtain the results for spikes. In particular, the above equations recover the behaviour of bubbles (or that of spikes) in an infinite-density-ratio system when $A=1$ (or $A=-1$).

3. Validation of the dual-source model

3.1. Rayleigh–Taylor instability

The dual-source theory is validated in comparison with the numerical simulations (2-D) carried out by Sohn (2004) for RTI, and the theoretical predictions given by Goncharov (2002) (2-D and 3-D) are also presented to highlight the geometric effect of shape evolution. According to the above references, the parameters and initial values are $g=1$, $k _{1}=1$, $\xi ^{0}=-0.25$ (2-D), $\xi ^{0}=-0.15$ (3-D) and $v^0=0$. Now ${a _{1}^{0}}$ can be solved by setting $F=0$ in (2.8), which comes from the constraint ${\dot {\xi }=0}$ at the initial time. For convenience of comparison, a scaled dimensionless time is introduced in the form ${\tau =\sqrt {Agk}t}$, which is the same as that suggested by Goncharov (2002). Shown in figures 3(a) and 3(b) are evolutions of the bubble fingertip velocity, figures 3(c) and 3(d) are the corresponding bubble curvature and figures 3(e) and 3(f) are the corresponding spike velocity for RTI at ${A=0.3}$ and ${A=0.7}$, respectively. Note that Goncharov (2002) only provides the theoretical prediction for bubbles. Therefore, the corresponding theoretical curves for spikes are absent in figures 3(e) and 3(f). Likewise, the results under the same conditions for 3-D cases are plotted in figures 4(a)–4(d).

Figure 3. Evolutions of the RTI bubble velocity (a,b), corresponding curvature (c,d) and spike velocity (e,f) given by the dual-source model for 2-D case (red solid lines) in comparison with the numerical simulations (Sohn 2004) (blue dashed lines) and the theory by Goncharov (2002) (black dotted lines).

Figure 4. Evolutions of the RTI bubble and spike velocity (a,b) and bubble curvature (c,d) given by the dual-source model for the 3-D case (red solid and dashed lines) in comparison with the theory by Goncharov (2002) (black dotted lines).

For bubbles, the comparisons show that the predicted ${v(t)}$ and ${\xi (t)}$ almost collapse onto the data of previous simulations in the linear, nonlinear and quasisteady stages. For spikes, the modelled result is consistent with the numerical results (Sohn 2004) at $A=-0.3$ as seen in figure 3(e), but a visible difference exists between the present theory and numerical simulation at $A=-0.7$ as observed in figure 3(f) (with the relative error up to $10.7\,\%$). On the one hand, it seems that the growth rate has not entered the quasisteady stage for the case $A=-0.7$ in comparison with the numerical simulation results for other values of $A$. Due to the difficulty encountered in the numerical simulation of high-(negative) Atwood number RTI, the prediction of spike evolution in the late stage may not be accurate. Nevertheless, the evolution trend offered by numerical simulation is consistent with the present theoretical results. On the other hand, it is easy to understand that, with the increase of density ratio, the growth rate on the spike side is much more affected by vorticity than on the bubble side, which shortens the quasisteady stage until it disappears. Therefore, the framework of LPM can be questionable on the spike side at high-(negative) Atwood numbers, which may also account for the deviation of present theory from numerical simulation. To the best of our knowledge, this is the first time for a theory to achieve a complete prediction of RTI growth rate for both ${z_0(t)}$ and ${\xi (t)}$ with $A<0$, though an obvious difference is observed between the modelled and numerical results at high-negative density ratios (e.g. $A=-0.7$).

It is clearly seen from figures 3 and 4 that there exist great differences in bubble velocity and shape curvature between the present model and the Goncharov model. Such differences are mainly caused by the discrepancy in how the curvature is constrained in different modelling strategies. In light of existing experiments and numerical simulations, it is suggested that the strong nonlinear interaction between the growth rate and the shape evolution of fingers should be taken into account in a theoretical modelling procedure. Specifically, the expansion terms ${x^2}$ and ${x^4}$ of the momentum equation (2.3) are used in the present model to solve the curvature coupled with the growth rate, which leads to a coupled system of equations $\dot {\xi } = f( \xi,v ,A )$ and $\dot {v} = g( \xi,v ,A )$. By contrast, the Goncharov model applies the expansion terms ${x^2}$ of the interfacial kinematics equation of a heavy fluid (2.2) to solving the curvature only with heavy fluid, which leads to an equation of the form $\dot {\xi } = f( \xi,v )$. In other words, the curvature evolution in the Goncharov model is determined only by the heavy fluid and is not really constrained by the fluids on both sides.

It is also noticed that there exist local maximum or minimum values of bubble velocity (figure 4a,b) and shape curvature (figure 4c,d), which correspond to their critical times. As seen in figure 4, the critical times for bubble velocity and shape curvature are different at the same $A$. This is because the bubble velocity is directly affected by the first derivative of curvature $(\dot {\xi })$ rather than curvature ($\xi$) itself, and there exists a time delay for such a nonlinear interaction to produce an effect on the evolution of bubble velocity. It can be seen from figure 4(a) that the curve for bubble velocity given by the Goncharov model intersects with that from the present model at $\tau \approx 1.6$, which approximately corresponds to the critical time of the shape curvature in figure 4(c). This also indicates that the solution process of curvature is decoupled with that of velocity in the Goncharov model. It is also observed that the sign change in the first derivative of curvature ($\dot {\xi }$) before and after its critical time corresponds to the alternation from suppression to promotion of bubble velocity development. The same mechanism also applies to the situations in the 2-D cases (see figure 3).

Meanwhile, according to Goncharov (2002), the adoption of higher-order modes may alter the asymptotic value of $\xi$ with respect to time, which will converge rapidly with the increase of mode order. The Goncharov model only uses the standard perturbation expansion with mode $k$, while the present model uses higher-order modes $3k, 5k, \ldots\,$. Therefore, it is suggested that the present model may provide more accurate asymptotic values for both bubble velocity and shape curvature than Goncharov model. As seen in figures 4(c) and 4(d), the asymptotic values of $\xi$ given by the present model are a little bit smaller than those reported for the Goncharov model in the 3-D case. However, this discrepancy is more pronounced in the 2-D cases (see figure 3a–d).

Finally, it should be pointed out that the current theory can successfully predict the possible overshot of bubble velocity before the quasisteady stage. Specifically, our theory indicates that the growth rate induced by acceleration will also be affected by the geometric effect of the shape evolution. Therefore, the overshot is essentially caused by the change of curvature. More importantly, according to the present theoretical solution, its amplitude is completely determined by the initial shape. Therefore, it is believed that consideration of the initial-value sensitivity and strong nonlinear coupling are the key to success of this new theory.

The quasisteady stage is extremely important for RTI flow (Clery 2015; Betti & Hurricane 2016; Casey et al. 2017; Matsuo et al. 2021; Sabet et al. 2021; Ding et al. 2017). In this stage, the fingertip velocity and curvature are nearly independent of time, i.e. ${\dot {v}=0}$, ${\dot {\xi } =0}$ and ${\dot {a}_1(v,\xi )=0}$. In such a case, (2.8)–(2.10) can be simplified as

(3.1)\begin{gather} F\left( {{a_1},{k_1},\xi^{qs} ,v^{qs}} \right) = 0, \end{gather}

(3.2)\begin{gather}{S_1}g + {S_2}{(v^{qs})^2} = 0, \end{gather}

(3.3)\begin{gather}{T_2} = 0. \end{gather}

Thus, the quasisteady velocity ${v^{qs}}$ and curvature ${\xi ^{qs}}$ of fingers can be algebraically solved from (3.1)–(3.3). Figures 5(a) and 5(c) compare the results within a wide range of $A$ for the current theory, Sohn's simulations and the Goncharov's theory in the 2-D and 3-D cases, respectively. It is shown that the quasisteady velocities of bubble/spike fingers predicted by the present theory exhibit excellent agreement with numerical simulations conducted by Sohn (2004), Ramaprabhu & Dimonte (2005) and Ramaprabhu et al. (2012), respectively. A functional form of ${v^{qs}}$ can be obtained by applying the least-squares method to the present data, which is formally similar to the formulation given by Goncharov (2002) and reads

(3.4)\begin{equation} {v^{qs}} \approx f(A)\sqrt {\frac{A}{{1 + A}}\frac{g}{k}}, \end{equation}

where $f_{2-D}(A)=0.814-0.027A$ and $f_{3-D}(A)=1.456 - 0.01A$ for the 2-D and 3-D cases, respectively. The new theory can also predict the variations of ${\xi ^{qs}}$ with $A$ very well in comparison with the numerical simulations by Sohn (2004), which are depicted in figures 5(b) and 5(d). Except for extreme cases, it is demonstrated that ${\xi ^{qs}}$ is almost independent of $A$, i.e. approximately ${-0.25}$ in the 2-D case and ${-0.14}$ in the 3-D one. In contrast, previous theories such as Goncharov's second-order theory suggests that $\xi ^{qs}(A)\equiv -1/6$ for the 2-D case and $-1/8$ for the 3-D one, which deviate approximately 33 % and 16.7 % from numerical simulations, respectively (see dashed lines). Although the theory developed by Goncharov (2002) can be further improved when expanded to infinite order (dotted lines in figure 5b,d), the correction of the result by the infinite-order theory as such is still limited.

Figure 5. (a,c) The RTI velocity and (b,d) bubble curvature at the quasisteady stage for all Atwood numbers in both 2-D and 3-D cases (red solid lines). The results from numerical simulations by (a,b) Sohn (2004) and (c,d) Ramaprabhu & Dimonte (2005) and Ramaprabhu et al. (2012) (blue solid circles) are included for comparison. The black dashed and short-dashed lines represent Goncharov's second-order and infinite-order theories, respectively.

3.2. Richtmyer–Meshkov instability

So far, the newly proposed theory can offer a near-perfect prediction of the single-mode RTI evolution at all density ratios in both the 2-D and 3-D cases. The new model is also suitable for the theoretical prediction of RMI. This is due to the fact that RMI is approximately equivalent to an incompressible RTI without an external force after the impingement of a shock wave (Richtmyer 1960; Hecht et al. 1994). Therefore, the evolutions of ${z_0(t)}$ and ${\xi (t)}$ in the RMI problem can be readily obtained by solving the same equations (2.8)–(2.10) if one sets the initial time to the shock arrival time and employs the same potential functions for RTI flow. In such a situation, (2.8)–(2.10) can be reduced to the following form:

(3.5)\begin{gather} \dot{\xi}= F\left( {{a_1},{k_1},\xi ,v} \right), \end{gather}

(3.6)\begin{gather}{S_2}{v^2} + {S_3}\dot{v} + {S_4}v{\dot{a}_1} = 0, \end{gather}

(3.7)\begin{gather}{T_2}{v^2} + {T_3}\dot{v} + {T_4}v{\dot{a}_1} = 0. \end{gather}

Nevertheless, there exist some subtle differences in the settings of $g$ and $a_{1}^{0}$ for RMI compared with those for RTI. In general, ${g}$ should be zero in the LPM theory. In a real situation of RMI, although no external force exists after the shock impact, there is still a small acceleration at the finger. Therefore, for comparability between theory, simulation and experiment, it is more reasonable to allow a small ${g}$ in practice. Moreover, ${a_{1}^{0}}$ should also be determined according to the relevant initial conditions after shock.

The proposed dual-source model is then validated for evolution of RMI in comparison with numerical simulations and experimental measurements. The results from a recently developed theoretical model by Zhang & Guo (2022) are also presented to gain insight into the geometric effect of shape evolution. When the present model is used to solve RMI problems, the required parameters are ${v^0}$, ${\xi ^0}$, ${a_{1}^{0}}$, which are essentially the initial values after the shock wave. In practice, it is suggested that $k _{1}=1$, ${v^0 = 1}$, $\xi ^{0}=-0.25$ (2-D) according to the references, and $\xi ^{0}=-0.15$ (3-D) according to the conclusion in the RTI problems (since no specific values are given in the literature for RMI flow). For comparison purposes, a scaled dimensionless time is introduced in the form ${kv^0t}$, and the initial velocity for numerical simulations and experimental measurements are given by introducing the scaled dimensionless form as ${v/v^{0}=1}$, which is consistent with the suggestion by Zhang & Guo (2022). The other settings are the same as those in the aforementioned RTI.

Shown in figures 6 and 7 are evolutions of the fingertip velocity and curvature for both 2-D and 3-D RMI at different Atwood numbers. Since the theoretical model given by Zhang & Guo (2022) assumes that the curvature is a time-invariant parameter, the corresponding results are not shown in figures 6(b) and 6(d). It should also be mentioned that only velocity evolutions are displayed for the 3-D cases due to the availability of reference data for comparison from experiment (Chapman & Jacobs 2006) and numerical simulation (Long et al. 2009). The values of ${a _{1}^{0}}$ and ${g}$ in our model are calibrated and provided in the figure. The comparisons show that the present theory provides very good predictions for the RMI finger velocity. As can be seen in figures 6(a) and 6(c), the theoretical results provided by both Zhang & Guo (2022) and the present paper agree well with the numerical results for bubbles. All the theoretical predictions and numerical results show that the bubble growth rate decreases monotonically with time. It is shown that the coupling effect between growth rate and curvature is relatively weak, and the coupling level proves to be different at various Atwood numbers. For spikes, the results given by Zhang & Guo (2022) and the present work are still in fairly good agreement with the numerical results. It is noted that the theoretical curves may underpredict or overpredict the numerical results within certain time ranges, especially for the small-Atwood number case. Specifically, both theories can capture the early-time peak phenomenon in spike growth rate for the large-Atwood number case, except that the peak value is apparently overestimated and the occurrence of such a peak is delayed in the theoretical model offered by Zhang & Guo (2022) (see figure 6c). The success in capturing the spike growth rate of the two models is attributed to the knowledge that the spike curvature has a strong impact on the growth rate. The difference is that the growth rate and curvature are coupled together in the present model, while an asymptotic spike curvature is employed in the scenario introduced by Zhang & Guo (2022). As for the finger curvature, slight deviations can be observed between the current predictions and numerical simulations by Sohn (2004), especially in the later stage. This is ascribed to the fact that the shock-induced interface acceleration in a realistic situation (also in simulations) should vary with time, while the present model takes a small constant ${g}$ instead for simplicity.

Figure 6. Evolutions of the RMI velocity (a,c) and bubble curvature (b,d) for the 2-D cases. The red solid and dashed lines denote the results for bubble and spike, respectively, given by the dual-source model. The black dashed and dash–dotted lines represent the 2-D simulations (Zhang & Guo 2022) for bubble and spike. The blue open circles and crosses represent the 2-D simulations (Sohn 2004) for bubble and spike, respectively. The specific values of ${a _{1}^{0}}$ and ${g}$ are given in panels (a,c).

Figure 7. Evolutions of the RMI velocity (a,b) for the 3-D cases. The red solid and dashed lines denote the results for bubble and spike, respectively, given by the dual-source model. The blue open circles and crosses denote the 3-D experiment (Chapman & Jacobs 2006) and numerical simulation (Long et al. 2009) for bubble and spike, respectively. The specific values of ${a _{1}^{0}}$ and ${g}$ are given in the figure.

4. Conclusion and discussion

In summary, a new analytical theory is developed for interfacial fluid instabilities by introducing a dual-source into the potential function, which is distinguished from previous scenarios within the LPM framework. The new theory is established qualitatively by introducing a dual-source model to mimic the pattern of real flow, and is further formulated quantitatively with the constraint of the geometric similarity principle. It is believed that the newly proposed methodology can better meet the requirement of the initial-value sensitivity in both RTI and RMI. In addition, the strong nonlinear interaction between the growth rate and shape evolution of fingers is also taken into account. These two challenges may be the major obstacles to develop an accurate theory for classical hydrodynamic instabilities. By introducing the double-source potential function into the governing equations, the derived ternary ODE system can simultaneously capture the growth rate and shape evolution. The predicted finger positions and curvatures agree well with the reference numerical and experimental data for both planar RTI and RMI. The theoretical predictions nearly collapse onto numerical simulations at the linear, nonlinear and quasisteady stages of RTI and the whole process of RMI. Moreover, it can be applied to both RTI and RMI at all density ratios in different dimensions, which confirms the validity of our theory. Furthermore, the slight overshot of finger velocity before the quasisteady stage is captured by our model for the first time, which benefits from its inclusion of the nonlinear coupling effect between shape evolution and growth rate.

We comment that the growth rate and shape curvature are equally important and should be involved in the theoretical modelling of classical hydrodynamic instabilities. The former represents the movement evolution of the fingertip, and the latter stands for the geometric changes of the fingertip and its neighbourhood. The present results show that the curvature plays a crucial role in the overshot phenomenon during the velocity evolution, which is consistent with existing numerical simulations. The traditional methods, however, cannot predict this phenomenon. Therefore, accurate capture of curvature evolution is crucially important (i.e. the geometric change of the fingertip) if one wants to predict the evolution of growth rate well. Specifically, ‘spatio’ in the article title indicates that the fingertip changes with the geometry of its neighbourhood, while ‘temporal’ implies that the fingertip also changes with time. These two parts must be satisfied simultaneously to accurately solve the fingertip motion.

Further, the mechanism discovered in this letter provides the possibility of controlling interfacial instabilities. It is also anticipated that the new theory should be extended to regimes with arbitrary acceleration (Dimonte, Ramaprabhu & Andrews 2007; Mikaelian 2009, 2010; Ramaprabhu et al. 2016). As a first step, one can replace the constant acceleration in the current theory with a time-varying acceleration. Thus, the new dual-source theory is expected to find applicability in better understanding the variable acceleration problems. For RTI and RMI problems, the dynamics in the mixing zone depends on the fingertip growth rate. The present work sheds light on the mechanisms of finger growth rate for single-mode RTI and RMI to some extent. If insight is further taken into the behaviour of finger interaction, it is possible to better understand the evolution of the resultant turbulent mixing.