The reflection of a planar impulsive shock wave at a liquid–gas interface

Abstract

The reflection of a shock pulse at a liquid–gas interface occurs in many applications, from lithotripsy to underwater explosions and additive manufacturing. In linear theory, reflection and transmission at an interface depend only on the impedance difference, but this does not hold for a nonlinear pulse. This work develops an analytical framework for computing the reflection and transmission coefficients for an impulsive shock wave at a liquid–gas interface. The problem is treated analytically by considering idealised pulses and solving a series of consecutive Riemann problems. These correspond to the initial interaction with the interface and important subsequent wave interactions that enable a complete description of the process to be obtained. Comparisons with numerical and existing analytical approaches are made for the case of a water–air interface. In the acoustic limit, the method produces results identical to those of linear acoustic theory. As the pulse strength increases, the proposed method agrees well with numerical simulation results, whereas existing analytical methods that consider only the interface fail. We detail how a reflecting pulse can put water into tension without any incident negative pressure. It is further shown that the magnitude of the reflection coefficient decreases with increasing incident shock pressure, and the reflected pulse widens. Reflections of pulses with positive and negative pressures temporarily create negative pressure regions with greater magnitude than the incident pulse. Finally, we consider non-idealised waves. Comparisons with simulations show that the reflection characteristics can be explained qualitatively using the analytical method, and the reflection coefficients are predicted accurately.

1. Introduction

The reflection of a shock wave at a liquid–gas interface is a phenomenon that plays a critical part in various applications that range from medical treatment procedures such as lithotripsy (Cleveland & McAteer 2012) and histotripsy (Maxwell et al. 2011; Xu et al. 2024), to underwater explosions (Holt 1977; Yu et al. 2024), additive manufacturing (Jalaal et al. 2019) and erosion (Obreschkow et al. 2011; Field et al. 2012), among others. More generally, reflections are ubiquitous in virtually all flows involving shock waves and interfaces; e.g. atomisation (Dekel et al. 1998; Avila & Ohl 2016; Stan et al. 2016; Rosselló et al. 2023) or bubble/cavity collapse (Bourne & Field 1992; Johnsen & Colonius 2009; Hawker & Ventikos 2012; Ohl, Klaseboer & Khoo 2015; Bempedelis & Ventikos 2020a; Bokman et al. 2023). These applications cover a large range of shock pressures and waveforms, including impulsive shock waves that have a short, finite width. In many cases, particularly in medical applications, having a precise understanding of the shock–interface interaction and the resulting waves and their amplitudes is of paramount importance.

In shock wave lithotripsy, the secondary cavitation resulting from shock reflections and interactions is known to promote stone fragmentation, but it is also the dominant cause of surrounding tissue injury (Cleveland & McAteer 2012); it thus needs to be controlled carefully. In cavitational histotripsy, mechanical tissue ablation is achieved through wave-incited cavitation cloud activity (Xu et al. 2024). The formation of the cavitation cloud, upon which the successful outcome of the treatment depends, relies strongly on the tension generated by the reflection of the incident waves (Maxwell et al. 2011). The formation of cavitation clouds upstream of a bubble is also observed in boiling histotripsy. Understanding and controlling these clouds is of paramount importance for this type of treatment as well. Pahk et al. (2021) used numerical modelling of the Westervelt equation to investigate how shock scattering off a bubble leads to a cavitation cloud forming upstream, which grows towards the shock source. The results supported the hypothesis that the formation of the cloud is related to the sum of the strength of the scattered shock wave and the incident pressure field, highlighting the importance of understanding not just the immediate reflection at the interface, but also the subsequent wave interactions.

Underwater explosions in shallow water lead to shock interactions with the ocean surface (Holt 1977). In an attempt to mitigate the impact of blast waves on marine life (but also of high-amplitude noise originating from other processes, such as pile driving for offshore wind farm installation), engineers make use of bubble curtains (Timofeev et al. 1985; Croci et al. 2014). Bubble curtains reduce wave transmission by partially reflecting and absorbing the energy of the incoming waves, with their efficiency depending strongly on the properties and distribution of the bubbles (Smith, Bempedelis & Grech La Rosa 2023).

In laser-induced forward transfer, an additive manufacturing technique that uses laser pulses to deposit material from a donor film to an acceptor substrate (Arnold, Serra & Piqué 2007), the reflection of shock waves at the free surface has been suggested as one of two mechanisms responsible for disrupting printing and resulting in irreproducibility and uncertainty (Jalaal et al. 2019).

Shock waves in liquid volumes reflect off the liquid–gas interface, leading to tension in the liquid. Focusing of the reflected waves by the curvature of the liquid volume increases the tension, resulting in the development of cavitation in the liquid (Obreschkow et al. 2011; Avila & Ohl 2016; Schmidmayer & Biasiori-Poulanges 2023). This can lead to jetting (Rosselló et al. 2023), fragmentation/atomisation (Avila & Ohl 2016), and eventually erosion (Obreschkow et al. 2011; Field et al. 2012).

The reflection/transmission of a shock wave at a fluid interface is often approximated by linear acoustic theory (Ohl 2002; Avila & Ohl 2016; Biasiori-Poulanges & El-Rabii 2021). In linear theory, the reflection and transmission of a wave at an interface can be described in terms of the impedance mismatch between the media at either side of the interface. Partial reflection occurs in the case of media with different impedance values. The problem can be formulated as a Riemann problem for the linear acoustic equations, and an exact solution can readily be found (LeVeque 2002). The solution to the Riemann problem is a left-going wave and a right-going wave, which allows for the following reflection and transmission coefficients to be derived:

(1.1)$$\begin{gather} C_\mathcal{R} = \frac{p_\mathcal{R}}{p_\mathcal{S}} = \frac{Z_R - Z_L}{Z_R + Z_L}, \end{gather}$$

(1.2)$$\begin{gather}C_\mathcal{T} = \frac{p_\mathcal{T}}{p_\mathcal{S}} = \frac{2Z_R}{Z_L + Z_R}, \end{gather}$$

where $p$ denotes the pressure, and $Z= \rho c$ is the impedance of each medium. In the above, $\rho$ is the density, and $c$ is the sound speed of the media on the left and right sides of the interface, denoted by the subscripts $L$ and $R$, respectively. The subscripts $\mathcal {R}$ and $\mathcal {T}$ denote reflection and transmission, respectively, and $\mathcal {S}$ denotes the incident shock. The coefficients derived in (1.1)–(1.2) depend only on the impedance mismatch, and are independent of the amplitude of the incident wave. Considering a one-dimensional problem involving water and air at atmospheric pressure, (1.1) yields $C_\mathcal {R} = -0.999$, indicating an almost perfect reflection of the incident wave. However, using a linear approach for a fundamentally nonlinear problem may neglect important elements of the underlying physics.

In the case of a shock wave with constant post-shock conditions behind its front, the solution can be found by solving a Riemann problem at the interface for a nonlinear conservation equation: usually either Burgers’ equation or the Euler equations. This problem was considered by Henderson (1989), who derived reflection and transmission coefficients based on the pressure resulting from the interaction of a shock front with a fluid interface. For an arbitrary equation of state, Henderson (1989) defined

(1.3)$$\begin{gather} C_\mathcal{R} = \frac{p_\mathcal{R} - p_\mathcal{S}}{p_\mathcal{S} - p_0}, \end{gather}$$

(1.4)$$\begin{gather}C_\mathcal{T} = \frac{p_\mathcal{T} - p_0}{p_\mathcal{S} - p_0}, \end{gather}$$

where $p_0$ is the ambient pressure. For an air–water interface, this again yields a coefficient close to $-$1, giving rise to the idea that a shock wave ‘inverts’ when it reflects off such an interface. However, by considering only what happens at the interface, it is implicitly assumed that the incident pulse and the reflected wave do not change as they interact, and that the fluid pressure can be defined by the ambient pressure plus the contributions of the relevant waves. This is true for a linear wave and for constant post-shock conditions, but – as will be shown in this work – a single Riemann problem at the interface is not sufficient to describe the shock–interface interaction for a finite-duration pulse unless the shock is very weak. The problem of a finite-duration shock pulse has received far less attention in the literature from an analytical perspective, but occurs in many practical applications, as described previously. This problem is also challenging for numerical methods because it involves multiple interactions with an interface with a high impedance difference. In many applications, for example those involving shock–bubble interactions, this may involve large-scale separations, with shock waves being significantly narrower than the bubble radius (see Bokman et al. (2023), for example). Analytical methods can therefore be used to help to validate numerical approaches by providing benchmark tests for such problems.

In this work, we develop an analytical method for computing the reflection and transmission of a shock pulse at a liquid–gas interface. The problem is treated by considering the one-dimensional Euler equations and two idealised waveforms: a positive square pulse, and a pulse with positive and negative parts. The solution is obtained by solving a series of Riemann problems corresponding to different wave interactions that take place during the reflection/transmission process. The solution is derived for a general fast–slow problem, and then examined in more detail for a water–air interface at atmospheric pressure. In the acoustic limit, $M\rightarrow 1^+$ (M denoting the Mach number), the method is shown to be exact. The validity of the method for a wide range of shock strengths as well as non-idealised waveforms is examined via comparisons with numerical simulations.

The two idealised waveforms are chosen to represent two different classes of shock pulses: purely compressive pulses and pulses with both compressive and tensile parts. These pulses are prototypical of the types of waves found in the applications mentioned above, but also in many others (see e.g. Bokman et al. 2023). Of practical interest are the maximum and minimum pressures that occur during and after the reflection/transmission process. Together with the complete solution structure, these are described in detail for the idealised waves and then considered both analytically and numerically for two non-square pulses: a modified Friedlander wave and a lithotripter pulse.

