1. Introduction
The Richtmyer–Meshkov instability (RMI) is a fundamental interfacial instability in hydrodynamics. It develops whenever a planar shock wave interacts with a material interface, separating two fluids with distinct thermodynamic properties. First predicted through numerical modelling by Richtmyer (Reference Richtmyer1960) and subsequently experimentally observed by Meshkov (Reference Meshkov1969), this phenomenon holds paramount significance across various realms of physics and astrophysics (see Inogamov Reference Inogamov1999; Brouillette Reference Brouillette2002; Aglitskiy et al. Reference Aglitskiy2010; Abarzhi, Gauthier & Sreenivasan Reference Abarzhi, Gauthier and Sreenivasan2013; Clavin & Searby Reference Clavin and Searby2016; Zhou Reference Zhou2017a,Reference Zhoub; Zhou et al. Reference Zhou2021; Zhou Reference Zhou2024; Zhou, Sadler & Hurricane Reference Zhou, Sadler and Hurricane2024) and references therein. In the domain of inertial confinement fusion and the associated area of high-energy-density physics (HEDP), the strong shocks propagate through condensed materials or high-pressure plasmas whose equations of state (EoS) under extreme conditions are typically non-ideal. Consequently, any theoretical framework aspiring to model HEDP conditions and replicate the RMI development must incorporate realistic fluid EoS.
Studies of interfacial instabilities begin with the small-amplitude, linear phase of their development, when the interfacial perturbations can be decomposed into Fourier series, and each Fourier component is dealt with separately. When the unperturbed state is time independent, which is the case for the Kelvin–Helmholtz (von Helmholtz Reference von Helmholtz1868; Thomson Reference Thomson1871) and the Rayleigh–Taylor (Rayleigh Reference Rayleigh1882; Taylor Reference Taylor1950) interfacial instabilities, it is possible to use the normal-mode decomposition to determine the exponential growth rate. The RMI analysis is more challenging because the unperturbed flow is not steady (it involves a transmitted shock wave and a reflected shock or rarefaction wave propagating away from the shocked material interface), and, therefore, ‘one cannot separate the time dependence from the linearised perturbation equations’ (Brouillette Reference Brouillette2002). This is why Richtmyer made his discovery by solving these equations numerically with the most powerful computer of the time, the MANIAC (mathematical analyser, numerical integrator and automatic computer) at Los Alamos. His linear analysis ‘places no restriction on the EoS, provided only that the values of several constants are known. But the only calculations performed with this code so far were for $\gamma$-law gases, initially cold, with the same value of $\gamma$’ (Richtmyer Reference Richtmyer1960). He was far ahead of his time: a code numerically solving the linear RMI perturbation equations was recreated only 40 years after Richtmyer's work by Yang, Zhang & Sharp (Reference Yang, Zhang and Sharp1994). It could handle arbitrary EoS, although its published output (see also Holmes et al. (Reference Holmes, Dimonte, Fryxell, Gittings, Grove, Schneider, Sharp, Velikovich, Weaver and Zhang1999)) still refers to ideal gases, generalising Richtmyer's work in several other ways (the case of reflected rarefaction wave, finite incident shock strength, different $\gamma$s.) The same applies to the analytical and semianalytical solutions of the linear RMI problem, starting from those published by Fraley (Reference Fraley1986) and Velikovich (Reference Velikovich1996) for the cases of reflected shock and reflected rarefaction, respectively.
Given the relative complexity of finding and using the exact solutions of the linear RMI perturbation problem, it is helpful, particularly for theoretical analysis (see Mikaelian (Reference Mikaelian1991, Reference Mikaelian1995), and references therein) to have an approximate but explicit formula for evaluating the constant asymptotic growth rate. Such a formula, also proposed by Richtmyer, is derived by treating the RMI as a limiting case of the Rayleigh–Taylor instability, with instantaneous, rather than constant, acceleration: $g=U_2\delta (t)$, where $U_2$ is the velocity instantly acquired by the shocked material interface, see figure 1, and $\delta (t)$ is the Dirac delta function. Then, the asymptotic growth rate is estimated as $\delta u_i^\infty =k \psi _i^0 U_2 \mathcal {A}_2$, where $\psi _i^0$ is the interfacial ripple amplitude, $\mathcal {A}_2$ is the Atwood number at the interface and $k$ is the interfacial ripple wavenumber. For an incompressible RMI flow triggered with a $\delta$-function acceleration pulse, this is indeed an exact result. Predictions of the linear compressible theory converge to it in the limits of low compressibility and/or weak shock, see the discussion by Velikovich, Herrmann & Abarzhi (Reference Velikovich, Herrmann and Abarzhi2014). However, suppose the compressibility of the materials separated by the shocked interface is relevant, which is the case in most strong-shock RMI situations. In that case, the preshock values of the Atwood number $\mathcal {A}_0$ and the ripple amplitude, $\psi _i^{0-}$, differ from the respective postshock values, $\mathcal {A}_2$ and $\psi _i^{0}$. Then, the governing incompressible equations do not help determine which values to substitute into the above formula, and the choice has to be made heuristically, resulting in so-called prescriptions. Richtmyer, comparing the above formula with his numerical results, formulated the prescription involving the postshock values. Meyer & Blewett (Reference Meyer and Blewett1972) similarly developed a slightly different prescription, featuring the average of preshock and postshock interfacial ripple amplitudes and the postshock Atwood number. Vandenboomgaerde, Mügler & Gauthier (Reference Vandenboomgaerde, Mügler and Gauthier1998) suggested yet another prescription, averaging the preshock and postshock products of the Atwood number and the ripple amplitude. Of course, in the incompressible and weak-shock limits all prescriptions converge to the same, exact value of the RMI growth rate.
Sometimes, the heuristic prescriptions provide an excellent approximation to the asymptotic RMI growth rate; sometimes, the approximation is not good at all, and, unfortunately, it is impossible to predict which way it is without solving the perturbation problem exactly. The physical reason for this uncertainty is that the fully compressible RMI growth rate is not settled at the instant when the rippled interface is shocked. As demonstrated by Wouchuk & Nishihara (Reference Wouchuk and Nishihara1997), the RMI growth rate is determined by two comparable contributions to the interfacial vorticity. One is the vorticity deposited at the instant of incident shock refraction. The other comes from the multiple refractions of sonic waves emitted by transmitted and reflected ripple shock fronts (or only transmitted shock front, in the case of reflected rarefaction). Both shock fronts are in causal interaction with the material interface, which plays the role of a rippled piston driving them both, and, through the piston, with each other. The sonic signals received and emitted by the interface become gradually weaker as the distance between the outgoing shock fronts and the interface increases, and this contribution to the RMI growth rates asymptotically approaches a constant value.
Due to the essential role of compressibility, this contribution cannot be captured by any of the prescriptions; one needs to solve the fully compressible perturbation problem. But when the EoS of shocked materials are not realistic, then the applicability of such analysis is still uncertain. Ideal-gas EoS, with $\gamma$s and the preshock pressure $p_0$ chosen as fitting parameters to match the one-dimensional (1-D) dynamics of observed shock flows, are often used for the compressible stability analysis of HEDP experiments with solid targets. For example, Ishizaki & Nishihara (Reference Ishizaki and Nishihara1997) selected $\gamma =3$ for shock-compressed plastic and a fitting preshock pressure of $p_0=0.7$ Mbar to theoretically describe the dynamics of the observed rippled-shock oscillations in the experiments by Endo et al. (Reference Endo1995). Similarly, Holmes et al. (Reference Holmes, Dimonte, Fryxell, Gittings, Grove, Schneider, Sharp, Velikovich, Weaver and Zhang1999) selected the values of $\gamma =1.45$ and $1.8$ for low-density foam and beryllium, respectively, with the preshock pressure $p_0=0.1$ Mbar roughly corresponded to the x-ray preheat for their analysis of the RMI experimental data by Dimonte & Remington (Reference Dimonte and Remington1993) and Dimonte et al. (Reference Dimonte, Frerking, Schneider and Remington1996). Both Ishizaki & Nishihara (Reference Ishizaki and Nishihara1997) and Holmes et al. (Reference Holmes, Dimonte, Fryxell, Gittings, Grove, Schneider, Sharp, Velikovich, Weaver and Zhang1999) successfully reproduced the results of their respective experiments during the small-amplitude stage of the perturbation development. However, since matching of the 1-D parameters does not necessarily translate into matching of the parameters determining the frequency and decay of the shock-fronts’ oscillations, which were observed by Endo et al. (Reference Endo1995) and affected the RMI growth measured by Dimonte & Remington (Reference Dimonte and Remington1993) and Dimonte et al. (Reference Dimonte, Frerking, Schneider and Remington1996), one cannot evaluate the accuracy of the compressible linear theory. This uncertainty is similar to that associated with the use of prescriptions. To eliminate it, the fully compressible linear RMI theory must be constructed for the realistic EoS, as envisioned by Richtmyer (Reference Richtmyer1960).
Numerical analysis of this problem was published by Ward & Pullin (Reference Ward and Pullin2011) for Mie–Grüneisen EoS; for details, see Ward (Reference Ward2011), and also by Tahir et al. (Reference Tahir, Shutov, Zharkov, Piriz and Stöhlker2011) for the case of aluminium and copper. Here, we report the analytic, fully compressible theory of the RMI for the reflected shock case with an arbitrary EoS, advancing the results of the earlier work by Wouchuk (Reference Wouchuk2001), Campos & Wouchuk (Reference Campos and Wouchuk2014) and Cobos-Campos & Wouchuk (Reference Cobos-Campos and Wouchuk2016). As examples of realistic EoS, we chose the venerable van der Waals (vdW) EoS for a non-ideal gas and the three-term, Grüneisen-type three-term EoS (Zel'dovich & Raizer Reference Zel'dovich and Raizer2002, Chapter XI, § 6), applicable to simple metals at pressures not exceeding $5$ Mbar. Of course, these results are readily modified for any EoS presented in analytical or tabular form. We demonstrate the specifics of the realistic-EoS, fully compressible RMI analysis, discuss various limiting cases and compare the theoretical predictions of the asymptotic growth rates with various prescriptions. The lack of Riemann problem similarity for different EoS (Quartapelle et al. Reference Quartapelle, Castelletti, Guardone and Quaranta2003) is addressed by finding the best fit for both EoS zero-order profiles. With this approach, the incident shock strength and the effective $\gamma$ values are selected to yield closely matching values for the target parameters: shock density and pressure jumps, postshock Atwood number and postshock dimensionless contact interface amplitude.
In this work, in addition to the extension of the formulation to be applicable to arbitrary EoS, we use a different resolution method compared with previous theoretical studies (Wouchuk Reference Wouchuk2001; Campos & Wouchuk Reference Campos and Wouchuk2014; Cobos-Campos & Wouchuk Reference Cobos-Campos and Wouchuk2016). Specifically, we address the initial-value problem by employing the direct inverse Laplace transform (ILT), which allows the final solution to be expressed as a closed integral. This approach overcomes issues related to the slow convergence of Fourier–Bessel coefficients and provides deeper insights into shock-to-shock acoustic coupling interactions, as demonstrated in Napieralski et al. (Reference Napieralski, Cobos, Sánchez-Sanz and Huete2024) for RMI developed in shock-flame interactions. Understanding the singularities in the complex plane is crucial for the ILT of the pressure field. For instance, the presence of non-removable singularities, as noted by Calvo-Rivera, Velikovich & Huete (Reference Calvo-Rivera, Velikovich and Huete2023), indicates the possibility of spontaneous acoustic emission (SAE) for one of the shocks. This allows the model to be extended to conditions where the interface is continuously excited by pressure disturbances.
The remainder of this paper is structured as follows. The problem formulation is presented in § 2 for the case of reflected and transmitted shocks. The following § 3 is devoted to the description of the pressure field comprised by the two postinteraction shocks. The analysis includes the possibility of SAE in the transmitted shock. Section 4 describes the velocity field split into acoustic and rotational perturbations. The asymptotic velocity field is also analytically derived. The RMI evolution of the interface is addressed in § 5 for different EoS and shock conditions. Results are compared with known prescribed formulae for the asymptotic growth rate and with $\gamma$-law EoS for metals. Lastly, conclusions are given in § 6. The resolution of the pressure field using hyperbolic coordinate transformation, followed by the Laplace transform and its inverse, is thoroughly detailed in Appendix A. The base-flow resolution particularised for vdW EoS, aluminium and copper is offered in the supplementary material available at https://doi.org/10.1017/jfm.2024.1010. Mathematica codes are also provided for both the base-flow properties and the perturbed RMI problem.
2. Problem formulation in Cartesian coordinates
2.1. Base-flow equilibrium condition for arbitrary EoS
The perturbation-free problem is readily given by the conservation equations across the shocks and the mechanical equilibrium condition at the interface. Because of there are five distinct flow regions, it is necessary to choose a clear notation criterion to facilitate their identification. In our case, as described in figure 1, numbers $0$, $1$ and $2$ are used to identify the preshock fluids, the flow behind the incident shock and the flow behind the transmitted and reflected shocks, respectively. An additional subindex is necessary to identify the position relative to the contact interface: $r$ and $t$ for the regions along where the reflected and transmitted shocks travel. The mechanical equilibrium applies at the interface. Before the interaction with the incident shock the equilibrium is dictated by the continuity of pressure $p_{0r}=p_{0t}$, provided that both fluids move at the same speed, i.e. null in the preshock materials’ reference frame.
In order to calculate the base flow equilibrium it is necessary to solve the conservation equations. In the absence of disturbances, the different regions separated by each shock wave and the contact interface are steady and uniform. Across each shock the flow variables changes are related to the Rankine–Hugoniot (RH) equations, namely
for mass, momentum and energy, respectively; where the square brackets $[\ldots ]$ denote the difference between the indicated quantities’ values ahead and behind the shock front. In (2.1a–c) $\rho$, $u$, $p$ and $E$ mean, respectively, density, velocity, pressure and internal energy evaluated in front and behind each shock wave, with the material velocity measured in the corresponding shock reference frame, so that $U=u+{\rm const}$., in relation to the material velocity indicated in figure 1.
The zero-order problem formulation is completed after providing the equilibrium conditions that must be imposed at the planar interface: $p_{2t}=p_{2r}$ and $U_2(x<0)=U_2(x>0)$. Furthermore, along with the isobaric condition, an additional parameter is needed to relate the two initial thermodynamic states of the fluids at both sides of the interface. Typically, the parameter employed is the preinteraction Atwood number
that is bounded $-1<\mathcal {A}_0<1$. The present study does not address the case when a rarefaction is reflected back after the incident shock refraction, which might manifest itself when $\mathcal {A}_0$ takes low (likely negative) values. However, the kind of the wave reflected (shock or rarefaction) cannot be predicted solely based on this parameter unless certain conditions are met, namely when the two materials share identical thermodynamic coefficients. So, only under these circumstances, a negative value for $\mathcal {A}_0$ ensures the occurrence of a reflected rarefaction. For an interface separating two ideal gases with similar adiabatic indices, the density jump $\rho _{0t}/\rho _{0r}$ is the only parameter that determines the problem configuration, rendering a shock-reflected problem whenever $\rho _{0t}/\rho _{0r}>1$. For the ideal-gas shock-interface interaction, the case of reflected rarefaction occurs if (5a)–(5b) in Yang et al. (Reference Yang, Zhang and Sharp1994) are satisfied. This requirement can be presented as a single inequality
as derived in Velikovich (Reference Velikovich1996), which is applicable to both $\gamma _r<\gamma _t$ and $\gamma _r>\gamma _t$ cases. In the present case, the boundary between the shock-reflected and rarefaction-reflected cases, typically referred as ‘total transmission’ since there are no waves reflected back, must be computed numerically for each specific EoS and flow conditions, as detailed in the supplementary material.
In dimensionless form, the pressure continuity at the contact interface, $p_{2t}=p_{2r}$, can be written as a function of the shock pressure jumps
where $\mathcal {P}_r=p_{2r}/p_{1r}$, $\mathcal {P}_i=p_{1r}/p_{0r}$ and $\mathcal {P}_t=p_{2t}/p_{0t}$ correspond to the pressure jump functions across the reflected, incident and transmitted shocks, respectively. On the other hand, the continuity in the streamwise velocity $U_2(x<0)=U_2(x>0)$ renders
that, as it can be seen, involves more parameters. The functions $\mathcal {M}_r=(D_{R}+U_1)/c_{1r}$, $\mathcal {M}_i=D_{I}/c_{0r}$ and $\mathcal {M}_t=D_{T}/c_{0t}$ stand for the reflected, incident and transmitted shock Mach numbers, respectively. The $c$ functions identify the corresponding sound speeds, that must be computed upon determination of the associated EoS. Equation (2.5) also depends on the corresponding density jump functions, which read, $\mathcal {R}_r=\rho _{2r}/\rho _{1r}$, $\mathcal {R}_i=\rho _{1r}/\rho _{0r}$ and $\mathcal {R}_t=\rho _{2t}/\rho _{0t}$. Additionally, the associated sonic coefficients, that in the general case are state functions, are defined as $\kappa = \rho c^2 / p$, where, again, $c$ represents the speed of sound. Finally, the post shock Mach numbers are defined as $\mathcal {M}_{2r}=(D_{R}+U_2)/c_{2r}$, $\mathcal {M}_{2i}=(D_{I}-U_1)/c_{1r}$ and $\mathcal {M}_{2t}=(D_{T}-U_2)/c_{2t}$.
Therefore, (2.4) and (2.5), along with the corresponding RH relations (2.1a–c), the EoS and the preshock parameters, can be solved for a given incident shock intensity. This process has been conducted in the supplemental material for the vdW EoS (including the specific case of the ideal gas) and for simple metals modelled with the three-term EoS (applied to aluminium and copper). Our analysis presumes the existence of a unique solution to this problem, although this is not always guaranteed. In most cases involving ‘regular’ EoS, such as those discussed in the supplementary material, a unique solution either exists for the transmitted and reflected shocks. The case that renders the reflected rarefaction wave being reflected is left for a separate analysis. However, since we are considering an arbitrary EoS, we must acknowledge the possibility of anomalous situations. For instance, an exponentially unstable shock scenario might arise, where, instead of a single transmitted shock front, two fronts are produced. Such situations are discussed in Menikoff & Plohr (Reference Menikoff and Plohr1989) and Kuznetsov (Reference Kuznetsov1989), and they are beyond the scope of this work.
2.2. Linearised perturbation equations
We consider a perturbed interface, separating two stagnant materials, that is weakly curved along the transverse coordinate $y$, before the incident shock first reaches it ($t\leqslant 0^{-}$), in the form of $\psi _i(t\leqslant 0^{-},y)=\psi _{i}^{0^-}\cos (ky)$, where the amplitude of the perturbation is $\psi _{i}^{0^-}$ and $k$ stands for the corresponding wavenumber. For linear theory to be applicable, the condition $\psi _{i}^0k\sim \varepsilon \ll 1$ must hold, with $\varepsilon$ being defined as the small parameter characterising the first-order perturbations magnitude. In what follows, all perturbation variables are of the same order $\varepsilon$. The dimensionless spatial coordinates are defined by $\bar {x}=kx$ and $\bar {y}=ky$, with the system of reference being fixed at the unperturbed interface, which moves with velocity $U_2 \,\hat {x}_0$ in the reference frame attached to the upstream stagnant gas.
When the incident shock wave has totally crossed the rippled interface, $t\geqslant 0^{+}$, two corrugated shock waves, the transmitted and the reflected one, are generated. The interfacial ripple amplitude, conveniently normalised with respect to its preshock value,
is initially compressed by the shock passage, thereby modifying its initial amplitude to $\xi _{i}^{0}=\xi _{i}(t=0^+)=1-U_2/D_{I}$, provided in (S.2.11) of the supplementary material. Note that the postshock amplitude differs from the preshock amplitude due to compressible effects of the shock. Since $\rho U ={\rm const}$, across the shocks, changes in density lead to corresponding adjustments in velocity. The corresponding shock front amplitudes are also scaled with $\psi _{i}^{0^-}$, thus, the transmitted and reflected shock ripples are equally reduced:
It is straightforward to see that, after incident shock refraction, the initial ripple amplitude of the transmitted and the reflected shocks can be determined by the relative distances between the most advanced parts of each front at the time the shock fully crosses the interface, thereby setting in their initial conditions through $\xi _{t}(t=0^+)=\xi _{t}^0=1-D_{T}/D_{I}$ and $\xi _{r}(t=0^+)=\xi _{r}^0=1+D_{R}/D_{I}$, where the relations $D_{T}/D_{I}$ and $D_{R}/D_{I}$ depend on the initial conditions and the corresponding EoS, as explicitly shown in (S.2.10) of the supplementary material.
All flow variables are disturbed behind the reflected and transmitted shocks, with the corresponding amplitude being of the same order as the perturbation parameter $\varepsilon$. This is used to define the dimensionless order-of-unity functions (bar symbols) for the pressure, density, longitudinal velocity and transverse velocity variables, namely
where subscript ‘$j$’ indicates whether we are identifying the perturbations behind the transmitted shock (denoted as ‘$j=t$’) or the reflected shock (denoted as ‘$j=r$’). Note that definitions in (2.8) utilise the periodic symmetry of the flow within the linear regime. Note that the spatial coordinates and velocities are measured in the reference frame of the interface, which moves together with the compressed material for $\tau _j>0$, where $\tau _{j}$ is dimensionless time constructed with the aid of the local speed of sound: $\tau _{j}=k c_{2j} t$. This dimensionless transformation proves to be advantageous when expressing the linearised Euler equations since they become devoid of parameters. As usual, the linear Euler equations have been obtained by retaining the first-order perturbations $O(\varepsilon )$ and neglecting those of higher order on either side of the interface, see
for the mass, streamwise momentum, transverse momentum and energy conservation equations, respectively, with the latter being reduced to the isentropic condition. The Euler equations can be rearranged to give the sound wave equation which, in this particular case, takes the form of the Klein–Gordon equation,
The noticeable advantage of working with two free-parameter sound wave equations comes with the price of handling two temporal scales, which are interconnected by the factor
This parameter plays a crucial role in the acoustic coupling between the reflected and transmitted shocks, as will be discussed in detail later.
The sound wave equations at both sides of the interface must be integrated, upon providing the corresponding initial and boundary conditions. The initial conditions are given by the initial shock ripple amplitudes, $\xi _t^0$ and $\xi _r^0$, provided in (2.7a,b) along with the isobaric condition ${\bar {p}}_{j}(\tau _j=0^+)=0$. The boundary conditions, on the other hand, are dictated by the perturbed RH equations and the equilibrium condition at the interface. For the reflected shock, at $\bar {x}_{rs}=-\mathcal {M}_{2r}\tau _r$, and the transmitted shock, at $\bar {x}_{ts}=\mathcal {M}_{2t}\tau _t$, respectively, they read as
where the subscript ‘$s$’ is used here to indicate that perturbations are evaluated at their respective reflected and transmitted shock position. The functions $\mathcal {M}_{2j}$, $\mathcal {R}_{j}$ and $h_{j}$, refer to the associated postshock Mach number, density jump and the D'yakov–Kontorovich (DK) parameter, named after the pioneering works of D'yakov (Reference D'yakov1954) and Kontorovich (Reference Kontorovich1957). For the sake of generality, they are included here as independent parameters, yet they are ultimately determined by the initial conditions, the EoS of the fluids and incident shock intensity, as demonstrated in the supplementary material for vdW EoS and three-term EoS for condensed materials.
The RH conditions (2.12) can be further reduced by taking the total derivative of the streamwise velocities $\bar {u}_j$ along the corresponding shock trajectories, together with the linear Euler equations (2.9). They read as
for the transmitted and reflected shocks, respectively. Equations (2.13a) and (2.13b), complemented with (2.12a) for the shock ripples evolution, form a closed system of equations for the perturbed shocks.
The problem formulation is completed after providing the mechanical equilibrium condition at the contact interface: pressure and velocity are continuous across the interface. In terms of pressure, they read as
where $\vartheta =\kappa _{2r}/\kappa _{2t}$ is the ratio of the sonic coefficients at both sides of the interface. For ideal gas EoS, this ratio reduces to $\gamma _r/\gamma _t={\rm const}.$, but this function depends on the corresponding thermodynamic state of the material, and therefore on the shock intensity, for more complex EoS, as provided in the supplementary material.
This problem can be readily numerically integrated using the method of characteristics, as done in Calvo-Rivera et al. (Reference Calvo-Rivera, Velikovich and Huete2023) for piston-driven shocks, yielding the spatial and temporal pressure distribution within the region bounded by the reflected and transmitted shocks. Nevertheless, the theoretical developments in the following sections, based on the formulation presented in Wouchuk (Reference Wouchuk2001), will provide analytical expressions for key quantities, including the asymptotic growth rate of the RMI.
3. The pressure field
3.1. The resolution of the pressure field
This section examines the physical properties of the pressure field, drawing on the mathematical properties of the corresponding Laplace transform detailed in Appendix A. Therefore, although the discussion primarily addresses physical concepts, it is necessary to explicitly reference some of the properties of the associated Laplace transform functions. In particular, the evolution of the pressure field is derived by applying the inverse transformation to the Laplace transform functions described in Appendix A.2).
In the transmitted-shock and reflected shock regions, it can be computed with the aid of
where $\mathbb {P}_{j}$ is the pressure Laplace transform with respect to the coordinate $\eta _j$, $s_j$ is the Laplace variable, $\mathrm {i}=\sqrt {-1}$ is the imaginary unit and $r$ is a real number to the right of the singularities of $\mathbb {P}_{j}(s_j)$. The independent variables $\eta _j$ and $\chi _j$ refer to the hyperbolic coordinates
defined conveniently to transform the triangular integration domain into a rectangular one (Zaidel’ Reference Zaidel’1960; Briscoe & Kovitz Reference Briscoe and Kovitz1968).
Keeping in mind the analysis of the singularities within the complex plane detailed in Appendix A.4, there are two contributions to the Bromwich integral, the first one is the integral around the integration contour in which we introduce the branch cut to avoid the multivalued nature of the complex functions (A10a) and (A10b), in particular from the line integrals along the imaginary axis, specifically between the branch points. The second contribution comes from the existence of imaginary poles inside the contour, which are calculated by the application of the Cauchy residue theorem. To clarify the analysis, we split the contributions in the form
where the first term of the right-hand side, obtained by the residue theorem, is associated with the non-decaying oscillations of a shock front producing SAE. The second term, which requires the evaluation of the integral along each path of the contour excluding the vertical line that pass through $r$, as sketched in Appendix A.4, is associated with the damped oscillations.
Beginning with the latter, it should be noted that the value of the ratio of the speeds of sound, $\beta$, determines the branch cut on each side. When $\beta >1$, the influence of the reflected shock oscillations involves a frequency that is higher than that of the transmitted shock oscillations. Therefore, the branch cut must be extended to range $\{-\mathrm {i}\beta,\mathrm {i}\beta \}$ for the transmitted side but not for the reflected one, see Appendix A.4. When $\beta <1$, the opposite happens, the oscillations of the transmitted shock exhibit a higher frequency than the reflected shock, signifying into an extension of the branch cut of the reflected side with the range $\{-\mathrm {i}\beta ^{-1},\mathrm {i}\beta ^{-1}\}$ due to the change in the dominant evanescent acoustic waves. It results in the closed integral expression for the pressure field which read as
where $\mathbb {P}_{j}$ is given by the solution of the functional system corresponding to each zone. The domain of validity of (3.4) goes $0\leqslant \chi _t \leqslant \tanh ^{-1}\mathcal {M}_{2t}$, for the transmitted side, and $-\tanh ^{-1}\mathcal {M}_{2r} \leqslant \chi _r \leqslant 0$ for the shock-reflected region within the domain $0\leqslant \chi _t \leqslant \tanh ^{-1}\mathcal {M}_{2t}$. Also $\zeta _j$ depends on the value of $\beta$, being
for the transmitted and reflected cases, respectively. The superscript ${dec}$ is used in (3.4) because the branch cut contribution always decays, similar to Bessel functions. These terms represent the typical transient growth contribution for the RMI and are consistently reported in the literature.
To obtain the contribution from the residues, it is necessary to consider the imaginary poles that arise in the complex plane, which are the isolated singularities of (A21a) and (A21b) for the transmitted and reflected cases, respectively. Such singularities are the complex roots of $D(s_{j})=0$ in (A22). The amplitude of the non-decaying oscillations is readily given by the residue theorem
where $s_{I,jk}$ is the $k$th imaginary pole in (A24).
We consider the case of the shock natural frequency, i.e. there are no secondary modes associated with reverberations with the interface and/or the reflected shock. In this case,
with use made of the symmetry property
Using the same decomposition of contributions as in (A21), the shock pressure amplitude associated with non-decaying oscillations can be expressed as
3.1.1. Pressure perturbations at the shock fronts
The evolution of the pressure field at the shocks is readily determined by the superposition of the decaying and non-decaying oscillating contributions, namely
when that SAE condition is met. Otherwise, only the decaying term contributes.
To illustrate the different possibility that may arise depending on the shock properties and the EoS parameters, figure 2 illustrates the pressure perturbation history at the shocks for four distinct cases. The transmitted shock coordinate corresponds to the positions $\bar {x}_{ts}=\mathcal {M}_{2t} \tau _t >0$ (right half-space) and the reflected shock to $\bar {x}_{rs}=-\mathcal {M}_{2r} \tau _r <0$ (left half-space). In figure 2(a), the pressure field is computed for the same case as in Yang et al. (Reference Yang, Zhang and Sharp1994), representing a shock wave propagating from air to ${\rm SF}_{6}$ with an incident Mach number $\mathcal {M}_{i}=1.24$. Figure 2(b) depicts a scenario with similar gas conditions but a stronger incident shock, $\mathcal {M}_{i}=3$. Here, the pressure field of the transmitted shock changes its initial slope sign and exhibits a sudden rise that subsides quickly. Figure 2(c) shows computations where SAE occurs at the transmitted shock side. The asymptotic non-decaying contribution is indicated with a dashed line. This case involves a highly compressible vdW gas with significant covolume and Coulomb contributions on the heavy gas side, similar to those described in Bates & Montgomery (Reference Bates and Montgomery1999), Huete et al. (Reference Huete, Velikovich, Martínez-Ruiz and Calvo-Rivera2021) and Calvo-Rivera et al. (Reference Calvo-Rivera, Velikovich and Huete2023). The conditions are those marked with a triangle in figure S2 of the supplementary material. Figure 2(d) presents a computation where the parameter $\beta$ is less than unity, necessitating a different contour path in its resolution via the ILT. For this case, it is noticeable the change in the amplitude of the oscillations of the reflected shock. Finally, figure 2(e,f) illustrate the results for a shock wave propagating through aluminium and impacting copper, leading to the reflection of a backward-moving shock wave. The incident Mach numbers are $\mathcal {M}_{i}=1.3$ and $\mathcal {M}_{i}=2.5$, corresponding approximately to the weak and strong test cases analysed in Tahir et al. (Reference Tahir, Shutov, Zharkov, Piriz and Stöhlker2011), with incident shock pressures of approximately $p_1= 0.28$ Mbar and $2.9$ Mbar, respectively. Base-flow conditions are described in figure S3 of the supplementary material.
The six panels in figure 2 demonstrate the robustness of the ILT method compared with the previously used separation of variables technique, which describes the pressure field in terms of Fourier–Bessel series. Unlike the separation of variables method, the use of which is limited by the requirement of convergence of its coefficients, the ILT method consistently provides accurate solutions even under varying conditions and long-time regimes, this is due to the boundary conditions being time dependent which makes it necessary to obtain a number of coefficients that tends to infinity for larger times, as seen in the first chapter of Duffy (Reference Duffy2004).
3.1.2. Influence of the pressure disturbances on the interface
Once we understand the evolution of the pressure field at the shock fronts and the potential for SAE, it is natural to consider the impact on the interface. Ultimately, this affects the development of the RMI, which is quantified in § 5. Starting from the decaying contribution, a direct variable change transformation in (3.4) yields the pressure field in Cartesian coordinates,
which is applicable within the domain enclosed by the corresponding shock and the interface.
To better understand the non-decaying oscillating contribution, we can isolate this effect by analysing a scenario where a shock oscillates with constant frequency and amplitude. First, we refer again to the dispersion relation (A22). From the analysis of the singularities in the complex plane performed before, purely imaginary poles arise when $\sigma _{cj}>\sigma _{bj}$, which ultimately translates into oscillations with frequency $\smash {\omega _{js}>(1-\mathcal {M}_{2j}^{2})^{1/2}}$. Using a normal mode decomposition for the acoustic frequency and wavenumber, referred to from now on as the acoustic eigenvalues, $\omega _{j}$ and ${k}_{j}$, the shock oscillation frequency is related to the acoustic eigenvalue through the sound wave equation $\smash {\omega _{j}^{2}={k}_{j}^{2}+1}$ and the compatibility condition $\smash {\omega _{js}=\omega _{j}-\mathcal {M}_{2j}{k}_{j}}$, to give
As the wrinkled shock oscillates, it generates pressure disturbances in the downstream fluid in the form of sound waves (see § 90 of Chapter IX in Landau & Lifshitz (Reference Landau and Lifshitz1987)). A simple inspection of ${k}_{j}$ in (3.12a,b) reveals that, depending on the shock oscillation frequency, these sound waves may or may not directly impact the contact interface. Therefore, three distinct regimes can be identified, as follows.
(i) When $\smash {\omega _{ts}^2 < 1-\mathcal {M}_{2j}^{2}}$, pressure disturbances at the shock decay according to the power-law $\smash {t^{-3/2}}$ and diminish exponentially with distance from the shock, resulting in evanescent waves in the longitudinal direction. This scenario corresponds to the typical case studied in RMI for both transmitted and reflected shock fronts moving in ideal gases. Therefore, the acoustic influence only affects the interface at early times, as the exponential decay quickly reduces its relative contribution.
(ii) When $\smash {1-\mathcal {M}_{2j}^{2}\leqslant \omega _{ts}^2<1}$, the stability behaviour of the shock front switches to a neutrally stable regime, and non-decaying acoustic waves are radiated by the shock. This regime, referred as SAE, is a manifestation of the DK instability (Calvo-Rivera et al. Reference Calvo-Rivera, Velikovich and Huete2023). In such case, the acoustic eigenvalues become real, with the wavevector pointing to the radiating shock direction. In the theory of dynamical systems, this corresponds to a transition through a Hopf-like bifurcation, see Clavin & Williams (Reference Clavin and Williams2009, Reference Clavin and Williams2012), to a limit cycle which will exhibit permanent oscillations, as explained in pp. 252–257 of Strogatz (Reference Strogatz2015). However, regarding the interface behaviour, the contact interface will not be affected by the acoustic waves, as they cannot reach the hydrodynamic interface. The permanently oscillating waves travel to the right in the compressed gas (and interface) reference frame but not reaching the shock front or the interface, as sketched in figure 3 for the transmitted shock case.
(iii) When $\omega _{js} \geqslant 1$, the behaviour of the corresponding shock front dynamics remains the same, but the acoustic influence changes as the longitudinal projection of the wavenumber in the compressed material reference frame points backwards. Consequently, sound waves propagate to the left, reaching the interface, as sketched in figure 4 for the transmitted shock.
Further scenarios can also be explored as $\omega _{js}$ increases beyond unity: the radiated sound waves can reach the reflected shock, causing it to oscillate with constant amplitude and frequency, and/or they can be reflected back towards the transmitted shock, imposing a second mode of oscillation, similar to what occurs with piston-driven shocks (Calvo-Rivera et al. Reference Calvo-Rivera, Velikovich and Huete2023). However, these cases are rare unless the shock oscillation frequency is externally imposed, as described in Velikovich et al. (Reference Velikovich, Wouchuk, Huete Ruiz de Lira, Metzler, Zalesak and Schmitt2007), where the shock travels across a non-uniform medium.
To facilitate the discussion, we consider the specific case where SAE occurs only in the transmitted shock, i.e. $\sigma _{bt}<\sigma _{ct}$, as depicted in figures 3 and 4 for cases (ii) and (iii), respectively. As per the former, the pressure field associated with the sound radiation by the transmitted shock, which also corresponds to the pressure field in the long-time regime in this case, can be expressed as a piecewise function of $\bar {x}/\tau _{t}$, namely
where the non-zero pressure field corresponds to the orange-shadowed region in figure 3. Notice the decaying contribution given in (3.11) is present within the whole domain. The two limiting cases that bound case (ii) are $\omega _{ts}\rightarrow (1-\mathcal {M}_{2j}^{2})^{1/2}$ and $\omega _{ts}\rightarrow 1$. Per the former, the value of $k_{t}/\omega _{t}\rightarrow \mathcal {M}_{2t}$, which corresponds to the trajectory of the transmitted shock wave, thereby eliminating the orange-shadowed region. Regarding the latter, the value of $k_{t}/\omega _{t}\rightarrow 0$, which corresponds to the trajectory of the entropy-vorticity waves (vertical lines in figure 3). Such a case renders a standing acoustic wave and the domain of influence for the orange-shadowed region is the space between the transmitted shock and the interface.
For the backwards running waves to reach the contact interface (case (iii)), the acoustic wave number ${k}_t$ must be negative, implying that $\omega _{ts} > 1$, which is a more stringent condition to satisfy if the shock oscillation frequency is not externally imposed. Nonetheless, it is worthwhile to outline the consequences if this condition is met, which can be readily anticipated by simple inspection of figure 4. It can be seen that four distinguished regions are identified, since the train of acoustic waves radiated by the transmitted shock are partially reflected from and transmitted through the contact interface, as detailed below. The following derivation is essentially the same as the one found in § 66 of Chapter 8 in Landau & Lifshitz (Reference Landau and Lifshitz1987).
As in (3.13), the hyperbolic character of the problem allows as to write the pressure field as a piecewise function of $\bar {x}/\tau _{t}$, namely
where the first region corresponds to the orange-shadowed zone in figure 4 and that is similar to that in figure 3, with the coefficient ${\rm \pi} _t$ being obtained from the residue theorem described before. The following green-shadowed region is composed of the coexisting shock-radiated waves moving backwards and the interface-reflected waves moving forwards, with the amplitude of the latter being readily given by the reflection coefficient for the dimensionless pressure amplitude, $\mathfrak {R}$. And the blue-shadowed region corresponds to the pressure field that crosses the interface with an amplitude that is determined by the transmission coefficient, $\mathfrak {T}$. Such are obtained by applying the boundary conditions at the material surface, (2.14), taking the form of
that involve the wavenumber at the transmitted and reflected sides, respectively, which are attached to the compatibility condition, i.e. both of them must fulfil the sound wave equation being subjected to the same frequency, $\omega _{t}$ giving ${k}_{r}^2=\omega _{t}^{2}/\beta ^{2}-1$ and ${k}_{t}^{2}=\omega _{t}^{2}-1$, and it can be checked that they have the same form of those given in problem 2 of § 84, Chapter IX of Landau & Lifshitz (Reference Landau and Lifshitz1987). Notice that, as ${k}_{t}={k}_{t}(\omega _{t})$ and ${k}_{r}={k}_{r}(\omega _{t})$, the acoustic impedance of the contact interface depends on the acoustic frequency because it changes the acoustic wavenumber, as expected. For instance, in the limiting case where $\omega _{ts}\rightarrow \infty$, the value of ${k}_{t}/\omega _{t}\rightarrow -1$ corresponds to the wavenumber of 1-D sound waves (represented by the grey lines). These scenarios extend beyond the current analysis because, with the modified wavenumber, the constant-amplitude sonic waves that are reflected and transmitted will reach both the transmitted and reflected shocks, respectively. This interaction induces at least secondary modes in the transmitted shock, as described by Calvo-Rivera et al. (Reference Calvo-Rivera, Velikovich and Huete2023) for piston-driven shocks. The present analysis remains consistent because, as $\beta \rightarrow 1$ and $\vartheta \rightarrow 1$, the coefficients $\mathfrak {R}\rightarrow 0$ and $\mathfrak {T}\rightarrow 1$, corresponding to the case of a weak discontinuity at the interface, addressing the problems already solved in Wouchuk, Huete Ruiz de Lira & Velikovich (Reference Wouchuk, Huete Ruiz de Lira and Velikovich2009) and Huete Ruiz de Lira, Velikovich & Wouchuk (Reference Huete Ruiz de Lira, Velikovich and Wouchuk2011).
The opposite limit, involving the total reflection of a sonic wave from a perfectly solid wall, cannot be directly addressed by simply taking limits of the parameters introduced here, i.e. $\mathfrak {T}\neq 0$. This is because the nature of the refracted sonic wave fundamentally changes. As explained by Haberman (Reference Haberman2013), see Chapter 4, § 6, the refraction–reflection problem is characterised by the angle of incidence of the sonic wave defined by the projection of the wavevector in the transverse direction. In our formulation, this angle is just $\sin ^{-1}(1/\omega _{t})$. When this angle is greater than the critical value, i.e. when $\sin ^{-1}(1/\omega _{t})>\sin ^{-1}(1/\beta )$, the refracted (transmitted acoustic wave in the reflected-shock region) sonic wave is no longer a running wave but it becomes evanescent. Instead, an evanescent wave is transmitted, which decays exponentially with the distance from the interface, as readily seen by the compatibility condition ${k}_{r}^2=\omega _{t}^{2}/\beta ^{2}-1$, since ${k}_{r}$ turns complex when $\omega _t<\beta$. By introducing this consideration and applying the same boundary conditions for the pressure and its gradient, the wavenumber for the refracted wave becomes complex. Consequently, one finds that $\mathfrak {R}=-1$, similar to the piston-driven shocks discussed in Calvo-Rivera et al. (Reference Calvo-Rivera, Velikovich and Huete2023).
For the sake of completeness, we briefly comment on the possibility of having permanent acoustic waves at the interface from the Laplace transform equations. With this purpose, it has proved convenient to write the pressure field in the transmitted shock side as functions of the coordinate $\chi _t$ and variable $q_t$,
where the equation particularised at the shock (A21a) is recovered if $\chi _{t}=\chi _{ts}$. This equation is particularly useful because it can easily adapt to different Riemann-type configurations. For instance, if we assume that no disturbance can propagate back to the transmitted shock from behind, making, $F_{t}^{+}={\rm const}.$, the expression simplifies to match (48) of Wouchuk et al. (Reference Wouchuk, Huete Ruiz de Lira and Velikovich2009) for isolated shocks. Despite the evident differences, a similar analysis can be applied. The first term in (3.16) introduces poles corresponding to the secondary DK modes. Conversely, the second term in (3.16) contains a common denominator, whose zeros are the same as those previously derived from the dispersion relation, but shifted relative to the shock reference frame,
where $\tilde {s}_{t}=\sinh {(q_{t}+\chi _{ts}-\chi _{t})}$. The poles given for each $\chi _{t}$ value, associated with different relative speeds, can be expressed in the form of
which, through manipulation and substitution of the hyperbolic change of variables, allows us to obtain the pressure field in the spatiotemporal domain,
where
as a proof of consistency between the Laplace transform and the normal mode analyses. An alternative scenario arises when the reflected shock satisfies $\omega _{rs} > 1$, with the acoustic radiation hitting on the interface and travels to the right from the reflected shock. Since the physical discussion is similar, this case will not be specifically considered in this work.
3.2. The infinitely dense limit for the heavy material: the sequential model
As commented on in Velikovich et al. (Reference Velikovich, Schmitt, Zulick, Aglitskiy, Karasik, Obenschain, Wouchuk and Cobos Campos2020), when the Atwood number is sufficiently high $1-\mathcal {A}_0\ll 1$, the RMI problem can be reduced, to some degree, to a shock-wall reflection. However, it is convenient to take this limit cautiously because what really distinguishes the interaction is the postshock Atwood number, in our formulation it is $\mathcal {A}_2= (\beta ^2-\vartheta )/ (\beta ^2+\vartheta )$. The condition $1-\mathcal {A}_2\ll 1$ may not be easily anticipated by preshock conditions for complex EoS, so the infinitely dense limit is conveniently rewritten in terms of $\beta ^2/\vartheta \gg 1$. In such a case, the mathematical description of the pressure field can be significantly simplified. In terms of the mathematical formulation, this distinguished limit allows us to obtain analytic expressions for the pressure fields at both sides. When $\beta ^2/\vartheta \gg 1$, expression (A4a–c) for the pressure gradient becomes $\bar {l}_{ri}\simeq 0$, that is equivalent to that with the problem of shock wave reflected by a weakly perturbed wall, as studied in Briscoe & Kovitz (Reference Briscoe and Kovitz1968) and extended to arbitrary EoS in Calvo-Rivera et al. (Reference Calvo-Rivera, Velikovich and Huete2023).
The physics behind the decoupling of the perturbation evolution of reflected and transmitted shocks can be understood as follows if dimensional variables are considered. Due to the significantly higher speed of sound on the reflected-shock side, the corrugated contact surface evolves more slowly compared with the dynamics of the reflected shock, whose oscillations decay before the interface shape undergoes significant changes. As a result, on the time scale of the interface shape variation, the reflected-shock side maintains a constant pressure from the low-density left-hand side of the interface. For the high-density right-hand side, the situation corresponds to a rippled surface of a compressible fluid is instantly loaded with a constant, uniform pressure, resulting in an RMI-like instability, as first analysed by Nikolaev (Reference Nikolaev1965) for the opposite limit $\beta ^2/\vartheta \ll 1$. In our dimensionless formulation, the reflected shock operates independently of any knowledge about the transmitted shock, while the latter is completely influenced by the former. Thus, the resolution occurs in a sequential manner: initially addressing the reflected shock, such as in the piston-reflection problem, and subsequently utilising the pressure field at the contact interface to resolve the transmitted shock. We refer to this as the sequential model. The mathematical details of the function decomposition in the Laplace transform are provided in Appendix A.5. By applying the same procedure used in the previously described ILT, the pressure field in the spatiotemporal domain can be computed.
For the sake of simplicity, the sequential model has been presented and computed in the absence of SAE, since the main purpose is the evaluation of the accuracy in regular conditions. For example, in figure 5(a), the pressure field is computed for the same case as in figure 2(a) corresponding to Yang et al. (Reference Yang, Zhang and Sharp1994). Figure 5(b) is a similar case but with higher Atwood number and shock strength. The orange line corresponds to the fully consistent model and the black-dashed line with the sequential model. It is observed that, even when the parameter $\beta ^{2}/\vartheta$ is not very large, the agreement in the shock pressure perturbations is fairly good. However, it should be noted that a relatively good agreement in the pressure evolution of the shocks does not necessarily guarantee a similar match in the evolution of the contact interface, since the simplification stems from a singular perturbation in neglecting the pressure gradient at the interface. Therefore, the validity of this method should be checked by comparing the growth rate given from the full formulation against the present limit, as done in § 5.
4. The acoustic, rotational and asymptotic velocity fields
This section is dedicated to the analysis of the velocity field. As previously established in Landau & Lifshitz (Reference Landau and Lifshitz1987) and Clavin & Williams (Reference Clavin and Williams2009, Reference Clavin and Williams2012), the velocity field can be decomposed into two distinct types of waves: an acoustic wave, which is generated by the pressure sound waves emitted by each front, and an entropy-vorticity wave, which, under linear stability analysis, represents an isobaric and incompressible flow convected by the unperturbed flow velocity. The general equation for such velocity field can be given by the Helmholtz decomposition in the form of
where $\boldsymbol {\nabla }\times (\boldsymbol {\nabla }\times {\boldsymbol {\bar {v}}}_{j})= \boldsymbol {\nabla }\times {\boldsymbol {\bar {w}}}_{j}$ is the curl of the vorticity field $\boldsymbol {\bar {w}}_j$. It is convenient to use Helmholtz decomposition of the velocity field by the superposition of its acoustic or potential part and its rotational or entropic contribution, namely ${\boldsymbol {\bar {v}}}_j(\boldsymbol {{\bar {x}}},\tau _j)={\boldsymbol {\bar {v}}}_j^{ac}(\boldsymbol {{\bar {x}}},\tau _j)+{\boldsymbol {\bar {v}}}_j^{rot}(\boldsymbol {{\bar {x}}})$. Note that the rotational contribution remains constant for each fluid particle, as derived from the linear Euler equations. Therefore, any changes in the rotational field (or vorticity) must be externally induced, as exemplified here by the corrugated oscillating shocks.
Similar to the pressure field that has already been calculated in the previous section, the velocity field can be expressed as a combination of the decaying and non-decaying components of that pressure field. First, we introduce the acoustic velocity field with contributions from the decaying part, followed by the permanent oscillatory component. This results in a piecewise function for the velocity field, depending on the domain of SAE influence.
4.1. The acoustic wave
The acoustic component of the velocity field can be calculated directly from Euler momentum equations (2.9). With use made of the hyperbolic transformation (3.2a,b), and the corresponding Laplace transforms, the streamwise component of the acoustic velocity field in the transmitted-shock and reflected-shock sides reads as
where ${\bar {v}}_{ji}^0$ is suited for both sides, namely, ${\bar {v}}_{ti}^0=-\mathcal {M}_{2t}(\mathcal {R}_{t}-1)\xi _t^0$ and ${\bar {v}}_{ri}^0=\mathcal {M}_{2r}(\mathcal {R}_{r}-1) \xi _r^0$ are the initial value for the transverse velocity at the interface, which are proportional to the initial shock corrugation amplitudes, defined after (2.7a,b) and explicitly given in (S.2.10) of the supplementary material. The spatiotemporal evolution can be obtained with the ILT, similarly as done in (3.3) for the pressure field. For the streamwise components, we obtain
for the transmitted and reflected sides, respectively. On the other hand, for the transverse velocity contribution, we find
In RMI, an important particular case involves the evolution of velocity profiles at the interface. This simplifies the aforementioned equations, as for $\bar {x}=\eta _{j}\sinh {0}=0$ and for the temporal variable $\tau _{j}=\eta _{j}\cosh {0}=\eta _{j}$, yielding
for the streamwise and transverse velocity contributions, respectively.
As previously mentioned, the expressions presented above are only valid when the shock front oscillates with a frequency lower than $(1-\mathcal {M}_{2j}^{2})^{1/2}$, which was referred as case (i). When the frequency of the shock exceeds this limit, permanent sound pressure waves are emitted in the compressed fluid, as identified previously as cases (ii) and (iii).
For clarity, we will restrict our consideration to $\omega _{st}>1$, associated with cases where the transmitted shock front can enter the SAE regime and these waves reach the interface, i.e. case (iii). In this scenario, the longitudinal velocity field can be expressed as $\bar {u}_{j}^{ac}(\tau _{j},{\bar {x}}) =\bar {u}_{j}^{dec}(\tau _{j}, {\bar {x}}) + \bar {u}_{j}^{osc}(\tau _{j},{\bar {x}})$, where the decaying component is given in (4.3) with the change of variables to Cartesian coordinates. The non-decaying component is straightforward to calculate, as it only involves differentiating and integrating corresponding harmonic functions in (3.14). Consequently, the transmitted- and reflected-side normal velocity components can be expressed as
where $Q_{t}={\rm \pi} _{t}{k}_{t}/\omega _{t}$ is the amplitude of the non-decaying contribution. The transverse part is equally simple to obtain and will be omitted here for the sake of conciseness. Upon close examination of these expressions, the influence of the SAE at the interface becomes clear. When the transmitted front oscillates at a frequency higher than one, the wavenumber turns negative. This causes the non-decaying component of the longitudinal interface velocity to oscillate with a mean value around zero, while the decaying part increases from zero to a final constant value. As a result, the oscillatory and decaying components compete over time without damping the oscillations.
4.2. The vortex-entropic wave
Due to the steady state of the rotational field, the governing equation (4.1) can be expressed as $\nabla ^{2}\boldsymbol {{\bar {v}}}_{j}^{rot}=-\boldsymbol {\nabla }\times \boldsymbol {\bar {w}}_{j}$. This vector equation can be decomposed into its streamwise and transverse components,
Recognising the incompressibility of the linear rotational velocity field, the latter equation can be replaced with $({\rm d}{\bar {u}^{rot}_j})/({\rm d}\kern0.7pt {\bar {x}})=-{\bar {v}}_{j}^{rot}$, simplifying the differential to be solved to only (4.7a) for each side of the contact interface. This fact allows us to calculate the value of ${\bar {w}}_{j}$ using RH conditions and Euler equations, namely
where the positive sign must be chosen for the transmitted shock $j=t$ and the negative sign for the reflected shock $j=r$. Notice that $\tau _t=\bar {x}/\sinh {\chi _{ts}}$ and $\tau _r=\bar {x}/\sinh {\chi _{rs}}$ represent the time when the transmitted and reflected shock fronts, respectively, have just arrived at the evaluated position $\bar {x}$, where $\sinh {\chi _{js}}=\mathcal {M}_{2j}/(1-\mathcal {M}_{2j}^2)^{1/2}$. The amplitude of the vorticity disturbances relative to the shock pressure perturbations is
Equation (4.7a) is a second-order linear differential equation, the solution of which involves a combination of the homogeneous and particular solutions. The particular solution can be determined using the method of undetermined coefficients, employing the integral equation as a proposed solution, as outlined in Wouchuk & Cobos-Campos (Reference Wouchuk and Cobos-Campos2018).
The homogeneous part can be easily solved by using an integrating factor of the form $A {\rm e}^{-\bar {x}}+B {\rm e}^{ \bar {x}}$ and neglecting the non-regular contributions of the general solution, in the case of the transmitted zone $B=0$ while for the reflected one $A=0$. In determining the integration constant, a boundary condition arises from the observation that the only contribution to the vorticity field stems from the oscillations of the shock. Hence, it becomes evident that at $\bar {x}=0$, the rotational contribution to the velocity field is zero. Thus, the ultimate expression for the normal component of the rotational velocity corresponds to the particular solution of (4.7a), which takes the form of
Given the incompressibility of the linear rotational field, the transverse contribution can be determined as the derivative with respect to ${\bar {x}}$ of (4.10), that renders
Notice that, as before the expressions are only valid when the frequency of the shock is lower than $\omega _{ts}<(1-\mathcal {M}^{2}_{2t})^{1/2}$. If one of the shock front exceeds that limit, the shock front enters in the SAE regime. In order to see what modification implicates such a change we consider the case presented in figure 2(c), i.e. the transmitted shock wave exhibits permanent oscillations. In terms of mathematical treatment, this only implies that the transmitted rotational velocity field should contains another particular solution proportional to $\sin {({k}^{rot}_{t} \bar {x})}$, where the rotational wavenumber is given by ${k}^{rot}_{t}=s_{I}(1-\mathcal {M}_{2t}^{2})/\mathcal {M}_{2t}$. The corresponding constant is obtained from the method of undetermined coefficients. Thus, it is possible to include the permanent oscillatory part of the rotational velocity field as
for the streamwise and lateral velocities associated with the transmitted-shock side, respectively. It is readily seen that SAE-induced rotational disturbances do not contribute to the velocity field at the interface $\bar {x}=0$, and therefore, they have no effect on the growth rate.
4.3. Asymptotic velocity field
The asymptotic velocity field depends on the shock oscillation frequency that ultimately determines the cases (i), (ii) and (iii) commented before. As per the former, as the shocks gradually attain planarity over extended periods, the acoustic contribution diminishes, resulting in an incompressible velocity field for $\tau _{j}\gg 1$. Consequently, (4.1) simplifies to
which is similar to that for the rotational field, albeit with the additional consideration of accounting for the history of the acoustic perturbation influence, which manifests as the homogeneous component of the solution to the differential equations. Consequently, the system of differential equations that must be solved for each side is written as follows:
The general solution for the streamwise and transverse asymptotic velocity field in the transmitted side can be written as
where the rotational contribution corresponds to the particular solution. The homogeneous part, involving the exponentially decaying term, can be calculated via the final value theorem,
Notice that this term, commonly referred to as the asymptotic growth rate in the RMI context, is pivotal and will be thoroughly examined in § 5.
Similarly, the reflected field can be articulated as
where the homogeneous term involves ${\bar {u}}_{ri}^{\infty }={\bar {u}}_{ti}^{\infty }/\beta$.
As discussed in Wouchuk (Reference Wouchuk2001) and Cobos-Campos & Wouchuk (Reference Cobos-Campos and Wouchuk2016), an alternative method exists for determining the asymptotic velocity, particularly the value of ${\bar {u}}_{ji}^{\infty }$. This involves substituting the independent variable ${\bar {x}}$ in (4.14a) with the shock position, denoted as ${\bar {x}}=\eta _{j}\sinh {\chi _{js}}$. Upon transformation in the Laplace variable, an expression in the Laplace variable for the asymptotic normal velocity can be found, namely
for both the transmitted and the reflected sides. By focusing on the second term in the aforementioned expression, it becomes evident that it corresponds to the Laplace transform of the second term in (4.15a) and (4.17a). However, attempting to ILT the first term in (4.18) yields non-physical results due to the poles at $s_{j}=\pm \sinh {\chi _{js}}$. Thus, to obtain a regular solution, this term must be nullified at those points. This requirement for regularity provides us with a set of equations in the form
where ${\bar {u}}_{ti}^{\infty }$, ${\bar {u}}_{ri}^{\infty }$, ${\bar {u}}_{ti}^{\infty }$ and ${\bar {v}}_{ri}^{\infty }$ are unknown quantities referring to the asymptotic normal and lateral velocities at both sides of the contact interface. The system of equations completes with use made of the continuity of normal velocity, ${\bar {u}}_{ti}^{\infty }=\beta \,{\bar {u}}_{ri}^{\infty }$, along with the conservation of tangential momentum at large times, $\beta ^2( \bar {v}_{ti}^{\infty }-\bar {v}_{ti}^{0})=\vartheta (\bar {v}_{ri}^{\infty }-\bar {v}_{ri}^{0})$. Simple manipulation allows us to express the asymptotic normal velocity at the interface with
that renders the same value as (4.16), the latter being similar to (40) of Wouchuk & Nishihara (Reference Wouchuk and Nishihara1997), which was further analysed in Wouchuk (Reference Wouchuk2001).
If we examine the two terms in (4.20), the first term represents the initial tangential velocities at the interface, which induce torque and subsequently enhance the interfacial normal velocity, thereby promoting the growth of ripple amplitude. The second term refers to the pressure wave induced by the oscillations of the shock fronts. This term can either augment ${\bar {u}}_{ti}^{\infty }$ or subtract from it, sometimes resulting in a freeze-out at the interface and halting the growth of interfacial instability (Fraley Reference Fraley1986; Mikaelian Reference Mikaelian1994; Wouchuk & Nishihara Reference Wouchuk and Nishihara2004; Wouchuk & Sano Reference Wouchuk and Sano2015). One might be tempted to assert that the calculation of the growth rate is only valid in the absence of SAE at the shock fronts, suggesting that the presence of acoustic disturbances in the long-term regime, as seen in cases (ii) and (iii), can affect the asymptotic growth rate. For case (ii), since non-decaying acoustic contributions do not reach the interface, they have no impact on the contact interface dynamics. In case (iii), where constant-amplitude acoustic disturbances do reach and cross the contact interface, the result remains the same. The mean value of this contribution is strictly zero, meaning there is no imbalance in the history of the acoustic influence for this non-decaying mode. Therefore, the asymptotic growth rate is also $\bar {u}_{ti}^{\infty }$, though the asymptotic dynamic response is $\bar {u}_{ti}^{\infty }+Q_{t}\sin (\omega _{t}\tau _{t})$. Thus, even when the interface oscillates without damping, oscillations have a mean value, which is $\bar {u}_{ti}^{\infty }$.
Figure 6 shows the asymptotic velocity fields, with figure 6(a) displaying the longitudinal components and figure 6(b) showing the transverse components. These are further divided into potential (dashed orange line) and rotational (dashed green line) contributions. Figure 6(c) presents a vector plot of the asymptotic velocity field overlaid on the two-dimensional vorticity field. The parameters are chosen to meet Yang et al. (Reference Yang, Zhang and Sharp1994) results. Dimensionless velocity disturbances on the reflected side are scaled by $\beta$ to match the normalisation of those on the transmitted shock side. This adjustment ensures that the continuity of the longitudinal contribution across the interface is easily identifiable in figure 6(a). As expected, the velocity field near the interface, which is crucial for the development of the RMI, is governed by the potential contributions, while the far-field domain is dominated by the rotational contributions. This is also noticed in figure 6(c), where the misalignment between the eddies and the vorticity peaks close to the interface appears due to the potential contribution to the velocity field.
In the linear regime, the evolution of the rippled interface is symmetric. However, as the distortion increases, weakly nonlinear effects become significant, leading to the emergence of distinct bubbles and spikes. The former refers to the lighter fluid penetrating into the heavier fluid, corresponding the vertical positions where the velocity peaks and $\bar {u}_{ti}>0$, namely $\bar {y}=0, \bar {y}= \pm 2{\rm \pi},\ldots$, in our reference frame. The term spike refers to the heavier fluid penetrating into the lighter fluid, corresponding the vertical positions where the condition $\bar {u}_{ti}<0$ peaks: $\bar {y}=\pm {\rm \pi}, \bar {y}=\pm 3 {\rm \pi},\ldots$. To evaluate whether the distinguished acoustic fields presented in figure 2 actually affect the evolution of the nonlinear interface, an estimation of the characteristic times must be done. From one side, linear theory can estimate the time at which nonlinear effects enter into play, $\xi _i(\tau _t)\sim \varepsilon ^{-1}$, which is found to occur for $\tau _t\sim 1/ (\varepsilon \bar {u}_{ti}^{\infty })$. On the other hand, the time when the acoustic influence vanishes is of the order of unity, see Velikovich et al. (Reference Velikovich, Herrmann and Abarzhi2014) and Probyn et al. (Reference Probyn, Williams, Thornber, Drikakis and Youngs2021). Therefore, if $\varepsilon$ is sufficiently small, nonlinearity – manifested as the distinction between spikes and bubbles – becomes significant only after the asymptotic RMI growth rate is reached. Regardless of the scenario, the combined rotational and acoustic contributions to the velocity field will affect the evolution of bubbles and spikes differently depending on the side of the interface, with SAE amplifying this difference.
5. The RMI at the contact interface
5.1. Prescriptions and approximated models
For many decades, various models have been developed to estimate the linear asymptotic growth rate. These models are referred to as prescriptions. They provide approximate but explicit formulae for evaluating the constant asymptotic growth rate. Such a formula, also proposed by Richtmyer (Reference Richtmyer1960), is derived by treating the RMI as a limiting case of the Rayleigh–Taylor instability, with instantaneous rather than constant acceleration, assuming that perturbation fields are incompressible. According to his model, Richtmyer's formula should be applied using preshock parameters, namely
where $\xi _i^{0-}=1$ with the normalisation chosen. However, to improve agreement with numerical calculations, Richtmyer himself suggested using postshock values instead of preshock ones, resulting in the widely known formula
where $\mathcal {A}_2$ is the postshock Atwood number and $\xi _i^0$ is the postshock interface ripple amplitude explicitly given in (S.2.11) of the supplementary material. For ideal-gas configurations, the impulsive models work relatively well for weak shocks and/or high adiabatic indexes, i.e. when compressible effects are negligible, see figure 1 in Velikovich et al. (Reference Velikovich, Herrmann and Abarzhi2014). For heavy-to-light accelerations, Meyer & Blewett (Reference Meyer and Blewett1972) empirically determined that $\xi _i^0$ should be replaced by the average of the initial unshocked and shocked amplitudes to align with their numerical results, $\bar {u}_{i}^{MB}= (\xi _i^{0-}+\xi _i^0)\mathcal {M}_{2t}\mathcal {A}_2/2$, where $\xi _i^{0-}=1$ with our variables normalisation. More accurate formulations have been proposed, as the prescription by Vandenboomgaerde et al. (Reference Vandenboomgaerde, Mügler and Gauthier1998) that combines both Richtmyer and Meyer–Blewett expressions,
For most cases tested in the literature, the second term of the right-hand side of (5.3) appears to be very small.
In addition to Richtmyer's impulsive model, another alternative impulsive model considers two different accelerations, as developed by Wouchuk & Nishihara (Reference Wouchuk and Nishihara1996). In this case, the growth rate corresponds to the first term of (4.20), which accounts only for the initial velocity shear deposited at the interface after the incident shock refraction,
where the contribution of transverse velocities at the interface is determined by the initial shock ripple amplitudes. This two-acceleration formula is an approximation that neglects the history of the transient acoustic disturbances, but it offers improved predictions for both reflected shocks and rarefaction cases compared with earlier incompressible models (Cobos-Campos & Wouchuk Reference Cobos-Campos and Wouchuk2016; Probyn et al. Reference Probyn, Williams, Thornber, Drikakis and Youngs2021; Li et al. Reference Li, Chen, Zhai and Luo2024). The irrotational Wouchuk–Nishihara model was later extended by the same authors to include bulk vorticity, as seen in (34) of Wouchuk & Nishihara (Reference Wouchuk and Nishihara1997), which is formally similar to (4.20) and (4.16) presented in this work. However, while these equations are analytical and exact, they cannot be expressed in a fully explicit form due to the impossibility to get an explicit expression for the corresponding Laplace transform: it requires the iterative evaluation of the functions $\mathbb {P}_{ts}$. Previously, Yang et al. (Reference Yang, Zhang and Sharp1994) developed a compressible linear model derived from the Euler equations, which allows for the calculation of the exact linear growth rate. However, unlike the posterior works by Wouchuk & Nishihara (Reference Wouchuk and Nishihara1997) and Wouchuk (Reference Wouchuk2001), the model developed by Yang et al. (Reference Yang, Zhang and Sharp1994) is not analytical and requires a numerical evaluation, except for the weak-shock limit that allows for explicitness. There is still potential for (4.20) to be further developed to achieve more accurate predictions in an explicit form. For example, Cobos-Campos & Wouchuk (Reference Cobos-Campos and Wouchuk2016) proposed an approximation of the second term on the right-hand side of (4.20) based on the first guess of the iterative process. Optimising this choice for the term that accounts for the compressible history of the interface can lead to very accurate predictions for the asymptotic growth rate. The functions involved in this second term can be expressed as in (A21), where $\alpha _{j}$ and $\mathbb {P}^{iso}_{js}$ are explicit, accounting for the pressure field generated by the reflected and transmitted shocks as if they were isolated. Therefore, the only non-explicit expressions are the acoustic-coupling functions $F_{t}^{+}$ and $F_{r}^{-}$, which require an iterative scheme to be fully calculated, see (A19). For weak shocks, only the first order of iteration is needed, as the acoustic coupling is negligible. However, since the effect of pressure disturbances generated by the oscillating shocks is only important very close to the shock fronts, except for conditions where SAE occurs, this approximate solution offers very good results even for finite strength shocks, as shown below.
The explicit approximated solution where the pressure field generated by the shocks is taken as in isolated-shock cases reads as
where $\mathbb {P}_{js}^{iso}$ are given in (A21), this expression gives a quick and more accurate evaluation of the asymptotic growth rate in most of the cases. Another explicit expression can be constructed that renders the most accurate results compared with all other approximations for the asymptotic growth rate. It is readily obtained by considering the first guess $[k=0]$ to the acoustic functions, so that ${\bar {u}}_{ti}^{[0]}= {\bar {u}}_{i}^{ISO} + \delta \bar {u}_i^{[0]}$, where
where $F_{t}^{+,[0]}$ and $F_{r}^{-,[0]}$ are given in (A20). The price to pay for this more complicated expression is compensated by the accuracy gained, as observed in figures 7 and 8, respectively. Nevertheless it is easy to see that the solution converges to the exact value of $\bar {u}_{ti}^\infty$ as the iteration index $k$ increases: ${\bar {u}}_{ti}^{\smash {[k]}}\rightarrow \bar {u}_{ti}^\infty$, for $k\gg 1$, so the expression converges to (4.20), or the equivalent (34) from Wouchuk & Nishihara (Reference Wouchuk and Nishihara1997), as the iteration step advances. The advantage of ${\bar {u}}_{ti}^{\smash {[k]}}$ is that it starts from a physically coherent scenario (isolated shocks), and the effects of acoustic coupling are incorporated iteratively, with the first correction being fully explicit.
Another approximate solution is the one given by the sequential model developed in Appendix A.5, which formally applies for conditions where $\beta ^2/\vartheta \gg 1$. For that case, the value of the asymptotic growth rate is that given by (4.20), where the pressure functions in the second term of the right-hand side are given by (A27) for ${\mathbb {P}}_{rs}$, and by the corresponding combination involving the acoustic function (A28) for ${\mathbb {P}}_{ts}$. This sequential-model solution will be denoted as ${\bar {u}}_{ti}^{SM}$.
Before examining the temporal evolution of the unstable contact interface in the linear regime, we conduct a direct comparative analysis of various prescriptions and approximated models for the long-time asymptotic growth rate, including the exact solution provided by (4.16) (or (4.20), equivalently). This analysis is illustrated in figure 7, where we present a shock-strength dependence analysis for certain preshock thermodynamic conditions of the two materials or materials. Figure 7(a) corresponds to the gas-to-air SF$_6$ case, as shown in figures 2(a) and 2(b). Results are aligned with Li et al. (Reference Li, Chen, Zhai and Luo2024), where they find, for the heavy-to-light case, that ‘the incompressible linear model loses validity when the incident shock is strong’. Figures 7(b) and 7(c) illustrate the shock-strength variation for the cases computed in figures 2(c) and 2(d), which are associated with the vdW EoS. Finally, figure 7(d) presents a shock-strength dependence study for Al–Cu, as depicted in figures 2(e) and 2(f), where the values for $\mathcal {M}_{i}=1.3$ and $\mathcal {M}_{i}=2.5$ correspond to the linear growth rate of the nonlinear ‘light-to-heavy’ configuration addressed numerically by Tahir et al. (Reference Tahir, Shutov, Zharkov, Piriz and Stöhlker2011); see § 5.3 for a more elaborated discussion.
Through direct inspection, we can draw both general and specific conclusions. First, as expected, all solutions converge to the linear shock-strength dependence, $\bar {u}_{ti}^\infty \sim \mathcal {M}_i$, in the weak-shock limit $\mathcal {M}_i-1\ll 1$. Second, the approximated solution ${\bar {u}}_{ti}^{\smash {[0]}}$ although complex in form, is fully explicit and is the most accurate across the entire range of parameters explored with various EoS. In some cases, the discrepancies are nearly imperceptible within the order-of-unity scale range. Third, the simplest impulsive-model prescription, based on preshock properties, consistently renders the least accurate predictions. Further conclusions are more specific. The approximation ${\bar {u}}_{ti}^{WN}$ is generally the second-best, except for ideal gases in relatively strong-shock conditions, where ${\bar {u}}_{ti}^{SM}$ and ${\bar {u}}_{ti}^{R}$ provide slightly better approximations. The prescription ${\bar {u}}_{ti}^{VMG}$, commonly applied to the reflected rarefaction scenario, is demonstrated here only under conditions where $\beta < 1$ (see figure 7b,c). The results indicate that the accuracy of prescriptions varies depending on the specific conditions, with no definitive trend favouring one model over others. This variability suggests that caution should be taken when applying them, particularly in the context of non-ideal EoS.
Upon close inspection of the inset in figure 7(d), which corresponds to the three-term EoS for Al–Cu, it can be observed that the slope of the weak-shock linear dependence varies depending on the model being assessed. This variation arises because the three-term EoS lacks sufficient accuracy in the weak-shock limit, where the state functions fail to converge to the imposed cold conditions, such as the reference sound speeds, $c_{0r}$ and $c_{0t}$. The sensitivity of the cold-condition approach depends on specific parameters, such as $\mathcal {R}$ or $\mathcal {M}_2$, resulting in each model being affected differently in this weak-shock regime. However, this is not a major concern, as the weak-shock limit holds minimal relevance in solid RMI.
5.2. The exact transient model
We have observed that prescriptions and approximated models perform relatively well under weak-shock conditions. However, significant deviations from the exact asymptotic solution occur at finite shock strengths. Since the RMI pertains to a transient growth where nonlinear effects can become significant as the contact interface amplitude increases, it is crucial to include the exact temporal evolution in the study.
The exact evolution of the ripple amplitude, (2.6), can be derived from the integration of the velocity field evaluated at the interface, since we are interested in the point where the amplitude is greater the growth rate of the interface ripple is given by
taking the Laplace transform in time of the above equation the solution take the form of
where $\xi _{i}^{0}$ is the dimensionless postshock amplitude at $\tau _t=0^+$, provided in (S.2.11) of the supplementary material. With use made of the ILT, along with the convolution theorem for the second term, we obtain
allowing us to write the temporal evolution of the ripple amplitude with an integral expression, since $\mathbb {L}_{ti}={\mathbb {L}}_{t}(s_{t},0)$ is already determined. However, for sufficiently large times $\tau _t\gg 1$, the asymptotic linear growth rate can be used to write
where the offset
is given through the Taylor series expansion of (5.8) in the variable $s_{t}$. The value of the second term is obtained by taking the derivative of the functional expression $\mathbb {L}_{ti}$ evaluated at $s_{t}=0$, as done in Cobos-Campos & Wouchuk (Reference Cobos-Campos and Wouchuk2016).
The above expressions are valid for cases (i) and (ii). In case (iii), where SAE disturbances reach and cross the contact interface, the result for the ripple evolution is similar but includes a new oscillatory term. Specifically, it can be readily found through straightforward construction that
i.e. an oscillatory term must be included. Likewise, we obtain
for the asymptotic long-time regime.
To analyse the transient evolution of the contact interface, we perform calculations of the interface dynamics under the same conditions as in figure 2. The results, presented in figure 8, indicate that in most cases, the known expressions are not suitable for predicting the RMI, in agreement with the discussion following figure 7. Additionally, it is worth noting the multifrequency character of the solution associated with the acoustic coupling with the oscillating reflected and transmitted shocks. For example, the solid line in figure 8(a) shows the interface evolution in the linear regime, along with the circle symbols corresponding to the case computed in Yang et al. (Reference Yang, Zhang and Sharp1994). Given that the application of the ILT in this specific case is novel, we believe the Yang results, based on the same physical assumptions as our model but integrated with a different numerical method, are well-suited for benchmarking our model. Note that in our figure, the vertical axis represents the dimensionless interface ripple growth rate and incorporates a factor of $\mathcal {M}_i c_{0r}/c_{2t}=3.075$ compared with that in figure 7 of Yang et al. (Reference Yang, Zhang and Sharp1994). Meanwhile, the horizontal axis, associated with the temporal scale, is the inverse. Figure 8(b), depicting a similar air-SF$_6$ configuration but with an stronger incident shock, demonstrates similar qualitative dynamics. In contrast, figure 8(c,d), which correspond to cases using the vdW EoS, exhibit noticeable differences in the frequency and amplitude of oscillations. These differences arise from the distinct characteristics of the two configurations: figure 8(c) depicts a neutrally stable shock with an oscillating regime that does not decay over time, while figure 8(d) illustrates a postshock situation where $\beta <1$. For this case, is the reflected shock the one producing stronger oscillations, as seen in figure 2(d). Please note that we have not included the line corresponding to Richtmyer's postshock formula, as it yields a negative value ($\bar {u}_{i}^{R_2} = -0.0297$). The curve $\bar {u}_{i}^{SM}$ has also been omitted for the similar reasons. Finally, figure 8(e,f) pertain to Al–Cu configurations under different shock strengths. The stronger-shock case results in a much smoother curve with significantly lower oscillations, consistent with the pressure field results shown in figures 2(e) and 2(f). A more detailed explanation of the effect of the EoS modelling for Al–Cu configurations is provided below.
5.3. The impact of the EoS choice
In the previous subsection we have analysed the accuracy of the approximated models or prescriptions when compared with the exact solution, regardless of the EoS employed in the analysis. The present work, which is accommodated to employ any EoS, can be used to analysed whether $\gamma$-like assumptions are sufficiently accurate, see Holmes et al. (Reference Holmes, Dimonte, Fryxell, Gittings, Grove, Schneider, Sharp, Velikovich, Weaver and Zhang1999), Ishizaki & Nishihara (Reference Ishizaki and Nishihara1997) and Endo et al. (Reference Endo1995), where a fitted or effective value of $\gamma$ and initial pressure $p_0$ in an ideal gas EoS are selected to match the post shock values of pressure $p_2$ and velocity $U_2$ observed in real experiments, or Ward & Pullin (Reference Ward and Pullin2011), where the authors carry out a dimensional analysis in order to find complete similarity between Mie–Grüneisen and perfect gas EoS modelled flows.
To conduct this analysis, we examine the case of RMI in an Al–Cu system, similar to the numerical results presented by Tahir et al. (Reference Tahir, Shutov, Zharkov, Piriz and Stöhlker2011) (case II), which also correspond to the cases shown in figure 2(e,f). The incident shock Mach numbers associated with incident shock pressures of approximately $p_1= 0.28$ Mbar and $2.9$ Mbar are $\mathcal {M}_{i}\sim 1.3$ and $\mathcal {M}_{i}\sim 2.5$, respectively, according to the three-term EoS model. Although the base-flow conditions are detailed in figure S3 of the supplementary material, here we analyse the ‘equivalent’ ideal gas problem. Due to the impossibility of finding an initial preshock combination of parameters for $\gamma$-like EoS that exactly matches the base-flow results with more complex EoS (as the three-term description used in this work), i.e. due to the lack of Riemann problem similitude for different EoS (Quartapelle et al. Reference Quartapelle, Castelletti, Guardone and Quaranta2003), we adopt an alternative approach. In particular, the procedure for setting up effective initial parameters for perfect-gas-like problems involves finding the best fit for both EoS zero-order profiles. The method is as follows: first, the preshock pressure $p_0$ and the adiabatic index $\gamma _r$ of the fluid through which the incident shock travels are determined by the density and pressure jumps across the shock. Then, the effective adiabatic indices $\gamma _t$, along with $\mathcal {M}_r$ and $\mathcal {M}_t$, are calculated to match the conservation of material velocity at the interface, given a fixed preshock Atwood number and the density jumps across the reflected and transmitted shocks. This procedure has been found to be very effective, at least for this Al–Cu combination considered in this work, as demonstrated in table 1, where all the base-flow density jumps across the shocks, and the postshock Atwood number and dimensionless interface corrugation are fairly similar between three-terms EoS and $\gamma$-like models. Specifically, for the Al–Cu case with $\mathcal {M}_i=1.3$, the equivalent postshock pressure conditions approximate those of an ideal gas with effective values of $\gamma _r=\gamma _{\rm Al}=3.64$, $\gamma _t=\gamma _{\rm Cu}=5.44$, with a preshock pressure $p_0=0.1366$ Mbar and the same initial density jump at the interface, characterised by $\mathcal {A}_0=0.53436$. Similarly, for the Al–Cu case with $\mathcal {M}_i=2.5$, the effective values are $\gamma _{\rm Al}=2.52$, $\gamma _{\rm Cu}=3.14$ and $p_0=0.3382$ Mbar.
The exact temporal evolution of the interface ripple velocity and its corresponding asymptotic value are shown in figure 9 using the three-term EoS (black) and the ideal-gas simplification (red) for the cases represented in table 1. Precisely, in the cases presented in figure 9(a,b), the ideal gas model underestimates the growth rate by 8 % and overestimates it by 28 %, respectively. Additionally, for figure 9(a), the ideal gas model predicts a slightly higher frequency and amplitude in the ripple velocity time evolution compared with the three-terms model. In figure 9(b), the frequency difference is minimal, but the amplitude increases because the temporal evolution of the three-terms model becomes quasiplanar after the velocity reaches its asymptotic value. Moreover, if a prescription for the asymptotic growth rate is used in addition to the simplification of the EoS, the discrepancies are expected to increase further, as demonstrated in figure 8.
Note that, although Tahir et al. (Reference Tahir, Shutov, Zharkov, Piriz and Stöhlker2011) examines a similar scenario of a shock incident on an Al–Cu interface, a direct comparison between the asymptotic growth rate and the numerical results there shown is not feasible due to differences in the regimes analysed. In particular, the numerical set-up in their study corresponds to a fully nonlinear regime where, in addition, elastic–plastic effects are considered. It is there shown that, for weak shocks, the interface grows from $\psi _i^{0^-}\sim 25\ \mathrm {\mu }{\rm m}$ to $\psi _i(t^*)\sim 77\ \mathrm {\mu }{\rm m}$, in a time scale that ranges $t^*\sim 800$ ns. In our dimensionless variables, this corresponds to $\xi (\tau ^*)\sim 3$, where the dimensionless time domain between frames corresponds to $\tau ^* \sim 270$, as $\lambda \sim 100\ \mathrm {\mu }{\rm m}$ and $c_{2t}\sim 5.5\ {\rm km}\ {\rm s}^{-1}$ for the weak shock case ($\mathcal {M}_i\sim 1.3$). Therefore, even if the initial amplitude were significantly smaller – which is not the case, as $\varepsilon \sim 1$ – the long-time behaviour would likely correspond to a nonlinear regime in practical conditions. The strong shock case considered in Tahir et al. (Reference Tahir, Shutov, Zharkov, Piriz and Stöhlker2011) pertains to the nonlinear regime too. Motivated by the lack small-amplitude studies, this work aims to address such shortage of exact verification solutions in RMI in solids. To facilitate this, we provide, not only an analytic description detailed here and in the supplementary material, but also Mathematica notebooks for both base-flow and perturbation problems. Note that, although the linear RMI problem presented here cannot be compared with the results in Tahir et al. (Reference Tahir, Shutov, Zharkov, Piriz and Stöhlker2011), the base-flow solutions do allow for some direct comparisons, as discussed in the supplementary material.
6. Conclusions
This work advances the well-studied problem of the RMI in the linear regime by presenting an analytical, fully compressible theory for RMI involving reflected shocks with arbitrary EoS. The key advancements are that on one hand, the initial value problem perturbation model is solved using the ILT. In contrast to the previously used separation of variables technique, which describes the pressure field in terms of Fourier–Bessel series, the method proposed in this work offers accurate solutions even under varying conditions and over long-time regimes. On the other hand, this study advances previous work by not being restricted to any particular EoS. Specific theoretical computations are conducted for an ideal gas, a vdW gas and three-term constitutive equations for simple metals. The problem formulation also includes the possibility of SAE for any of the reflected or transmitted shocks and the study of its influence on the compressed material and, eventually, on the contact interface. Additionally, an explicit expression for the asymptotic growth rate is presented, demonstrating excellent agreement across a wide variety of conditions, including variations in incident shock intensity and types of EoS. Two major conclusions are drawn, as follows.
(i) When comparing heuristic prescriptions with the exact model for the linear evolution of the interface, it is found that the former have a very narrow limit of validity when the shock is not sufficiently weak, regardless of the choice of EoS.
(ii) When comparing realistic EoS with $\gamma$-like EoS, artificially constructed to mimic zeroth-order conditions, it is found that such a EoS simplification offers low accuracy concerning the asymptotic growth rate and the characteristic dynamical properties of the transient evolution. Differences increase with the shock strength. Such comparison has been done with use made of the exact linear model.
(iii) The combined use of heuristic simplifications and $\gamma$-like EoS models rarely yields accurate results for the linear evolution of the interface. To reliably predict the contact interface dynamics and growth rates for the RMI, it is crucial to develop precise analytical solutions and realistic EoS models. As shown in this work, achieving this requires extensive algebraic calculations.
This study opens several opportunities for extensions and potential improvements that the authors consider highly significant. First, the heavy-to-light configuration should be addressed in similar terms. Second, for both the light-to-heavy and heavy-to-light cases, the effect of a non-ideal EoS in the nonlinear regime can be investigated using the incompressible model (see Velikovich & Dimonte (Reference Velikovich and Dimonte1996), Zhang & Sohn (Reference Zhang and Sohn1999), Matsuoka, Nishihara & Fukuda (Reference Matsuoka, Nishihara and Fukuda2003b), Matsuoka, Nishihara & Fukuda (Reference Matsuoka, Nishihara and Fukuda2003a), Velikovich et al. (Reference Velikovich, Herrmann and Abarzhi2014), Matsuoka & Nishihara (Reference Matsuoka and Nishihara2020) and references therein). Third, the use of more complex materials as those undergoing elastic–plastic regimes (Piriz et al. Reference Piriz, López Cela, Tahir and Hoffmann2008). Finally, multimode and non-periodically symmetric configurations should be explored to better simulate realistic scenarios (see Groom & Thornber (Reference Groom and Thornber2019), Thornber et al. (Reference Thornber, Griffond, Bigdelou, Boureima, Ramaprabhu, Schilling and Williams2019), Probyn et al. (Reference Probyn, Williams, Thornber, Drikakis and Youngs2021) and Dimonte et al. (Reference Dimonte, Nagori, Ramaprabhu and Boureima2024)).
Supplementary materials
Supplementary material, including detailed mathematical derivations of the base-flow solution for the van der Waals gas and simple metals equations of state, as well as Mathematica codes for calculating base-flow parameters for various EoS and a Mathematica script for the linear evolution of the Richtmyer–Meshkov instability, are available at https://doi.org/10.1017/jfm.2024.1010. Further data can be provided upon request.
Acknowledgements
The authors wish to thank Professor J.G. Wouchuk for carefully reading the manuscript and for his valuable comments.
Funding
This work has been also supported by projects PID2022-139082NB-C51 and TED2021-129446B-C41 funded by MCIN/AEI and also by the Madrid Government (Comunidad de Madrid-Spain) under the Multiannual Agreement with UC3M (H2SAFE-CM-UC3M). F.C. has received support from MCIN/AEI under grant no. PID2021-125550OB-I00. A.L.V. was supported by the National Nuclear Security Administration of the U.S. Department of Energy.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Resolution via the Laplace transform
This appendix presents the resolution of the initial value problem enunciated before through the Laplace transform method. As in previous works (Zaidel’ Reference Zaidel’1960; Wouchuk & Nishihara Reference Wouchuk and Nishihara1996; Cobos-Campos & Wouchuk Reference Cobos-Campos and Wouchuk2016; Huete & Vera Reference Huete and Vera2019; Calvo-Rivera et al. Reference Calvo-Rivera, Velikovich and Huete2023), it is useful to apply a change of variables that transforms the triangular domain of the initial value problem into a square-shaped domain through a hyperbolic stretching at the origin. Subsequently, the Laplace transform is applied.
A.1. The hyperbolic transformation
The relation between the Cartesian coordinates with the hyperbolic coordinates is given by (3.2a,b) that can be applied to the sound wave equations and the corresponding boundary and initial conditions. As per the former, (2.10) takes the form
where the function
has been conveniently introduced to facilitate the manipulation of the equations, as demonstrated later on. Likewise, the boundary conditions at the transmitted shock, see (2.13a) and (2.13b) for the pressure, and (2.12a) for the ripples of the transmitted and reflected side, respectively, and
respectively, where the coefficients $\sigma _{1j}$, $\sigma _{2j}$ and $\sigma _{3j}$ have similar form as those in Huete & Vera (Reference Huete and Vera2019) and Calvo-Rivera et al. (Reference Calvo-Rivera, Velikovich and Huete2023), specifically
The mechanical equilibrium at the interface is
Notice that, although the relation between the $\eta _r$ and $\eta _t$ depends on the corresponding trajectories $\chi _r$ and $\chi _t$, they can be uniquely related at the interface, where $\cosh (\chi _{r}=0)=\cosh (\chi _{t}=0)=1$, resulting in the relationship $\eta _{r}=\beta \eta _{t}$.
A.2. The Laplace transform
The analysis proceeds by employing the Laplace transform on the variable $\eta _{j}$, which corresponds to the scaled temporal variable $\tau _j$ along a specific fixed trajectory $\chi _j$. In particular, the Laplace transforms over the functions ${\bar {p}}_{j}$ and ${\bar l}_{j}$ at both sides of the contact interface are
Following the same structure as before, we first present the Laplace transform of the sound wave equation, followed by the shock and interface boundary conditions. Regarding the former, and using the change of variable $\sinh {q_{j}}=s_{j}$, one can write the expressions (A1) and (A2) as a dynamical system with the form (Briscoe & Kovitz Reference Briscoe and Kovitz1968)
which applies at both sides of the interface. The RH conditions for the transmitted and reflected shocks are, respectively,
Notice that the time scale factor $\beta$ is used here to link the frequency variables in the form $\sinh {q_{t}}=\beta \sinh {q_{r}}$. The same transformation can be applied for the continuity of pressure and streamwise velocity at the interface (A4a–c), which is written in terms of the variable $q_j$ reads as
Equations (A7)–(A9a,b) form a self-contained system that can be numerically integrated to yield the Laplace transform of the pressure field. Nevertheless, there is room for additional analytical derivation, as the solution to the system described in (A7) can be expressed as a combination of acoustic functions, much akin to the d'Alembert solution (see Briscoe & Kovitz (Reference Briscoe and Kovitz1968)) for the 1-D wave equation, namely
where $F_{j}^{-}$ and $F_{j}^{+}$ denote unknown functions that characterise sound wave propagation in the negative and positive directions of $\chi _j$, respectively. To clarify, $F_{j}^{-}$ represents left-travelling sound waves originating from the transmitted shock and heading towards the reflected shock, while $F_{j}^{+}$ corresponds to right-travelling sonic waves. The main advantage of considering (A10) for $F_{j}^{-}$ and $F_{j}^{+}$, over (A7) for $\mathbb {P}_{j}$ and $\mathbb {L}_{j}$, is that the system of differential equations is converted into a system of functional equations (Wouchuk Reference Wouchuk2001).
Referring to (A8), we can express the equations in relation to the acoustic function through (A21), with use made of (A10) to write $\mathbb {P}_j$ as function of the two pair of acoustic functions $F_{j}^{\pm }$,
where $\alpha _{j}(q_{j})$ are auxiliary functions that represent the downstream acoustic influence, incorporating the corresponding Doppler-shift effects between the two shocks and the interface,
and where $\mathbb {P}_{js}^{iso}$ are the isolated-shock pressure contributions,
with the minus sign in (A21b) being given by the change in the system of reference with respect to isolated shock configurations. In addition, the following definitions have been included: $\sigma _{bj}=\sigma _{2j}$, $\sigma _{cj}=\sigma _{1j}\sigma _{3j}$ and $\sigma _{dj}=\sigma _{3j}\xi _{j0}$, in accordance with the notation employed in Calvo-Rivera et al. (Reference Calvo-Rivera, Velikovich and Huete2023).
On the other hand, the equilibrium conditions at the contact interface, (A9a,b), can be also rearranged as a function of the pair of acoustic functions in the following form:
with use made of the auxiliary functions
Therefore, the expressions (A11) and (A14) form a closed functional system of four equations that can be easily reduced to a system of two equations,
with the auxiliary functions $\varTheta _{j1}$, $\varTheta _{j2}$ and $\varTheta _{j3}$ defined as follows:
and
for the reflected and transmitted shocks, respectively. In writing the functions $\varTheta _{j1}$, $\varTheta _{j2}$, we have made use of the definitions $\mathbb {P}_{js}^{iso}$ and $\alpha _{j}(q_{j})$, defined in (A12) and (A13). The resolution of the system of equations is far from trivial, as functions depend on themselves shifted with respect to the variable $\chi _j$, which is the mathematical manifestation of the Doppler shift effect between the shocks and the contact interface. It is worth emphasising that functional equations, including linear ones, often defy analytical solutions. The details of the resolution of this functional system are explained below.
A.3. The iteration scheme
The iteration scheme, based on fixed point method from Wouchuk (Reference Wouchuk2001), renders
where the superscript $[k]$ indicates the iteration steep and where the dependence of $\varTheta$ functions on $q$ is deliberately omitted for clarity purposes. Highlighting the non-trivial nature of selecting an iteration scheme, it is important to emphasise that choosing the correct functional expression is pivotal for guaranteeing solution convergence. By referencing the recurrence functions, it becomes evident that an effective iteration path should consider both its own terms and the contribution from its counterpart shock front. Keeping this concept in mind, the iteration scheme (A19) was proposed, as demonstrated in Wouchuk (Reference Wouchuk2001). Additionally, to effectively resolve the problem, an initial estimate for (A19) needs to be supplied. One possibility is the high-frequency limit, where the initial guess is directly given by the limit $q_{j}\gg |\chi _{js}|$ that neglects Doppler shift contributions, obtaining
where the superscript $[0]$ refers to initial guess. With use made of (A20) in (A19), a good solution convergence accuracy is reached within a few iterations, only demanding around five iterations for the most exigent cases associated with high compressible situations (Cobos-Campos & Wouchuk Reference Cobos-Campos and Wouchuk2016).
A.4. Singularities in the complex plane and the possibility of neutral stability
Once the functions $F_{t}^{-}$ and $F_{r}^{+}$ are determined from (A16a) and (A16b), the other functions $F_{t}^{+}$ and $F_{r}^{-}$ can be easily derived from the equilibrium the RH conditions, (A11), which can be gather together to get $\mathbb {P}_{j}(q_{j},\chi _{j})$, with use made of the wavefunction solution (A10a). Therefore, the pressure field at both sides of the interface is fully determined in the Laplace-frequency variables, either $q_{j}$ or $s_j=\sinh q_{j}$. The pressure field within the domain $-\mathcal {M}_{2r}\tau _r\leqslant \bar {x}\leqslant \mathcal {M}_{2t} \tau _t$ is ultimately determined by the corresponding inverse of the Laplace transform of $\mathbb {P}_{j}(q_{j},\chi _{j})$. Likewise, the pressure gradient function $\mathbb {L}_{j}(q_{j},\chi _{j})$ is given by (A10b). In order to obtain the pressure field description in the spatiotemporal domain an analysis of the singularities arisen in the complex plane is needed. Basically, there are two types of singularities that can appear in the complex plane: isolated singularities and branch points which are related with complex functions that are multivalued. The former ones split into three different types depending on the behaviour of the Laurent expansion for the corresponding complex function (Duffy Reference Duffy2004, Reference Duffy2016): essential singularities which are characterised by the fact of have an infinite number of negative-order coefficients; removable singularities which have zero residue; and poles, whose largest inverse (negative) power of the Laurent expansion is finite. In the present case, the only isolated singularities that appear are poles, which can be easily anticipated by the study of the EoS in the same way as in previous works (D'yakov Reference D'yakov1954; Bates & Montgomery Reference Bates and Montgomery1999; Wouchuk & Cavada Reference Wouchuk and Cavada2004; Calvo-Rivera et al. Reference Calvo-Rivera, Velikovich and Huete2023). For example, simple rearrangement of the above Laplace-variable equations renders
From (A21) it is clear that the roots in the denominator are similar to those occurring in the isolated shock cases, since the roots implicitly associated with the functions $F_t^+$ and $F_r^-$ involve lower-frequency roots. The existence of imaginary poles that lie out of the previously defined branch cut affects the stability of the shock fronts, rendering the regime of SAE, which means that the shock fronts exhibit a permanent oscillating behaviour entering in a neutrally stable regime. It is easy to anticipate the existence of poles inside of the Bromwich contour by analysing the denominator of the expressions (A12) and (A13), which in terms of the variable $s_{j}$ is
where $D(s_{j})=0$ is the dispersion relation that one gets by performing the normal mode analysis as D'yakov (Reference D'yakov1954), obtaining the same as Clavin & Searby (Reference Clavin and Searby2016) show in their § 4.4, with the substitution $s={\rm i}\omega$. Direct inspection in the dispersion relation provides limit below which damped oscillations occurs, namely $s_j=\pm {\rm i}$, since imaginary poles are placed below than the branch point value so that they are omitted from the integration path. For $\sigma _{cj}<\sigma _{bj}<\sigma _{cj}+1/(4\sigma _{cj})$ (a condition that shock waves in an ideal gas meet) the corresponding decay follows the following temporal power law $\smash {t^{-3/2}}$. For $\sigma _{cj}+1/(4\sigma _{cj})<\sigma _{bj}$, the long-time power law is the same but the initial damping is primarily exponential (Bates Reference Bates2004). Under this condition, $\sigma _{bj}>\sigma _{cj}$, the only singularities present in this problem are exclusively branch points, mathematically associated with the multivalued character of the functions involving the term $\cosh q_j$ (or $\sqrt {\smash [b]{s_j^2+1}}$), and physically related to the self-oscillation frequency of each shock. When $\sigma _{bj}=\sigma _{cj}$, decaying oscillation changes it envelope to $t^{-1/2}$, at this point one get the critical value of the DK parameter $h_j=h_{cj}$, where
with $j=t,r$ being transmitted or reflected front density jump and postshock Mach number, respectively. Finally, for $\sigma _{bj}<\sigma _{cj}$ in the shock oscillates with constant amplitude in the long-time regime. If this condition applies, the roots of the dispersion relationship
can be split into the imaginary and real parts, namely $s_{I}= \pm [n(1-m)]^{\smash {1/2}}$ and $s_{R}= \pm [n(1+m)]^{\smash {1/2}}$, respectively. Here, the real part, $s_{R}$, provides spurious roots of the dispersion relationship.
Due to the acoustic coupling between the two fluids, an additional pair of branch points will arise in both sides. For example, in the transmitted pressure field the term $\cosh {q_{r}}$ appears, which in terms of $s_{t}$ variable means $[\smash [b]{(s_t/\beta )^{2}+1}]^{\smash {1/2}}$. Then, if $\beta >1$, the frequency induced by the reflected shock is going to exceed the oscillation frequency of the transmitted one, signifying in a extension of the branch cut, again we refer to scheme given in figure 10.
A.5. The sequential model
Therefore, under the limit $\beta ^2/\vartheta \gg 1$, simple manipulation of the Laplace pressure functions of the reflected shock renders
involving the following auxiliary functions:
which can be solved with a recurrence expression
as shown in Wouchuk & Cavada (Reference Wouchuk and Cavada2004), Calvo-Rivera et al. (Reference Calvo-Rivera, Velikovich and Huete2023) and repeated here for convenience. From the above equation, that describes the reflected-shock pressure field, the function $F_{r}^{-}$ can be obtained by direct substitution into the RH condition (A21). Besides, through the pressure relationships at the interface in the asymptotic limit $\beta ^2/\vartheta \gg 1$, namely $F_{r}^{+}(q_{r})=F_{r}^{-}(q_{r})$, the complete pressure field can be constructed within the reflected-shock side. As per the transmitted-shock side, the knowledge of the pressure field at the contact interface can be used to fully determine the transmitted shock evolution. For example, a single functional expression for $F_{t}^{+}$, i.e.
must be resolved, with the auxiliary functions
being expressed in terms of the already known pressure function $F_{r}^{-}$, and also $\mathbb {P}_{ts}^{iso}$ being the translated function given in (A13), and $\alpha _{t}$ being determined in (A12) after the corresponding variable shift. Once $F_{t}^{+}$ is established, the remaining acoustic functions $F_{t}^{-}$, $\mathbb {P}_{t}$ and $\mathbb {L}_{t}$ can be readily determined through straightforward manipulation.