The purpose of this work is threefold. First, it presents a complete analytical framework for computing the reflection and transmission of a finite-duration shock pulse at an arbitrary fast–slow interface. Second, the approach is applied to the frequently encountered problem of a finite-duration shock interacting with a water–air interface. Expressions for the reflection and transmission coefficients across a wide range of shock pressures are derived and compared to alternative approaches used in the literature. The discrepancies between simpler but often-used approaches and the one presented here highlight the benefits of this new approach. Third, the approach developed can be used to assess the performance of numerical methods for multiphase flow problems.

The paper is organised as follows. The governing equations and general solution to the multi-fluid Riemann problem are presented in § 2. Section 3 contains the derivation of the complete solution for the reflection and transmission of the idealised shock pulses. The predictions of the developed analytical method are discussed and compared with numerical simulations in § 4. The reflection of non-idealised shock pulses is then considered in § 5, before conclusions are presented in § 6.

2. Physical and mathematical background

2.1. Governing equations and thermodynamics

In this study, we use the compressible Euler equations as a model for the fluid dynamics of both the liquid and gas phases. In one dimension, these are shown in their conservation form as

(2.1) \begin{equation} \frac{\partial \boldsymbol{U}}{\partial t} + \frac{\partial}{\partial x}\boldsymbol{F}(\boldsymbol{U}) = 0, \end{equation}

where $\boldsymbol {U} = (\rho, \rho u, E )^{\rm T}$ denotes the vector of conserved variables, and $\boldsymbol {F}(\boldsymbol {U}) = (\rho u,$ $\rho u^2 + p, u(E + p ) )^{\rm T}$ denotes the flux vector. In the above, $\rho$ denotes the density, $u$ is the velocity, and $p$ is the pressure. The total energy $E$, defined as the sum of the kinetic and internal energies, is

(2.2) \begin{equation} E = \rho(\tfrac{1}{2}u^2 + e), \end{equation}

where $e$ denotes the specific internal energy. To close the equations, we also require an equation of state. In this study, the stiffened gas equation of state (SG EoS) is adopted (Menikoff & Plohr 1989). This is a simplified form of the Grüneisen equation of state, and essentially treats the liquid as an ideal gas that is under high pressure, as defined by a constant $p_\infty$. The specific internal energy is defined as

(2.3)\begin{equation} e = \frac{p+\gamma p_\infty}{\rho(\gamma - 1)}. \end{equation}

For the gas phase, we may still use the SG EoS since it reduces to the ideal gas law when $p_\infty = 0\,\text {Pa}$. For the SG EoS, the speed of sound is defined as

(2.4)\begin{equation} c = \sqrt{\frac{\gamma(p+p_\infty)}{\rho}}. \end{equation}

When considering the propagation of shock waves, the Rankine–Hugoniot relations are most useful for describing the relationship between the states immediately either side of a shock. These are described in a number of texts (e.g. Sochet 2017) and are presented here for a stiffened gas. For the case of a shock wave connecting a ‘shocked’ region $\varOmega _s$ to an ‘unshocked’ region $\varOmega _0$, we have

(2.5)$$\begin{gather} p_s = (p_0+p_\infty)\,\frac{2\gamma M^2 - (\gamma - 1)}{\gamma+1} - p_\infty, \end{gather}$$

(2.6)$$\begin{gather}\rho_s = \rho_0\left(\frac{(\gamma+1)M^2}{(\gamma-1)M^2 + 2}\right), \end{gather}$$

(2.7)$$\begin{gather}u_s = \frac{2c_0}{\gamma+1}\left(\frac{M^2-1}{M}\right), \end{gather}$$

where $M$ denotes the Mach number of the shock wave, and $c_0$ denotes the speed of sound in the unshocked fluid.

2.2. The exact solution for a multi-fluid Riemann problem

The Riemann problem is a well established concept in fluid mechanics that forms the basis of many discretisation schemes, particularly for hyperbolic equations such as the Euler equations (LeVeque 2002; Toro 2013). In this context, we seek an exact solution to the Euler equations with piecewise discontinuous initial data, defined as the left and right initial states. These states are connected by either shock waves or rarefaction waves, and a contact discontinuity depending on the initial conditions. A detailed exposition of the problem can be found in a number of textbooks, e.g. LeVeque (2002) and Toro (2013), and will thus not be given here. A typical solution to the problem in $x\unicode{x2013}t$ space is shown in figure 1. The primitive variables $\boldsymbol {w} = (\rho, u, p )^{\rm T}$ vary across both rarefaction and shock waves. However, it can be shown that pressure and velocity are invariant across the contact discontinuity (Toro 2013), and it is this invariance that enables us to find a solution to the problem.

Figure 1. Typical solution to the Riemann problem showing a right-going shock wave and a left-going rarefaction wave.

While many works consider this problem assuming that the fluid is the same either side of the initial discontinuity, there is nothing in the formulation of the Riemann problem that precludes having different fluids or even different equations of state either side of the interface, provided that the appropriate equation of state is adopted for waves travelling into each fluid (LeVeque 2002). In this study, it is necessary to solve the Riemann problem for both single- and multi-fluid cases, so we seek a general solution that permits the left and right initial states to be different fluids. To this end, the exact solution is derived for the SG EoS, noting that it reduces to the ideal gas law for $p_\infty =0$ Pa. The Riemann problem with the SG EoS was considered by Haller, Ventikos & Poulikakos (2003) to study the wave structure in a droplet following a high-speed impact with a surface. Here, this approach is extended to the more general multi-fluid problem where the fluids either side of the initial interface can have different equations of state.

For a typical Riemann problem, the solution consists of a 1-wave, a contact discontinuity, and a 3-wave. The 1- and 3-waves can be either a rarefaction wave (denoted with the subscript $r$) or a shock wave (denoted with the subscript $s$) depending on the initial data. We seek a solution for the flow variables in the so-called star region that lies between the initial data in the phase plane (see figure 1). To do this, we note that the pressure and velocity are invariant across the contact discontinuity. Therefore, by deriving expressions for $u^*$ in terms of $p^*$, we can obtain an algebraic expression for computing these variables. That is, depending on the wave structure, we seek a solution to

(2.8)\begin{equation} \phi_L(p^*) = u^* = \phi_R(p^*). \end{equation}

For a 1- or 3-rarefaction wave, we can derive functions $\phi (p^*)$ using the Riemann invariants. These are quantities that remain constant across a particular wave. In the case of a rarefaction wave, the following quantities remain constant (Haller et al. 2003):

(2.9a,b)\begin{equation} \frac{p+p_\infty}{\rho^\gamma},\quad u\pm\frac{2c}{\gamma-1}, \end{equation}

where $+$ is for a 1-rarefaction, and $-$ is for a 3-rarefaction. Therefore, for a 1-rarefaction,

(2.10)\begin{equation} \phi_L(p^*) = u^* = u_L + \frac{2(c_L-c_L^*)}{\gamma_L - 1}, \end{equation}

where $c_L$ is known from the left state data, but $c_L^*$ is not. However, we can use the isentropic assumption to express $c_L^*$ in terms of a single unknown $p^*$:

(2.11)\begin{equation} \frac{p_L+p_{\infty,L}}{\rho_L^{\gamma_L}} = \frac{p^* + p_\infty^*}{\rho_L^{*\gamma_L}}. \end{equation}

Rearranging this expression yields

(2.12)\begin{equation} \rho_L^* = \rho_L\left(\frac{p_L+p_{\infty,L}}{p^*+p_{\infty,L}}\right)^{1/\gamma_L}. \end{equation}

Hence

(2.13)\begin{equation} c_L^* = \sqrt{\frac{\gamma_L (p^*+p_{\infty,L})}{\rho_L \left(\dfrac{p_L+p_{\infty,L}}{p^*+p_{\infty,L}}\right)^{1/\gamma_L}}}. \end{equation}

The same approach can be used to derive an expression for a 3-rarefaction wave, yielding

(2.14)\begin{equation} \phi_R(p^*) = u^* = u_R - \frac{2(c_R-c_R^*)}{\gamma_R - 1}, \end{equation}

where

(2.15)\begin{equation} c_R^* = \sqrt{\frac{\gamma_R (p^*+p_{\infty,R})}{\rho_R \left(\dfrac{p_R+p_{\infty,R}}{p^*+p_{\infty,R}}\right)^{1/\gamma_R}}}. \end{equation}

For a 1- or 3-shock, the entropy does change across the wave, so we cannot use the same approach for computing $\phi$. Instead, following the approaches of Haller et al. (2003) and Toro (2013), we consider the conservation laws across the shock wave in the frame of reference moving with the shock:

(2.16a,b)\begin{equation} \hat{u}_L = u_L - c_s,\quad \hat{u}^* = u^* - c_s. \end{equation}

Then the Rankine–Hugoniot relations across the shock can be written as

(2.17)$$\begin{gather} \rho_L\hat{u}_L = \rho_L^*\hat{u}, \end{gather}$$

(2.18)$$\begin{gather}\rho_L\hat{u}_L^2 + p_L = \rho_L^*\hat{u}^{*^2} + p^*, \end{gather}$$

(2.19)$$\begin{gather}\hat{u}_L(\hat{E}+p_L) = \hat{u}^*(\hat{E}_L^* + p^*). \end{gather}$$

Defining the mass flux as

(2.20)\begin{equation} Q_L=\rho_L\hat{u}_L = \rho_L^*\hat{u}^*, \end{equation}

and substituting into (2.18), we obtain

(2.21)\begin{equation} Q_L\hat{u}_L + p_L = Q_L\hat{u}^* + p^* \Rightarrow Q_L ={-}\frac{p^* - p_L}{\hat{u}^* - \hat{u}_L}. \end{equation}

From (2.16a,b), we have that $\hat {u}_L - \hat {u}^* = u_L - u^*$, so

(2.22)\begin{equation} Q_L ={-}\frac{p^* - p_L}{u^* - u_L} \Rightarrow u^* = u_L - \frac{p^*-p_L}{Q_L}. \end{equation}

We now seek an expression for $Q_L$ in terms of only the known initial data and $p^*$. To do this, we must rewrite (2.21) in terms of pressure and density, and then use the equation of state. By noting $\hat {u}_L = Q_L/\rho _L$ and $\hat {u}^*=Q_L/\rho _L^*$ from (2.20), we have

(2.23)\begin{equation} Q_L^2 ={-}\frac{p^* - p_L}{\dfrac{1}{\rho_L^*} - \dfrac{1}{\rho_L}}. \end{equation}

To find an expression for the unknown $\rho _L^*$, we expand (2.19) as

(2.24) \begin{equation} \hat{u}_L\biggl(\frac{\rho_L\hat{u}_L^2}{2} + \rho_Le_L + p_L\biggr) = \hat{u}^*\biggl(\frac{\rho_L^*\hat{u}^{*^2}}{2} + \rho_L^*e_L^* + p^*\biggr), \end{equation}

and by rearranging and using (2.20), we obtain

(2.25)\begin{equation} \frac{\hat{u}_L^2}{2} + e_L + \frac{p_L}{\rho_L} = \frac{\hat{u}^{*^2}}{2} + e_L^*+\frac{p^*}{\rho_L^*}. \end{equation}

It is now convenient to introduce the free enthalpy, defined as $h = e + {p}/{\rho }$. Hence

(2.26)\begin{equation} \frac{\hat{u}_L^2}{2} + h_L = \frac{\hat{u}^{*^2}}{2} + h_L^*. \end{equation}

To obtain expressions for $\hat {u}_L^2$ and $\hat {u}^{*^2}$ in terms of pressure and density, we combine (2.17) and (2.18) and rearrange to obtain

(2.27)$$\begin{gather} \hat{u}^{*^2} = \left(\frac{\rho_L}{\rho_L^*}\right)\left(\frac{p_L - p^*}{\rho_L - \rho_L^*}\right), \end{gather}$$

(2.28)$$\begin{gather}\hat{u}_L^2 = \left(\frac{\rho^*}{\rho_L}\right)\left(\frac{p_L - p^*}{\rho_L - \rho_L^*}\right). \end{gather}$$

Rewriting (2.26) with these expressions yields

(2.29)\begin{equation} h_L - h^* = \frac{1}{2}(p_L - p^*)\left(\frac{\rho_L + \rho_L^*}{\rho_L\rho^*}\right). \end{equation}

Rewriting in terms of specific internal energy yields

(2.30)\begin{equation} e_L - e^* = h_L - h^* + \frac{p^*}{\rho_L^*} - \frac{p_L}{\rho_L} = \frac{1}{2}(p_L + p^*)\left(\frac{\rho_L-\rho_L^*}{\rho_L^*\rho_L}\right). \end{equation}

Substituting in the expression for $e$ from the stiffened gas equation (2.3) and manipulating gives

(2.31)\begin{equation} \rho_L^* = \rho_L\left(\frac{p_L(\gamma_L-1)+p^*(\gamma_L+1)+ 2\gamma_Lp_{\infty,L}}{p_L(\gamma_L+1)+p^*(\gamma_L-1)+2\gamma_Lp_{\infty,L}}\right). \end{equation}

Substituting back into (2.23) and tidying up yields

(2.32)\begin{equation} Q_L^2 = \tfrac{1}{2}\rho_L(p_L(\gamma_L - 1)+p^*(\gamma_L+1)+2\gamma_Lp_{\infty,L}). \end{equation}

Finally, we can write an expression for $\phi _L(p^*)$ by substituting this into (2.22) to give

(2.33)\begin{equation} \phi_L(p^*) = u^* = u_L - \frac{p^* - p_L}{\sqrt{\tfrac{1}{2}\rho_L \left(p_L(\gamma_L - 1)+p^*(\gamma_L+1)+2\gamma_Lp_{\infty,L}\right)}}. \end{equation}

The same process can be followed to obtain the expression for a 3-shock:

(2.34)\begin{equation} \phi_R(p^*) = u^* = u_R + \frac{p^* - p_R}{{\sqrt{\tfrac{1}{2}\rho_R \left(p_R(\gamma_R - 1)+p^*(\gamma_R+1)+2\gamma_Rp_{\infty,R}\right)}}}. \end{equation}

Therefore, using (2.10), (2.14), (2.33) and (2.34), we can solve the algebraic expression

(2.35)\begin{equation} \phi_L(p^*)-\phi_R(p^*) = 0 \end{equation}

to obtain $p^*$. It should be noted that at this stage, we do not know the structure of the solution, so the choice of which functions to use may at first appear ambiguous. Fortunately, it can be shown that any of the expressions for the left and right functions can be used provided that those for the correct phase are applied. This is because there is only one root of (2.35) (Menikoff & Plohr 1989). The resulting expression is a nonlinear algebraic equation that must be solved iteratively. The solution is therefore not, strictly speaking, exact, but given that we can find the solution to an arbitrary degree of accuracy, it is often considered as such.

Once $p^*$ has been found, the wave structure becomes obvious. If $p_L < p^*$, then the left-going wave is a shock. If $p_L > p^*$, then the left-going wave is a rarefaction. The same holds for the right-going waves. The velocity in the star (*) region can be obtained by evaluating any of the individual expressions used to obtain $p^*$, e.g. (2.10). The density in the left and right star regions can then be found, and these depend on the wave structure for the specific problem. If the left- or right-going wave is a rarefaction, then we can compute the density by rearranging the equation for the speed of sound in the left or right star region, which is given by (2.13) and (2.15). Hence for a left-going rarefaction we have

(2.36)\begin{equation} \rho_L^* = \frac{\gamma_L (p^* + p_{\infty,L})}{c_L^{*^2}}. \end{equation}

In the case of a shock wave, the density in the star region can be computed using (2.31).

Now that all of the variables in the star region have been calculated, we need to obtain expressions for the evolution of the waves. For a shock wave, we can use the Rankine–Hugoniot relations. Since we know the pressure either side of the shock, we can rearrange (2.5) to obtain the Mach number and hence the shock speed. As $u^*$ is constant across the contact discontinuity, this must propagate at $u=u^*$. Finally, we must specify the structure of any rarefaction waves. Using the approach set out by LeVeque (2002) but with the SG EoS, we define

(2.37)\begin{equation} \xi = x/t. \end{equation}

For a left-going rarefaction, we know that at each point in the rarefaction $\xi = u - c$, so

(2.38)\begin{equation} c=u-\xi. \end{equation}

For a left-going rarefaction, substituting this into the equation for the Riemann invariant yields

(2.39)\begin{equation} u(\xi) = \frac{u_L(\gamma_i - 1) + 2(c_L + \xi)}{\gamma_i + 1}, \end{equation}

where $i=l,g$ denotes the fluid in which the rarefaction wave is propagating. We can use this expression to obtain an expression for density:

(2.40)\begin{equation} \rho(\xi) = \left(\frac{\rho_L(u(\xi)-\xi)^2}{\gamma_i(p_L + p_{\infty,i})}\right)^{{1}/{(\gamma_i -1)}}. \end{equation}

Finally, considering that entropy is constant across a rarefaction, we can obtain the following expression for the pressure:

(2.41)\begin{equation} p(\xi) = \left(\frac{\rho(\xi)}{\rho_L}\right)^{\gamma_i}(p_L + p_{\infty,i}) - p_{\infty,i}. \end{equation}

The above expressions define the velocity, density and pressure in the rarefaction wave. The domain of this wave is bounded by the speeds of the rarefaction front and tail, which are denoted by the subscripts $F$ and $T$, respectively. The front travels at

(2.42)\begin{equation} \zeta_{r,F} = u_L - c_L. \end{equation}

The speed of the tail is found by solving the velocity equation for $u(\xi )=u^*$, which yields

(2.43)\begin{equation} \zeta_{r,T} = \tfrac{1}{2}(u^*(\gamma_i +1)-u_L(\gamma_i -1)-2c_L). \end{equation}

Because the tail will always move more slowly than the head, the rarefaction wave spreads out over time. A similar process can be followed to obtain the structure of a right-going rarefaction, to obtain

(2.44)$$\begin{gather} u(\xi) = \frac{u_R(\gamma_i - 1) - 2(c_R - \xi)}{\gamma_i + 1}, \end{gather}$$

(2.45)$$\begin{gather}\rho(\xi) = \left(\frac{\rho_R(u(\xi)-\xi)^2}{\gamma_i(p_R + p_{\infty,i})}\right)^{{1}/{(\gamma_i -1)}}, \end{gather}$$

(2.46)$$\begin{gather}p(\xi) = \left(\frac{\rho(\xi)}{\rho_R}\right)^{\gamma_i}(p_R + p_{\infty,i}) - p_{\infty,i}. \end{gather}$$

The front of the right-going rarefaction wave travels at a speed

(2.47)\begin{equation} \zeta_{r,F} = u_R + c_R, \end{equation}

and the tail speed is

(2.48)\begin{equation} \zeta_{r,T} = \tfrac{1}{2}(u^*(\gamma_i +1)-u_R(\gamma_i -1)+2c_R). \end{equation}

This provides the complete structure of the wave pattern as a function of space and time for a Riemann problem where the fluids on either side of the initial discontinuity may be different but where either the ideal or stiffened gas equations of state apply.

3. Construction of the solution for the reflection of an impulsive shock

Now that the solution to the Riemann problem has been derived for a multi-fluid problem, we can construct the complete solution for an idealised shock pulse interacting with a liquid–gas interface. Two pulse types are considered, with the second being an extension of the first. An illustration of the two pulses is given in figure 2. The first pulse type is a positive square pulse, denoted a P-type pulse, where a region of shocked fluid is bounded by a shock and a rarefaction. The initial location of the shock front is $x=0$, and the tail is at $x=x_T$.

Figure 2. Illustration of the two idealised pulses: (a) P-type pulse, (b) PN-type pulse. The dash-dotted vertical line indicates the initial location of the interface.

The second pulse consists of both positive and negative regions, and is denoted a PN-type pulse. This is motivated by practical shock pulses in liquids that consist of a shock front followed by an exponentially decaying pressure, with the pressure going negative before returning to ambient pressure. An example of this is a lithotripter wave (Church 1989; Johnsen & Colonius 2008). The reflection/transmission problem for such a wave can again be treated analytically by considering the idealised version of this pulse. This consists of a shock front at $x=0$, an (initially) infinitely thin rarefaction wave ($x=x_M$) and a shock wave at the tail ($x=x_T$) separating the ambient liquid from the rarefied liquid. In this work, we derive the solution for a pulse where the positive and negative parts can have different widths and pressures.

3.1. P-type pulse

The P-type pulse consists of a shock wave and an infinitely thin rarefaction wave separated by some distance. This distance determines the pulse width $\varDelta _\mathcal {S}$. The region in the middle is said to be shocked (denoted by the subscript $\mathcal {S}$), and the variables in this region can be determined via the Rankine–Hugoniot relations provided that either the shock pressure or Mach number is given.

The solution to the reflection/transmission problem comprises three consecutive Riemann problems, which are shown in figure 3 and overviewed below. We begin with an initial state with a liquid on the left and a gas on the right, separated by an interface at $x=0$. Part of the liquid is shocked, and we begin with the shock front arriving at the liquid–gas interface.

1. (i) $t=t_1$. The right-going shock reaches the liquid–gas interface. This is the first Riemann problem, with the left state being $\boldsymbol {w}_\mathcal {S}$, and the right state $\boldsymbol {w}_{g,0}$. The solution to this Riemann problem is a left-going rarefaction wave and a right-going shock wave that propagates into the gas phase. We now have two new states corresponding to the rarefied liquid $\boldsymbol {w}^*_{L,1}$ and the shocked gas $\boldsymbol {w}^*_{R,1}$.

2. (ii) $t=t_2$. The left-going rarefaction wave reaches the right-going rarefaction at the tail of the pulse, creating the second Riemann problem, with left state $\boldsymbol {w}_{l,0}$ and right state $\boldsymbol {w}^*_{L,1}$. The solution to the second Riemann problem is two rarefaction waves, with the pressure in the star region being lower than the pressures in the left and right states. This creates a single rarefied state $\boldsymbol {w}^*_2$.

3. (iii) $t=t_3$. The right-going rarefaction wave from the second Riemann problem reaches the fluid interface, which has been moving rightwards at $u_1^*$ for $t>t_1$. This creates the third Riemann problem, with left state $\boldsymbol {w}^*_2$ and right state $\boldsymbol {w}^*_{R,1}$. The solution to the final Riemann problem is a shock wave travelling left into the rarefied liquid, and a rarefaction wave travelling right into the shocked gas. Two additional states are created, corresponding to the shocked liquid $\boldsymbol {w}^*_{L,3}$ and the rarefied gas $\boldsymbol {w}^*_{R,3}$.

Figure 3. Illustration of the stages of the solution for the reflection/transmission of a P-type pulse. The Riemann problems occurring at $t=t_1, t_2, t_3$ are shown on the right-hand side, with their left and right initial states. Regions of shocked and rarefied fluid are shown in red and blue, respectively. Shock and rarefaction waves are denoted $s$ and $r$, respectively.

For a fast–slow interface such as water–air, the solution structure shown in figure 3 remains the same irrespective of the shock strength. This is shown in Appendix A for a wide range of shock pressures.

The time at which each Riemann problem takes place is a function of the Mach number and the width of the shock, and can be used to compute the widths of the transmitted and reflected pulses. Before considering this, we must first acknowledge a key assumption made in this work. In practice, a square shock will not remain square because the rarefaction wave spreads out over time, with its front travelling faster than its tail. In order to have correctly defined Riemann problems, an assumption is made that, over certain time scales, the rarefaction wave remains infinitely thin. The validity of this assumption is central to the analytical approach derived in this work, so is considered in more detail here. Suppose that we have a right-going rarefaction wave connecting a lower-pressure state $\boldsymbol {w}_L$ to a high-pressure right state $\boldsymbol {w}_R$. The rate at which this wave spreads out can be obtained by considering the speeds of the front and tail of the rarefaction wave. Defining the width of the rarefaction wave as $\varDelta _r$, the time rate of change of width is

(3.1)\begin{equation} \frac{{\rm d}}{{\rm d} t}(\varDelta_r) = \varDelta_r' = \zeta_{r,F}-\zeta_{r,T} = \frac{1}{2}(\gamma_i+1)(u_R - u_L). \end{equation}

In the case of a rarefaction connecting a quiescent fluid to a shocked fluid (as in the case of a square shock), we have, using the Rankine–Hugoniot relations,

(3.2)\begin{equation} \varDelta_r' = c_0\left(\frac{M^2-1}{M}\right). \end{equation}

Clearly, $\varDelta _r' \rightarrow 0$ as $M\rightarrow 1$, so the rarefaction will remain infinitely thin in this limit. Figure 4 shows how this changes as we move away from the limit. As the Mach number increases, the rate of spread of the wave increases, but because of the stiff nature of water, the shock pressure can rise quite high before this becomes significant. This is why the method set out in this work can readily be applied only to fast–slow problems and not slow–fast problems.

Figure 4. Time rate of change of width of a rarefaction wave moving into a region of shocked water.

This analysis is used to make the following assumptions in the solution process outlined in figure 3. First, the right-going rarefaction wave at the tail of the square wave remains infinitely thin over $t_1 < t \leq t_2$. Second, the rarefaction created at time $t=t_1$ is infinitely thin over $t_1< t \leq t_2$, and finally, the right-going rarefaction created at $t_2$ remains infinitely thin over $t_2< t\leq t_3$. These assumptions are required only when a rarefaction comes into contact with another wave or discontinuity. Thus for the left-going rarefaction created at $t=t_2$ and the right-going rarefaction created at $t_3$, we do not make this assumption, and the wave structures are computed as per § 2.2. It is noted that Grove & Menikoff (1990) describe the rarefaction wave in the case of the fast–slow interaction of a shock wave in water with a gas bubble as a ‘rarefaction shock’ on account of how thin it is, so these assumptions do have some physical basis. The validity of these assumptions is tested later in this work by comparing with high-order-accurate numerical simulations.

We can now derive the times and locations at which each Riemann problem will take place, which subsequently allows for the widths of the transmitted and reflected pulses to be obtained. The time and location of the second Riemann problem can be computed exactly based on the wave speeds derived from the solution to the first Riemann problem. Note that the assumption during this time period is that the rarefaction does not spread out. The time $t=t_2$ is

(3.3)\begin{equation} t_2 = \frac{-x_T}{2c_\mathcal{S}}, \end{equation}

and the location is

(3.4)\begin{equation} x_2 = \frac{(c_\mathcal{S} - u_\mathcal{S})x_T}{2c_\mathcal{S}}. \end{equation}

The third Riemann problem takes place when the right-going rarefaction reaches the interface:

(3.5)\begin{equation} t_3 = x_T\biggl(\frac{u_s - c_\mathcal{S} - u_1^* - c_{L,1}^*}{2c_\mathcal{S} c_{L,1}^*}\biggr). \end{equation}

The location $x_3$ is then the position of the interface at this time:

(3.6)\begin{equation} x_3 = x_{lg}(t_3)= u_1^*t_3. \end{equation}

This time marks the end of the interaction, leading to a transmitted pulse and a reflected pulse. The width of the reflected pulse at $t=t_3$ is the distance between the contact and the left-going rarefaction wave created at time $t=t_2$. Therefore,

(3.7)\begin{align} \varDelta_{\mathcal{R}} & = x_3 - (x_2+(u_{l,0}-c_{l,0})(t_3-t_2)) \nonumber\\ & ={-}x_T\biggl(\frac{(u_1^* + c_{L,1}^* - u_{l,0} + c_{l,0})(u_1^* + c_\mathcal{S} - u_\mathcal{S})}{2c_{L,1}^*c_\mathcal{S}}\biggr)\geq \varDelta_\mathcal{S}. \end{align}

The width of the reflected pulse is therefore at least the width of the incident shock, with the width increasing as the shock pressure increases due to the higher speed of the interface following the first Riemann problem. The shocked gas is bounded by the portion of the incident shock that was transmitted into the gas at $t=t_1$ and the rarefaction wave created at $t=t_3$. For times $t>t_1$, the position of the right-going shock wave is $x_{\mathcal {S},g}(t) = \zeta _{\mathcal {S},g}t$. Therefore, the width of the shocked region in the gas phase at time $t=t_3$ is

(3.8)\begin{equation} \varDelta_{\mathcal{T}} = \zeta_{\mathcal{S},g}t_3 - x_3 = x_T \biggl(\frac{(u_1^* + c_{L,1}^* - u_\mathcal{S} + c_\mathcal{S}) (u_1^*-\zeta_{\mathcal{S},g})}{2c_{L,1}^*c_\mathcal{S}}\biggr). \end{equation}

3.2. PN-type pulse

The solution strategy for this problem is similar to that of the P-type pulse, with Riemann problems being solved each time two waves interact or a wave interacts with the interface. The solution structure is shown in figure 5 and consists of six Riemann problems. The same assumptions are made that rarefaction waves remain infinitely thin during the intermediate stages. This applies only to waves that will subsequently interact with another wave.

1. (i) $t=t_1$. The right-going shock reaches the liquid–gas interface. This is the first Riemann problem, with the left state being $\boldsymbol {w}_{\mathcal {S},+ve}$, and the right state $\boldsymbol {w}_{g,0}$. As with the positive square wave, the solution consists of a left-going rarefaction and a right-going shock. The interface will move with the velocity determined from this Riemann problem.

2. (ii) $t=t_2$. The left-going rarefaction created at $t=t_1$ reaches the right-going rarefaction that separates the positive and negative parts of the incident pulse. The solution is two rarefaction waves with a constant pressure between them.

3. (iii) $t=t_3$. The left-going rarefaction created at $t_2$ reaches the shock wave at the tail of the pulse. The solution is a left-going rarefaction and a right-going shock. Strictly speaking, this leads to two states, but as will be shown, the density across the discontinuity is extremely small, so the discontinuity may be considered degenerate.

4. (iv) $t=t_4$. The right-going rarefaction created at $t_2$ reaches the interface. If the widths of the two parts of the pulse are similar, then this will likely take place after $t_3$, but this order does not affect the solution. The solution is a left-going shock and a right-going rarefaction. The velocity in the star region is now negative, so the interface will move leftwards.

5. (v) $t=t_5$. The right-going shock created at $t_3$ and the left-going shock created at $t_4$ meet, creating the fifth Riemann problem. The solution is two new shocks propagating in opposite directions: one going left, and one going right, towards the interface.

6. (vi) $t=t_6$. The right-going shock created at $t_5$ reaches the interface, which has been travelling at $u_{L,4}^*$ for $t>t_4$. The solution to this Riemann problem is a left-going rarefaction and a right-going shock.

Figure 5. Illustration of the stages of the solution for the reflection/transmission of a PN-type pulse. The Riemann problems occurring at $t=t_1,\ldots,t_6$ are shown on the right-hand side, with their left and right initial states. Regions of shocked and rarefied fluid are shown in red and blue, respectively. Shock and rarefaction waves are denoted $s$ and $r$, respectively.

As with the P-type pulse, it can be shown that this solution structure is always the same for a water–air interface, but only for certain pressures. In Appendix B, this is shown for pulses where the magnitude of the pressure in both the negative and positive parts of the pulse is up to 500 MPa. Beyond this, it will be shown in § 4.5 that a vacuum state can be created to the right of the interface, and the present model cannot accommodate such a case. However, this depends on the fluids considered and the value of $p_\infty$. Clearly, a higher value of $p_\infty$ will accommodate lower pressures in the liquid, but this may lead to the gas phase being put into vacuum earlier. Depending on the widths of the two parts of the pulse, $t_4$ may occur before or after $t_3$. For example, if the positive part is very thin compared to the negative part, then it is likely that $t_4$ will occur first. However, Appendix B shows that this will not alter the overall solution. The times and locations of each Riemann problem can again be derived by considering the speeds at which the different waves and the interface travel during the process. These are also given in Appendix B, and can be used to derive the widths of the reflected and transmitted pulses. The reflected pulse is bounded by the rarefaction wave created at time $t=t_3$ and the rarefaction wave created at $t=t_6$. This leads to the expression

(3.9)\begin{equation} \varDelta_\mathcal{R} = x_6 - x_3 - (t_6 - t_3)(u_{l,0}-c_{l,0}). \end{equation}

The transmitted pulse is bound by the shock front that is transmitted into the gas at $t=t_1$ and the shock created at $t_6$. Therefore,

(3.10)\begin{equation} \varDelta_\mathcal{T} = \zeta_{R,1}t_6 - x_6. \end{equation}

The widths of the parts are discussed in subsequent sections in the context of the different regimes: the acoustic limit, weak shocks and strong shocks.

4. Interaction of an impulsive shock with a water–air interface

4.1. Comparison with numerical simulations

Thus far, a description of the process of reflection and transmission of two idealised pulses at a liquid–gas interface has been presented. In this subsection, the analytical model is compared to numerical simulations to confirm that it is indeed able to describe the reflection and transmission process. To this end, we employ a well-established front-tracking framework, FronTier (Glimm et al. 1981, 1998; Chern et al. 1986; Du et al. 2006; Bempedelis & Ventikos 2020b), that solves the compressible Euler equations closed with the SG EoS (see § 2.1). The equations are solved on a uniform grid using a fifth-order WENO scheme and a third-order Runge–Kutta method. In the numerical results presented, a grid convergence study (see Appendix C) has confirmed that numerical diffusion effects are negligible. The employed front-tracking solver has been validated against several experimental and numerical data, as well as benchmark tests, in a wide range of multiphase flow problems involving shock waves and interfaces (Bempedelis & Ventikos 2020a,b, 2022; Bempedelis et al. 2021).

Suppose that we have a water–air interface initially at $x=0\,\text {m}$. The primitive variables in the ambient liquid and gas are

(4.1a,b)\begin{equation} \boldsymbol{w}_{l,0} = \begin{pmatrix} \rho_{l,0} \\ u_{l,0} \\ p_{l,0} \end{pmatrix} = \begin{pmatrix} 998 \\ 0.0 \\ 1\times10^5 \end{pmatrix},\quad \boldsymbol{w}_{g,0} = \begin{pmatrix} \rho_{g,0} \\ u_{g,0} \\ p_{g,0} \end{pmatrix} = \begin{pmatrix} 1.16 \\ 0.0 \\ 1\times10^5 \end{pmatrix}. \end{equation}

At time $t=t_1$, a right-going square pulse with width $\varDelta _\mathcal {S}=2\times 10^{-4}\,\text {m}$ comes into contact with the interface. For both pulse types considered, the pressure in the positive part is $p_\mathcal {S} = p_{\mathcal {S},+ve} = 10^7$ Pa. For the PN-type pulse, $p_{\mathcal {S},-ve}-p_{l,0} = -(p_{\mathcal {S},+ve} - p_{l,0})$ and $\varDelta _{\mathcal {S},+ve}=\varDelta _{\mathcal {S},-ve}=1\times 10^{-4}\,\text {m}$. The velocity and density in the pulses are determined using the Rankine–Hugoniot relations. Figures 6 and 7 show the solutions for both cases. The initial conditions are shown along with intermediate states, and at time $t=2\times 10^{-7}\,\text {s}$, which is shortly after the process completes in both cases. For the P-type pulse, all intermediate conditions are shown, and the states are labelled to show that the solution does indeed follow the process set out in § 3. For the PN-type pulse, a single intermediate condition is shown for brevity. The agreement between the analytical and numerical predictions is excellent at all stages, showing that the pulse–interface interaction process can indeed be approximated via a small number of wave interactions.

Figure 6. Analytical and numerical predictions for the reflection/transmission of a P-type pulse ($p_\mathcal {S} = 10\,\text {MPa}$) at a water–air interface: (a) initial conditions; (b) $t=5\times 10^{-8}\,\text {s}$ ($t_1< t< t_2$); (c) $t=1\times 10^{-7}\,\text {s}$ ($t_2< t< t_3$); (d) $t=2\times 10^{-7}\,\text {s}$ ($t_3< t$).

Figure 7. Analytical and numerical predictions for the reflection/transmission of a PN-type pulse ($\kern0.7pt p_{\mathcal {S},+ve} = 10\,\text {MPa}$) at a water–air interface: (a) initial conditions; (b) $t=8\times 10^{-8}\,\text {s}$ ($t_4 < t < t_5$); (c) $t=2\times 10^{-7}\,\text {s}$ ($t_6< t$).

For the P-type pulse, the reflected pulse puts the liquid into tension. It can also be seen that this happens only when the rarefaction wave meets the tail of the incident pulse, corresponding to the second Riemann problem. Both the numerical and analytical models predict a slightly lower pressure magnitude in the reflected pulse than is predicted by linear theory. Linear theory predicts pressure $-9.89$ MPa, whereas the numerical and analytical approaches lead to reflected pressure $-9.71$ MPa. Whilst this difference is relatively small, it will be shown later that the present theory diverges significantly from linear theory for larger shock pressures. This also means that the pressure in the reflected pulse is not simply a sum of the strength of the left-going rarefaction wave created at $t_1$ and the ambient pressure. The reflected pulse is slightly wider than the incident pulse, and the interface has moved slightly rightwards. This pulse is considered to be ‘weak’, with $u_\mathcal {S} \ll c_{l,0}$, and it can be seen from the numerical results that the rarefaction waves do remain extremely thin over the relevant time scales. This suggests that for weak shock waves at least, the assumptions made in order to derive the analytical solution are valid.

For the PN-type pulse, a reflection of the wave can be seen and the magnitude of the negative pressure is the same as for the P-type pulse. The region in tension is slightly wider that the initial shocked region, again giving the same result as for the P-type pulse. The negative part of the incident pulse reflects as a positive part, and the pressure in this region is higher than would be predicted by linear theory, and also higher than the incident shock pressure. Linear theory predicts reflected pressure 9.89 MPa, whereas both the analytical and numerical models predict 10.09 MPa. Furthermore, while the tension region in the reflected pulse is wider than the incident shocked region, the shocked region in the reflected pulse is slightly narrower. Figure 7(b) also shows that a region with pressure significantly below that of the incident pulse exists in the intermediate state. Referring to figure 5, this is $\boldsymbol {w}_2^*$, and it will exist only temporarily ($t_2< t< t_5$). It comes about when the rarefaction wave created at $t_1$ reaches the rarefaction wave that separates the positive and negative parts of the pulse. The pressure in this state will always be more negative than the pressure in either the incident or reflected parts of the pulse, and should be considered carefully in practical applications as this could lead to significantly higher levels of cavitation depending on the duration for which it exists.

These results show that the proposed analytical approach is valid, with the results closely matching those from high-order numerical simulations of the Euler equations. It is also clear that there are differences between these results and the predictions of linear theory. Given that linear theory allows for the derivation of reflection and transmission coefficients, we will now extend these concepts to the nonlinear problems considered in this work.

4.2. Reflection and transmission coefficients

In the context of this work, the simplest definition of a reflection and transmission coefficient is the ratio of the peak pressure in the reflected and transmitted waves relative to the incident shock pressure. For the P-type pulse, these coefficients are defined as

(4.2)$$\begin{gather} C_\mathcal{R} = \frac{p_2^*-p_{l,0}}{p_\mathcal{S}-p_{l,0}}, \end{gather}$$

(4.3)$$\begin{gather}C_\mathcal{T} = \frac{p_1^*-p_{g,0}}{p_\mathcal{S}-p_{l,0}}. \end{gather}$$

There are two comparisons that we can make here: first with linear theory, and second with a solution derived from a single Riemann problem solved at the interface. As discussed in the Introduction, linear theory results in reflection and transmission coefficients that depend only on the differences in impedance between the media across the interface. An alternative approach is to solve a Riemann problem only at the interface to determine the strength of the right-going shock and the left-going rarefaction. This is what would be obtained if the approach of Henderson (1989) were to be applied directly to a finite-width pulse. In this case, the pressure in the reflected pulse is the sum of the strength of the rarefaction wave and the ambient pressure in the liquid. Figure 8 shows the reflection and transmission coefficients for a range of shock pressures using the three different approaches. Numerical solutions are also included for comparison. The results are presented for the peak over-pressure $p'_\mathcal {S} = p_\mathcal {S} - p_{l,0}$.

Figure 8. Reflection and transmission coefficients based on peak pressures for a P-type pulse.

First, it is clear that in the acoustic limit, all methods tend to the same value. Away from this, significant differences are observed. For the reflection coefficient, the current analytical approach and the numerical simulations predict a decrease in the magnitude of the reflection coefficient as the shock pressure increases. This represents a divergence from linear theory as one might expect for what is a highly nonlinear problem. The reflection coefficient derived by solving a single Riemann problem also differs significantly from the present approach, suggesting that the pulse–interface interaction problem cannot be modelled accurately by only considering what happens at the interface. This is because the assumption that wave interactions do not alter the wave properties is not valid. Thus we must consider how the rarefaction wave is altered as it propagates into the region of lower pressure behind the shock. For the transmission coefficient, there are no further interactions, so solving a single Riemann problem at the interface is sufficient. However, it can again be seen that linear theory fails at large shock pressures. The numerical results agree very well with the present analytical approach, although some divergence is seen in the transmission coefficient at higher pressures. This is likely due to the assumptions in solving the interface interaction, i.e. the coupling of the two fluids at the interface, in the numerical methods that we use, that include linearisations, isentropic extrapolations and ‘overheating’ effects (for more details, see Fedkiw, Marquina & Merriman 1999; Hu & Khoo 2004).

For a pulse with both positive and negative parts, it was shown in the previous section that the two parts reflect differently. It therefore makes sense to consider reflection and transmission coefficients for the two parts separately to understand how and why these vary. The coefficients are defined using specific solutions to the Riemann problems that lead to the reflected and transmitted pulses

(4.4)$$\begin{gather} C_{\mathcal{R}_{{+}ve}} = \frac{p_3^* - p_{l,0}}{p_{\mathcal{S},+ve} - p_{l,0}}, \end{gather}$$

(4.5)$$\begin{gather}C_{\mathcal{R}_{{-}ve}} = \frac{p_5^* - p_{l,0}}{p_{\mathcal{S},-ve} - p_{l,0}}, \end{gather}$$

(4.6)$$\begin{gather}C_{\mathcal{T}_{{+}ve}} = \frac{p_1^* - p_{g,0}}{p_{\mathcal{S},+ve} - p_{l,0}}, \end{gather}$$

(4.7)$$\begin{gather}C_{\mathcal{T}_{{-}ve}} = \frac{p_1^* - p_{g,0}}{p_{\mathcal{S},-ve} - p_{l,0}}. \end{gather}$$

Thus $C_{\mathcal {R},+ve}$ refers to the reflection of the positive part of the incident pulse. We first consider the case where the magnitudes of the over- and under-pressures in the incident pulse are the same. The coefficients are shown in figure 9, along with numerical results. Agreement between the analytical and numerical results is again very good. The reflection of the positive part of the incident wave is the same as for the P-type pulse. The magnitude of the reflection of the negative part increases with shock pressure, meaning that the maximum pressure of the reflected pulse will be greater than the maximum pressure in the incident pulse. This behaviour is mirrored in the transmission coefficients, and again the transmission of the positive part of the pulse is the same as for the positive square wave.

Figure 9. Reflection and transmission coefficients based on peak pressures for a PN-type pulse where the widths of the two incident regions are constant and $p'_{\mathcal {S},-ve} = -p'_{\mathcal {S},+ve}$. Numerical results are shown with crosses.

Coefficients (4.4)–(4.7) assume that the reflection and transmission of each part of the pulse depend only on the magnitude of that particular part. In many practical applications, the magnitude of the negative part is different from that of the positive part. It is therefore worth considering how these coefficients change when the peak pressure in the negative part is different from that in the positive part. This is shown in figure 10. This tells us two things. First, the reflection of the positive part is independent of the magnitude of the negative part. This has already partly been shown because the reflection coefficient for a P-type wave was the same as that derived for the positive part of the PN-type wave. However, the reflection coefficient of the negative part is dependent on the strength of both the positive and negative parts of the incident pulse.

Figure 10. Reflection coefficients for a PN-type pulse where the negative and positive pressure magnitudes differ.

For the PN-type pulse, a region of strongly rarefied liquid ($\boldsymbol {w}_2^*$) results when the left-going rarefaction wave created at $t=t_1$ reaches the tail of the positive part of the incident pulse ($t=t_2$). This state exists only for $t_2< t< t_5$. As with the negative part of the reflected pulse, the magnitude of this negative pressure is also a decreasing function in the incident shock pressure. By non-dimensionalising by the difference between the positive and negative pressures in the incident wave so that

(4.8)\begin{equation} C_{\mathcal{R},{tmp}}=\frac{p_2^* - p_{l,0}}{p_{S,+ve} - p_{S,-ve}}, \end{equation}

the pressures in this temporary region all collapse onto a single line for weaker shocks, as shown in figure 11. For shocks where $p_{\mathcal {S},+ve}>100$ MPa, small deviations are seen, with the magnitude of the coefficient reducing for pulses where the relative magnitude of the negative pressure is greater. This pressure will exist only for a time determined by the width of the negative part of the incident pulse. Depending on the width and the minimum pressure that exists during this time, this may be of considerable importance as it could lead to higher levels of cavitation than would be predicted by only considering the final pressures in the reflected pulse.

Figure 11. Coefficient for the minimum pressure that will exist during the reflection process, as defined by $C_{\mathcal {R},{tmp}}$.

The solution (§ 3) and reflection/transmission coefficients have been derived for an arbitrary set of liquid–gas media modelled with the SG EoS, but their exact values depend on the EoS parameters. In this work, we consider the values $\gamma =2.955$ and $p_\infty =722\,\text {MPa}$ for liquid water, but different values are used in other works, depending on factors such as the temperature. For example, Haller et al. (2003) use $\gamma =3$, $p_\infty =618$ MPa, and Métayer, Massoni & Saurel (2004) use $\gamma =2.35$, $p_\infty =1000$ MPa. Irrespective of the particular liquid under consideration, the values of $p_\infty$ do have an effect on the reflection and transmission of the pulse. This is exemplified in figure 12 for a P-type pulse, suggesting two things. First, the level of ‘stiffness’ of the water influences the results significantly only for higher pressures, which is expected as compressibility plays a greater role here. Second, the ‘stiffer’ the water is assumed to be, the smaller the change in the reflection coefficient as the shock pressure increases. The change in the transmission coefficients can be explained more readily using either linear theory or an approximate Riemann solver (e.g. Hu & Khoo 2004) because the change in $p_\infty$ alters the speed of sound, which in turn changes the impedance.

Figure 12. (a) Reflection and (b) transmission coefficients for a P-type pulse at a water–air interface with varying values of $p_\infty$, for $\gamma _l=2.955$ and $\gamma _g=1.4$.

As well as considering the pressure in the reflected and transmitted pulses, it is also of interest to consider the velocity. This is shown in figure 13 for both pulse types following the solution to each Riemann problem. For the PN-type pulse, $p'_{\mathcal {S},-ve} = p'_{\mathcal {S},+ve}$ in this example. The particle velocities are remarkably constant as a function of the particle velocity in the incident pulse. Following the initial interaction of the shock front with the interface, the interface will move rightwards at $u^*_1 \approx 2u_\mathcal {S}$ for all shock pressures. It is interesting to note that this is the same as the particle velocity obtained for a linear problem by considering the limit as the difference in impedance becomes very large. This relationship was also noted by Grove & Menikoff (1990) when considering a similar fast–slow problem for a shock wave. For the P-type pulse, $u_1^*$ is the velocity in the transmitted pulse, and $u_2^*$ is the velocity in the reflected pulse. Also, $u_3^*$ is the velocity of the fluid in between the reflected and transmitted pulses, and it can be seen that this is zero for all shock pressures, again demonstrating that the entire process results in only a reflected and transmitted pulse with no further changes to either fluid. For the PN-type pulse, as with the P-type pulse, the interface will initially move at $u_1^* \approx 2u_\mathcal {S}$, this time over the interval $t_1 < t < t_4$. The solution to the fourth Riemann problem then reverses the interface, with the magnitude being almost exactly the same for all but the largest shocks. Finally, the interface is stopped at $t_6$, signifying the end of the process.

Figure 13. Particle velocity in the star regions following each Riemann problem: (a) P-type pulse; (b) PN-type pulse.

4.3. The acoustic limit

To provide further insights into the properties of the reflected and transmitted waves, it is useful to consider the problem in three different regimes: the acoustic limit, the weak shock regime and the strong shock regime. Formally, the acoustic limit is defined by $M\rightarrow 1^+$. The particle velocities tend to zero, and the wave speeds all tend to the sound speed in the relevant fluid. It has already been shown that any rarefaction waves will remain infinitely thin in the acoustic limit. Since this is the only assumption made in the derivation in § 3, the method is exact in this limit.

The reflection and transmission coefficients were all shown to tend to the values predicted by linear theory in figures 8, 9 and 10. Furthermore, we also have that

(4.9)\begin{equation} \lim_{M\rightarrow 1^+} u_\mathcal{S} = 0, \end{equation}

so $u_i^* \rightarrow 0$ for both pulses. The sound speeds and wave speeds will all tend to the ambient speeds in the respective fluids. For a P-type pulse, the width of the reflected wave, $\varDelta _{\mathcal {R}}$, is obtained by taking the limit of (3.7):

(4.10)\begin{equation} \lim_{M\to 1^+} \varDelta_\mathcal{R} ={-}x_T = \varDelta_\mathcal{S} \end{equation}

and

(4.11)\begin{equation} \lim_{M\to 1^+} t_3 ={-}x_T/c_{l,0}. \end{equation}

In other words, the width of the reflected pulse is exactly the same as the width of the incident pulse. The width of the transmitted pulse can be found by taking the limit of (3.8):

(4.12)\begin{equation} \lim_{M\to 1^+} \varDelta_\mathcal{T} ={-}x_T\,\frac{c_{g,0}}{c_{l,0}}. \end{equation}

That is, the width of the transmitted pulse is equal to the width of the incident pulse multiplied by the ratio of the sound speeds of the two fluids. This is analogous to the linear case, where the wavelength reduces proportional to the ratio of sound speeds of the two fluids.

Using the same argument for a PN-type pulse,

(4.13)\begin{equation} \lim_{M\to 1^+} x_6 = 0, \end{equation}

so the interface will also remain in the same place in the acoustic limit. The reflected and transmitted pulse widths are found by taking the limit of (3.9) and (3.10) as $M\rightarrow 1^+$, and it is straightforward to show that these lead to the same expressions as for the P-type pulse.

4.4. Weak shocks

Moving away from the acoustic limit, the assumptions that the particle velocities are zero and the Mach number is $M = 1$ are no longer valid. This represents a departure from where linear theory is formally valid. For a weak shock, it is assumed that $u_\mathcal {S}$ remains small compared to $c_{l,0}$. Therefore, the rarefaction waves will remain very thin over the relevant time scales. For the P-type pulse, the negative pressure in the reflected pulse comes about following the second Riemann problem. For weak shocks, the pressure behind the rarefaction wave produced by the initial interaction of the shock and the interface remains roughly constant at levels slightly above the ambient pressure. Therefore, the solution to the second Riemann problem is determined primarily by the difference in the particle velocities. The increasing positive velocity in the right initial state as the incident shock pressure increases leads to a rarefied region with an increasingly negative pressure.

For weak shocks, because the solution to the first Riemann problem yields a pressure that is very close to the ambient pressure, it is possible to derive a closed-form approximate expression for the reflected and transmitted wave using the solution derived in § 3. The interface moves with $u_1^* \approx 2u_S$ for $t_1 \leq t \leq t_3$, corresponding to the liquid becoming rarefied and the gas becoming compressed. For a weak incident shock, $p_1^*$ is only slightly larger than $p_{l,0}$, so $c_{L,1}^* \approx c_{l,0}$. Therefore, using these two approximations, we can rewrite $x_3$ as

(4.14)\begin{equation} x_3 \approx{-}u_\mathcal{S} x_T \left(\frac{c_{l,0} + c_\mathcal{S} + u_\mathcal{S}}{c_\mathcal{S} c_{l,0}}\right). \end{equation}

Using the same assumptions and assuming that the ambient liquid is at rest, the width of the rarefied region at $t=t_3$ is

(4.15)\begin{equation} \varDelta_\mathcal{R} \approx{-}x_T \left(\frac{(c_{l,0}+u_\mathcal{S}) (c_\mathcal{S}+u_\mathcal{S})}{c_{l,0}c_\mathcal{S}}\right). \end{equation}

We now call upon the conservation of mass in the liquid. Comparing the initial state and that at $t=t_3$, we have

(4.16)\begin{equation} \int_{-\infty}^{x_T}\rho_{l,0}\,{{\rm d}\kern0.7pt x} + \int_{x_T}^{0}\rho_\mathcal{S}\,{\rm d}x = \int_{-\infty}^{x_3 - \varDelta_\mathcal{R}} \rho_{l,0}\,{{\rm d}\kern0.7pt x} + \int_{x_3-\varDelta_\mathcal{R}}^{x_3} \rho_2^*\,{{\rm d}\kern0.7pt x}. \end{equation}

Evaluating this and rearranging gives

(4.17)\begin{equation} \rho_2^* = \frac{\rho_{l,0}(x_T + \varDelta_\mathcal{R}-x_3) -\rho_\mathcal{S} x_T}{\varDelta_\mathcal{R}} = \frac{c_{l,0} c_\mathcal{S} \rho_\mathcal{S}}{(c_0+u_\mathcal{S})(c_\mathcal{S}+u_\mathcal{S})}. \end{equation}

The pressure in the rarefied region can be obtained by rearranging the Rankine–Hugoniot relations (2.5) and (2.6):

(4.18)\begin{equation} p_2^* = \frac{(p_{l,0}+p_\infty)(1+\gamma_l)}{1+\gamma_l \left(\dfrac{4\rho_{l,0}}{\rho_{l,0}(1-\gamma_l)+\rho_2^*(1+\gamma_l)}-1\right)} - p_\infty. \end{equation}

Similarly, in the gas we have that

(4.19)\begin{equation} \int_0^\infty \rho_{g,0}\,{{\rm d}\kern0.7pt x} = \int_{x_3}^{x_3+\varDelta_\mathcal{T}} \rho_{R,1}^*\,{{\rm d}\kern0.7pt x} + \int_{x_3+\varDelta_\mathcal{T}}^\infty p_{g,0}\,{{\rm d}\kern0.7pt x}. \end{equation}

The shock wave that propagates into the gas is very weak in comparison to that in the liquid, so $\zeta _{\mathcal {S},g} \approx c_{g,0}$. Using this, we can obtain the following expression for $\rho _{R,1}^*$, which again depends only on the initial data:

(4.20)\begin{equation} \rho_{R,1}^* = \frac{c_{g,0}\rho_{g,0}}{c_{g,0}-2u_\mathcal{S}}. \end{equation}

The pressure in the shocked gas can again be obtained by way of the Rankine–Hugoniot relations. These expressions enable the reflection–transmission of a shock wave at a liquid–gas interface to be computed without solving a Riemann problem, potentially allowing for the approach to be more easily integrated into other applications.

This also has a useful physical interpretation, explaining how a region of rarefied liquid (where $p< p_{l,0}$ and $\rho <\rho _{l,0}$) can result from a region of shocked liquid. For a solution consisting of only a left-going rarefied region (in the liquid), the only way for mass in the liquid to be conserved is for the interface to move rightwards to precisely compensate for the difference in the mass of the incident and reflected pulses. Should this not happen (e.g. if the interface velocity was altered by some other mechanism), then additional waves would be created in order to balance out the mass in the liquid. However, in the case of a weak shock, these differences are balanced by the movement of the interface. For the PN-type pulse, the interface initially moves rightwards before moving leftwards for $t_4 < t < t_6$. This explains why the negative part of the reflected pulse is wider, and the positive part is narrower. Noting these changes in the widths of the pulse, it is also useful to consider the total ‘strength’ of the reflected and transmitted pulses by integrating the pressure over their respective widths. For the P-type pulse, modified reflection and transmission coefficients can therefore be defined as

(4.21)$$\begin{gather} \hat{C}_\mathcal{R} = \frac{\displaystyle\int_{\varDelta_\mathcal{R}}(p - p_{l,0})\,{{\rm d}\kern0.7pt x}}{ \displaystyle\int_{\varDelta_\mathcal{S}}(p_\mathcal{S}-p_{l,0})\,{{\rm d}\kern0.7pt x}} , \end{gather}$$

(4.22)$$\begin{gather}\hat{C}_\mathcal{T} = \frac{\displaystyle\int_{\varDelta_\mathcal{T}}(p - p_{g,0})\,{{\rm d}\kern0.7pt x}}{ \displaystyle\int_{\varDelta_\mathcal{S}}(p_\mathcal{S}-p_{l,0})\,{{\rm d}\kern0.7pt x}} . \end{gather}$$

For the PN-type pulse, we again consider the reflection of the positive and negative parts separately:

(4.23)$$\begin{gather} \hat{C}_{\mathcal{R},+ve} = \frac{\displaystyle\int_{\varDelta_{\mathcal{R},-ve}}(p - p_{l,0})\, {{\rm d}\kern0.7pt x}}{\displaystyle\int_{\varDelta_{\mathcal{S},+ve}}(p_{\mathcal{S},+ve}-p_{l,0})\,{{\rm d}\kern0.7pt x}}, \end{gather}$$

(4.24)$$\begin{gather}\hat{C}_{\mathcal{R},-ve} = \frac{\displaystyle\int_{\varDelta_{\mathcal{R},+ve}}(p - p_{l,0})\,{{\rm d}\kern0.7pt x}}{\displaystyle\int_{\varDelta_{\mathcal{S},-ve}}(p_{\mathcal{S},-ve}-p_{l,0})\,{{\rm d}\kern0.7pt x}}. \end{gather}$$

Similar expressions can be written for the modified transmission coefficients. For the P-type pulse, the modified reflection coefficient is less sensitive to the incident shock pressure than the coefficient presented earlier. This is shown in figure 14(a) and is because the drop in the pressure magnitude in the reflected pulse is partially offset by the increase in its width. A similar pattern is observed for the PN-type pulse. The wider reflected pulse offsets the reduction in the pressure magnitude, although by less than for the P-type pulse. In the weak regime where the rarefaction waves are assumed to not spread out, the modified reflection and transmission coefficients for the P-type pulse are independent of the incident pulse width. For the PN-type pulse, this is true only when $\varDelta _{\mathcal {S},+ve}/\varDelta _{\mathcal {S},-ve}$ is constant. If the relative widths vary, then differences are seen in the modified coefficients even for weak shocks, as shown in figure 14(b). This shows the two modified reflection coefficients for a range of width ratios. For a water–air interface, the magnitude of $\hat {C}_{\mathcal {R},+ve}$ is always less than would be predicted by linear theory. Whilst this coefficient is less dependent on the negative part of the pulse than the negative coefficient is on the positive part, there is still a weak dependency, because the rarefaction wave at the front of the reflected pulse has interacted with the negative portion of the incident wave. The reflection coefficient of the negative portion of the pulse is more sensitive to changes in the pulse parameters, which is to be expected because it depends more strongly on the positive portion as well as the negative portion. This again highlights the nonlinearity of the process even for weak shocks.

Figure 14. Modified reflection coefficients for: (a) both P-type and PN-type pulses (for the PN-type pulse, the widths and magnitudes of over- and under-pressures are equal); (b) the PN-type pulse with varying width and pressure ratios, where $p'_{\mathcal {S},+ve}=10$ MPa.

4.5. Strong shocks

As we move away from the weak shock regime, the assumption made in this study regarding rarefaction waves starts to break down. Figure 15 shows the pressure, velocity and density at time $t=2\times 10^{-7}\,\text {s}$ for a P-type pulse with $p_\mathcal {S}=5\times 10^8\,\text {Pa}$. The agreement between the analytical approach and the numerical simulations for the minimum/maximum pressure in the reflected/transmitted pulse is again very good, but there are discrepancies with the overall shape. This is because the assumption that the rarefaction waves remain infinitely thin over the relevant time scales becomes increasingly less valid as the shock pressure increases. The particle velocity in the shock and reflected waves is $284\,\text {m}\,\text {s}^{-1}$, and this leads to the rarefaction waves spreading out, even over short time scales. The analytical solution allows for the leftmost rarefaction wave to spread out for $t> t_2$, so it is not infinitely thin; but because it is not allowed to spread out for $t_1\leq t \leq t_2$, it appears thinner than is predicted by the simulation. The same effect is seen in the left-going shock that propagates into the rarefied liquid. In the analytical approach, this is a single shock resulting from the infinitely thin rarefaction arriving at the contact at $t=t_3$. In the numerical simulation, this wave has spread out slightly, effectively leading to the formation of a series of shocks that look like an inverse rarefaction wave. It should be noted that this problem could be ameliorated by splitting up the rarefaction waves into segments and then treating each segment as an individual wave. This would increase the number of wave interactions taking place, and hence the number of Riemann problems to be solved. Such an approach is permissible in this framework, but is beyond the scope of what we consider here.

Figure 15. Pressure, density, and velocity at $t=2\times 10^7\,\text {s}$ following interaction of a P-type pulse with $p_\mathcal {S}=5 \times 10^8\,\text {Pa}$ with a water–air interface.

Because the magnitude of the reflection coefficient reduces with increasing shock pressure (see figure 8), positive shock pressures greater than $p_\infty$ are permissible. From a mathematical standpoint, a solution can be obtained provided that the pressure does not drop below $-p_\infty$. When using the approach set out in this study, we suggest $p_\infty$ as an appropriate limit for the maximum shock pressure.

For the PN-type pulse, another limiting factor must be considered. At time $t=t_4$, the negative part of the incident pulse starts to transmit into the gas phase. The solution to the Riemann problem at this time determines the pressure and density behind the rarefaction wave that propagates into the gas phase. A pressure or density of zero would imply a vacuum, and this presents a limit to the model presented in this work. Figure 16 shows the pressure in the negative portion of the transmitted wave following the solution to the fourth Riemann problem as a function of the pressure in the negative portion of the incident pulse. As can be seen, for a water–gas interface at atmospheric pressure, the gas will be put into a vacuum state once the negative pressure in the incident shock drops below approximately $-$500 MPa. However, it can also be seen that this is dependent on the pressure in the positive part of the incident pulse as well as the negative part. Higher positive shock pressures require lower negative pressures to put the gas into vacuum. For $p'_{\mathcal {S},-ve}/p'_{\mathcal {S},+ve}=-0.5$, the model predicts that the pressure in the liquid will drop below $p_\infty$ before the pressure in the gas reaches zero. This further highlights the nonlinear behaviour that arises from the different wave interactions, particularly as the pressure magnitudes increase.

Figure 16. Pressure in the negative portion of the transmitted pulse in the gas phase.

It is important to note that the reflection of many of these stronger shocks would likely induce cavitation, which would radically alter the dynamics of the fluid, and it is beyond the scope of this work to consider this problem.

5. Reflection of non-square waves

To demonstrate the applicability of the approach set out in this work, we consider the reflection of two non-square, i.e. non-idealised, pulses. A modified Friedlander wave and a lithotripter wave are used for this as they represent practical realisations of the P-type and PN-type pulses.

The modified Friedlander equation describes a waveform that closely matches experimental measurements of blast waves (Dewey 2010). Based on the original Friedlander model (Friedlander 1946), the modified equation was introduced to provide better agreement for blast waves with larger peak over-pressures. It has the form

(5.1)\begin{equation} p'(t) = p'_\mathcal{S}\left(1-\frac{t}{t_I}\right)\exp({-bt/t_I}), \end{equation}

where $t_I$ is the over-pressure duration, and $b$ is a constant. Besides blast waves, the Friedlander model was also found to provide a good description for laser-induced shock waves in water by Bokman et al. (2023), who used it to investigate shock-induced bubble collapse. The model parameters considered in this work are based on the experiments of Bokman et al. (2023), with $t_I = 8.265 \times 10^{-7}$ s, $b=10$, and $p'_\mathcal {S}=3$, $30$ and $300$ MPa.

In the context of medical applications, shock waves generated by a lithotripter are often modelled (Church 1989; Johnsen & Colonius 2008) as

(5.2)\begin{equation} p'(t) = 2 p'_\mathcal{S} \,{\rm e}^{{-}a t} \cos{(2 {\rm \pi}f t + {\rm \pi}/3)}. \end{equation}

This pulse consists of a shock front followed by an exponentially decaying pressure, and includes a tensile part before returning to atmospheric pressure. The parameters considered in this work are $a=9.1 \times 10^5$ s$^{-1}$, $f = 83.3$ kHz (Shams, Bidi & Gavaises 2024), and $p'_\mathcal {S} = 1$, 10 and 100 MPa.

Figure 17(a) shows the non-dimensional pressure of the incident and reflected Friedlander waves for each over-pressure level. For the lowest over-pressure, the reflected waveform is similar to that of the incident wave. At higher over-pressures, the front of the reflected pulse spreads out and the pressure magnitude reduces, as also predicted by the analytical model. The reflection coefficients from the numerical simulations are $C_\mathcal {R} = -0.997, -0.970, -0.778$ compared with $C_\mathcal {R} = -0.997, -0.973, -0.780$ predicted by the analytical approach. This agreement suggests something important about the nature of the reflection process. From an analytical perspective, the reflection of this type of wave actually consists of an infinite number of Riemann problems. However, it is possible to predict the peak pressure of the reflected wave with good accuracy by considering only three.

Figure 17. Reflection of (a) a modified Friedlander pulse and (b) a lithotripter pulse, at a water–air interface for three different peak over-pressures.

Similar observations can be made for the reflection of the positive part of the lithotripter wave (figure 17b). Here, the coefficients predicted by the simulations are $C_{\mathcal {R},+ve} = -0.999, -0.986, -0.910$ compared with analytical predictions $C_{\mathcal {R},+ve} = -0.999, -0.990, -0.914$. Again, there is good agreement between the two approaches, and a reducing magnitude of the coefficient is predicted with increasing shock pressure. The reflection coefficients for the negative part are $C_{\mathcal {R},-ve} = -0.999, -0.992, -0.992$ for the numerical simulation, and $C_{\mathcal {R},-ve} = -0.999, -1.000, -1.001$ for the analytical model. Whilst the agreement is still fairly good, with the magnitude being higher than that for the positive part, there is a subtle difference. Because the negative part is not bounded by a discontinuity, the reflection behaves more like a linear process, and we do not see a significant change in the coefficient as a function of pressure.

6. Conclusions

This work has set out an approach for computing the reflection and transmission coefficients for an impulsive shock wave interacting with a liquid–gas interface. The problem is described as a series of Riemann problems for an idealised pulse, which allows for the pressure, density and velocity in the reflected and transmitted pulses to be computed analytically. The idealised pulses considered are a positive square wave (P-type pulse) and a square wave consisting of a positive portion ahead of a negative portion (PN-type pulse). Comparisons with numerical simulations show that the analytical approach accurately predicts the solution structure and also the pressure, density and velocity values for a wide range of incident shock pressures. Good agreement with numerical simulations was also observed for the reflection of non-idealised waves.

The proposed approach was used to explore the sensitivity of the reflection and transmission coefficients to changes in shock pressure, pulse width and, for the PN-type pulse, the ratios of widths and pressure magnitudes. In the acoustic limit, the method yields the same results as linear theory. Away from the acoustic limit, it was shown that defining a reflection coefficient using only the solution to a Riemann problem at the interface will lead to an over-prediction of the negative pressure. At least one additional Riemann problem must be solved to compute the pressure in the reflected pulse to account for the additional interactions between the reflected wave and the tail of the pulse.

For both square and non-square waves, the shape of the reflected wave changes as the shock pressure increases. The most notable change is the shape of the front: the shock front is a single discontinuity, whereas the front of the reflected wave is a rarefaction that spreads out as it propagates. This becomes more pronounced as the shock pressure increases due to the larger particle velocities. Furthermore, the magnitude of the reflection coefficient decreases as the incident shock pressure increases, and it is shown that this effect can be captured only by treating the interaction of the initially reflected rarefaction with the tail of the shock as a Riemann problem. Applying linear superposition of the reflected rarefaction wave and the ambient pressure will lead to erroneous predictions for all but the weakest shocks. For a PN-type pulse, it has also been shown how a region of negative pressure with a magnitude greater than that of the incident shock can exist, but only for a finite period of time, depending on the width of the negative part of the pulse. Accounting for this is important for applications where knowing the exact magnitude of the negative pressure is essential.

Future work should consider in more detail the case of strong shock interactions. In particular, an understanding of the dynamics where vacuum states are likely is needed. Non-planar waves should also be considered, as should a non-planar interface such as a bubble. Additionally, future research may involve analysing the performance of numerical methods for compressible multiphase flows, with a specific focus on the interaction between short shocks and interfaces.