2.1. Physical model

The previous study (Igra & Sun Reference Igra and Sun2010) has shown that the flow characteristics inside the two-dimensional water column are similar to those inside the spherical droplet, and that three-dimensional numerical simulation comes at a huge CPU time cost. Hence, the two-dimensional water column is chosen in the present study to save computing resources as much as possible and improve computing efficiency. A schematic diagram of the interaction of a cylindrical incident shock wave with a water column is shown in figure 1, including a converged and diverged one. We use $R_0$ to represent the radius of the cylindrical shock wave when it just touches the water column and $R_D$ to represent the radius of the water column. Moreover, the dimensionless radius ${\omega (=R_0/R_D)}$ is used to normalise the curvature effect of the incident shock wave. Referring to the experiment (Igra & Takayama Reference Igra and Takayama2001a) and numerical simulation (Xiang & Wang Reference Xiang and Wang2017), the initial value of $R_D$ is taken as 2.4 mm. Moreover, we use $C$ and $O$ to represent the water column centre and the origin of coordinates, respectively. The water column and the air ahead of the shock wave are initially in equilibrium with a temperature of 300 K and a pressure of 101 325 Pa. The Weber numbers in the present numerical simulations are higher than 1000, and the corresponding Reynolds numbers are over 40 000. Therefore, the viscosity and the surface tension can be neglected in the present study (Meng & Colonius Reference Meng and Colonius2015).

Figure 1.

Schematic diagram of the interaction of a cylindrical shock with a water column. (a) The interaction of a converged shock with a water column; (b) the interaction of a diverged shock with a water column.

In the present study, the generation of the cylindrical shock wave is based on the theory of shock dynamics, which can characterise the propagation of shock waves with an arbitrary profile. The Chester–Chisnell–Whitham (CCW) relation is the basis of shock dynamics for a uniform quiescent gas ahead of shock, which is referred from Chester (Reference Chester1954), Chisnell (Reference Chisnell1957) and Whitham (Reference Whitham1957, Reference Whitham1958, Reference Whitham1959). The CCW relation describes how the shock wave Mach number $M$ varies with the shock-front area $A$ ($= 2{\rm \pi} R$ in two-dimensional cases), which can be written as

(2.1)\begin{equation} \frac{2M\,\mathrm{d}M}{({{M}^{2}}-1)K( M )}+\frac{\mathrm{d}A}{A}=0, \end{equation}

where

(2.2)\begin{align} K( M )=2{{\left[ 2\mu +1+\frac{1}{{{M}^{2}}} \right]}^{{-}1}}{{\left[ 1+\frac{2}{\gamma +1}\left( \frac{1}{\mu }-\mu \right) \right]}^{{-}1}} ,\quad \mu =\sqrt{\frac{( \gamma -1 ){{M}^{2}}+2}{2\gamma {{M}^{2}}-( \gamma -1 )}}. \end{align}

Zhai et al. (Reference Zhai, Liu, Qin, Yang and Luo2010) used this shock dynamics theory to design a curved wall profile (V-shaped geometry) and realised the transformation from a planar shock wave to a cylindrical converged shock wave, as shown in figure 2. Both numerical and experimental results show a perfect circular shock front. However, this method comes at enormous computing resource cost due to the large computation domain and usually can only obtain a cylindrical shock wave with a small converging angle ($= 2\theta _0$), which leads to a small range of radius of the water column that can be investigated in the present study. Hence, a $90^{\circ }$ computational domain is used in this paper, as shown in figure 1($a$), which is widely used in the literature (Lombardini, Pullin & Meiron Reference Lombardini, Pullin and Meiron2014; Zhai et al. Reference Zhai, Zhang, Zhou, Ding and Wen2019; Wu, Liu & Xiao Reference Wu, Liu and Xiao2021a). The fluid variables behind the cylindrical shock wave are calculated by the following algorithm: given the cylindrical shock radius of $R_0$, the cylindrical shock wave intensity $M_0$ and fluid variables ahead of the incident cylindrical shock wave $p_0$ and $T_0$. When the radius of the cylindrical shock wave is $R$, the intensity of the cylindrical shock wave $M$ satisfies

(2.3) \begin{equation} \int_{{{M}_{0}}}^{M}{\frac{2m\,\mathrm{d}m}{({{m}^{2}}-1)K( m )}}={-}\int_{{{R}_{0}}}^{R}{\frac{\mathrm{d}r}{r}}. \end{equation}

Figure 2.

The schematic drawing of the wall profile, transforming the planar shock into a cylindrical one. Here, $H$ is the half-height of the V-shaped geometry; $l$ the length of the V-shaped geometry; $Ma_0$ the Mach number of the incident planar shock waves; $Ma_D$ the Mach number of the cylindrical shock wave; $\theta _0$ the half-converging angle.

Integrating (2.3) by numerical iteration, the value of $M$ is obtained. Further, using the Rankine–Hugoniot conditions, the fluid variables behind the cylindrical shock wave can be calculated, including pressure $p$, density $\rho$, temperature $T$ and velocity magnitude $V_m$. The accuracy verification of CCW theory is presented in detail in Appendix A.

2.2. Method of the ray analysis

In this sub-section, the ray analysis method (Heymann Reference Heymann1969; Haller et al. Reference Haller, Poulikakos, Ventikos and Monkewitz2003; Wu et al. Reference Wu, Xiang and Wang2018), which is based on the Huygens principle used in the theoretical analysis of wave configurations, is presented. For convenience, two moving point disturbances, which are denoted as $S_1$ and $S_2$, and a straight material interface are used to briefly introduce the ray analysis method, as shown in figure 3. At the initial instant, two moving point disturbances coincide at point $S_0$, located on the straight interface. Then these two point disturbances move away from $S_0$ in opposite directions with the same speed, which gradually decreases with distance s from $S_1$ ($S_2$) to $S_0$ and is denoted as $V(s)$. Based on the Huygens principle, an individual wavelet will be emitted at each instant at $S_1(s)$ and these individual wavelets propagate at a constant speed of $W$. It can be seen that the envelope of all individual wavelets generated before $t$ forms the wavefront induced by these two moving point disturbances. According to the relationship between $V(s)$ and $W$, there are four envelope shapes, as shown in figure 3(a–d). Note that, for the case corresponding to figure 3(b), when the propagation speed of the point disturbance $V(s)$ is less than $W$, the newer individual wavelets cannot catch up with the envelope of wavelets, as will be explained in detail in § 3.1.

Figure 3.

Schematic diagram of the Huygens principle.

For three cases except for the case corresponding to figure 3(b), it is not difficult to understand the position and shape of the envelope of wavelets, while it is still confusing to understand the evolution characteristics of the envelope for the rest of the cases. Hence, the ray analysis method is used, as shown in figure 4. The propagation of each wavelet can be equivalent to the propagation of infinite rays emitted at its origin, and the length of the ray is equal to the propagation distance of the wavelet. It is obviously found that not all of the rays emitted from the same wavelet can effectively contribute to the envelope of wavelets. Therefore, to reveal the physical mechanism of the motion of the wavefront induced by point disturbances, it is necessary to find these special rays which have contributed to the envelope of wavelets. Here, the two wavelets emitted from $S_s$ and $S_{s+\Delta s}$ are selected for detailed discussion, as shown in figure 4(b), and these two emission points are very close to each other. The radii of these two wavelets are $r(t,s)$ and $r(t, s+\Delta s)$, respectively,

(2.4)\begin{equation} r( {t,s} ) = W( {t - {t_s}} ) ,\quad {r( {t,s + \Delta s} ) = W( {t - {t_{s + \Delta s}}} )}, \end{equation}

where $t_s$ represents the time for moving point disturbance $S_1$ from $S_0$ to $S_1(s)$.

Figure 4.

Schematic diagram of the ray analysis method.

The intersection point of these two wavelets is $G_s$, and the angle between the vector $\overrightarrow {{S_s}{S_{s+\Delta s}}}$ and $\overrightarrow {{S_s}{G_{s}}}$ is $\alpha _s$, as shown in figure 4(b). The expression $\alpha _s$ is derived from the law of cosines in the triangle $\Delta {S_s}{S_{s+\Delta s}}{G_s}$

(2.5) \begin{equation} \cos {\alpha_s} = \frac{{( {2t - {t_s} - {t_{s + \Delta s}}} )}}{{2( {t - {t_s}} )}} \frac{{( {{t_{s + \Delta s}} - {t_s}} )W}}{{\Delta s}} + \frac{{\Delta s}}{{2( {t - {t_s}} )W}}. \end{equation}

When $S_{s+\Delta s}$ is infinitely close to $S_s$ ($\Delta s\to {0^+}$), $\alpha _s$ is the angle between the vector of the emitted ray and the material interface. The endpoint of this ray is the unique contribution of the wavelet emitted from $S_s$ to the envelope of wavelets

(2.6)\begin{align} \cos {\alpha_s}&=\lim_{\Delta s \to {0^ + }} \left[ {\frac{{( {2t - {t_s} - {t_{s + \Delta s}}} )}}{{2( {t - {t_s}} )}}\frac{{( {{t_{s + \Delta s}} - {t_s}} )W}}{{\Delta s}} + \frac{{\Delta s}}{{2( {t - {t_s}} )W}}} \right] \nonumber\\ &= \lim_{\Delta s \to {0^ + }} \frac{{( {{t_{s + \Delta s}} - {t_s}} )W}}{{\Delta s}} = W{t'_s}. \end{align}

Equation (2.6) shows that the emission angle of these particular rays is only determined by the propagation speed of the wavelets and the kinematic characteristics of the moving point disturbance along the material interface. Note that, if the material interface were a closed curve, these rays would reflect on the material interface after a period of propagation. This reflection property will be investigated in detail in § 3.1.

2.3. Numerical models

The interaction between the shock wave and the droplet is a strong compressible multiphase hydrodynamics problem, which involves complex factors such as a large density ratio and strong shock waves. In this paper, the numerical simulation is carried out by the in-house software (SCP-tran), which was previously applied to study a variety of compressible multiphase flow problems (Xiang & Wang Reference Xiang and Wang2017; Wang, Xiang & Hu Reference Wang, Xiang and Hu2018; Wu et al. Reference Wu, Xiang and Wang2018). The five-equation model (Allaire, Clerc & Samuel Reference Allaire, Clerc and Samuel2002; Johnsen & Colonius Reference Johnsen and Colonius2006) is used to solve the gas–liquid hydrodynamic system, and the governing equations consist of two continuity equations for each phase, a mixture momentum equation, a mixture energy equation and a volume fraction advection equation of the liquid phase

(2.7)\begin{equation} \left.\begin{gathered} \frac{{\partial ( {{\alpha_l}{\rho_l}} )}}{{\partial t}} + \frac{{\partial ( {{\alpha_l}{\rho_l}u} )}}{{\partial x}} + \frac{{\partial ( {{\alpha _l}{\rho_l}v} )}}{{\partial y}} = 0,\\ \frac{{\partial ( {{\alpha_g}{\rho_g}} )}}{{\partial t}} + \frac{{\partial ( {{\alpha_g}{\rho_g}u} )}}{{\partial x}} + \frac{{\partial ( {{\alpha _g}{\rho_g}v})}}{{\partial y}} = 0,\\ \frac{{\partial ( {\rho u} )}}{{\partial t}} + \frac{{\partial ( {\rho {u^2} + p} )}}{{\partial x}} + \frac{{\partial ( {\rho uv} )}}{{\partial y}} = 0,\\ \frac{{\partial ( {\rho v} )}}{{\partial t}} + \frac{{\partial ( {\rho uv} )}}{{\partial x}} + \frac{{\partial ( {\rho {v^2} + p})}}{{\partial y}} = 0,\\ \frac{{\partial \rho E}}{{\partial t}} + \frac{{\partial [ {( {\rho E + p} )u} ]}}{{\partial x}} + \frac{{\partial [ {( {\rho E + p})v} ]}}{{\partial y}} = 0,\\ \frac{{\partial {\alpha_l}}}{{\partial t}} + u\frac{{\partial {\alpha_l}}}{{\partial x}} + v\frac{{\partial {\alpha_l}}}{{\partial y}} = 0, \end{gathered}\right\} \end{equation}

where $\rho _l$ and $\rho _g$ represent the density of the liquid and gas phases, respectively, $\alpha _l$ and $\alpha _g$ represent the volume fraction of the liquid and gas phases, respectively, $\rho$, $u$, $v$, $p$ and $E$ represent the mixture density, $x$-velocity, $y$-velocity, pressure and specific total energy, respectively. The numerical diffusion appears significant at the two-phase interface after several time steps. In this diffuse region, the mixture variables are given as (Saurel, Petitpas & Abgrall Reference Saurel, Petitpas and Abgrall2008)

(2.8)$$\begin{gather} \rho = {\alpha_l}{\rho_l} + {\alpha_g}{\rho_g}, \end{gather}$$

(2.9)$$\begin{gather}\rho E = {\alpha_l}{\rho_l}{e_l} + {\alpha_g}{\rho _g}{e_g} + \tfrac{1}{2}\rho ( {{u^2} + {v^2}} ), \end{gather}$$

(2.10)$$\begin{gather}\rho {c^2} = {\alpha_l}{\rho_l}c_l^2 + {\alpha_g}{\rho _g}c_g^2, \end{gather}$$

where $e_l$ and $e_g$ represent the specific internal energy of the liquid and gas phases, respectively, $c_l$ and $c_g$ represent the sound speed of the liquid and gas phases, respectively. In this study, the stiffened gas equation of state is used to close the governing equations

(2.11)$$\begin{gather} {{\rho_k}{e_k} = \frac{{p + {\gamma_k}{p_{\infty ,k}}}}{{{\gamma_k} - 1}}} ,\quad {k = l,g}, \end{gather}$$

(2.12)$$\begin{gather}{{c_k} = \sqrt {\frac{{\left( {p + {p_{\infty ,k}}} \right){\gamma_k}}}{{{\rho_k}}}} },\quad {k = l,g}, \end{gather}$$

where $\gamma$ is the specific heat ratio and $p_\infty$ is the reference pressure. For air, $\gamma = 1.4$ and $p_\infty = 0$, and the stiffened gas equation of state reduces to the ideal gas equation. Referring to the works of Kondo & Ando (Reference Kondo and Ando2016) and Xiang & Wang (Reference Xiang and Wang2017), the parameters for water are taken to be $\gamma = 6.12$ and $p_\infty = 343.44$ MPa.

2.4. Numerical treatments

The SCP-tran fluid dynamic software uses a finite volume method (Titarev & Toro Reference Titarev and Toro2004) to discretise the above governing equations in a uniform Cartesian grid system. The component-wise fifth-order incremental weighted essentially non-oscillatory reconstruction is applied, as previously proposed by the present author (Wang et al. Reference Wang, Xiang and Hu2018). The Harten–Lax–van Leer contact approximate Riemann solver (Toro Reference Toro2009) is employed to solve the numerical flux at the cell face. The third-order total variation diminishing Runge–Kutta method (Gottlieb & Shu Reference Gottlieb and Shu1998) is chosen to advance the solutions over time. Considering that part of the boundaries does not coincide with the interface of the grid cell, such as the immersed non-reflecting boundary as shown in figure 1, a ghost-cell immersed boundary method (IBM) for distinguishing geometrically complex boundaries is used to realise the non-reflecting boundary condition (Thompson Reference Thompson1987, Reference Thompson1990), and for the detail of IBM the reader is referred to the works of Choung et al. (Reference Choung, Saravanan, Lee and Cho2021) and Saravanan, Choung & Lee (Reference Saravanan, Choung and Lee2021). Since the evolution characteristics of the interaction between shock waves and the water column are symmetric, to improve the calculation efficiency, the symmetric boundary condition at the axis of the liquid column is adapted to carry out numerical simulations. The numerical verification of the grid sensitivity is presented in Appendix B.

2.5. Validation of the solver

Thanks to the experiments in the literature, the interaction of a planar shock wave with a water column is chosen as a validation case to validate the SCP-tran fluid dynamic software. For a qualitative comparison, the dimensionless time $t^*$ is used, which is the ratio of the physical time over the characteristic time $\tau$ ($\tau = 2R_D/ V_{ts}$, $V_{ts}$ is the propagation speed of the transmitted shock wave). For convenience, the zero instant is marked as $t_0$ ($t^* = 0.0$) when the incident shock wave just touches the water column. According to the Rankine–Hugoniot relation (Haller et al. Reference Haller, Ventikos, Poulikakos and Monkewitz2002, Reference Haller, Poulikakos, Ventikos and Monkewitz2003; Nagayama et al. Reference Nagayama, Mori, Motegi and Nakahara2006), the transmitted shock wave speed $V_{ts}$ can be estimated as

(2.13)\begin{equation} {V_{ts}} = \frac{{{\gamma_l} + 1}}{4}\left( {{u_l} + \sqrt {u_l^2 + 16\frac{1}{{{{( {{\gamma_l} + 1} )}^2}}}c_{l,0}^2} } \right), \end{equation}

where $\gamma _l$ represents the specific heat ratio of water, $c_{l,0}$ ($\sim$1500.0 m s$^{-1}$) represents the sound speed of the water at the initial state (300 K and 101 325 Pa) and $u_l$ represents the velocity of the liquid inside the water column behind the transmitted shock wave. Due to the acoustic impedance of water being much higher than that of air, most of the energy of the incident shock wave is reflected. Hence, the velocity change of the water is almost zero, no bigger than 5.0 m s$^{-1}$, before and after the impingement of the incident shock wave, and the sound speed of water $c_{l,0}$ can be chosen as an equivalent of $V_{ts}$.

Figure 5 shows the comparisons between the numerical result based on the present numerical methods and the experimental results from Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016) of the interaction of a planar shock wave with a water column for $M_0 = 2.4$. When the planar incident shock wave impinges on the water column, it is reflected off the column surface and transmitted into the water column. As the interaction continues, the Mach stem and the slip line appear subsequently. Due to the specific water column surface, the reflected rarefaction wave (Sembian et al. Reference Sembian, Liverts, Tillmark and Apazidis2016) (also called the reflected expansion wave) will focus inside the water column. These flow structures have been observed both in experiments and simulations. Besides the qualitative analysis, the quantitative comparison of the pressure profiles of the two sensors inside the water column between the experimental result and numerical simulations is also presented in figure 6. It is found that our numerical method (SCP-tran) can effectively capture the pressure evolutions inside the water column, and is in good agreement with the numerical result and is also approximately in agreement with the experimental result. Due to the SCP-tran software having higher accuracy of the time advance and the space discretisation causing lower numerical dissipation, compared with the work of Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016), the present numerical result obtains stronger focus pressures, especially for the focus of the second reflected waves, as shown in figure 6(a). Note that the sensing geometry is compressed; hence, the sensor cannot quantitatively measure the negative pressure in experimental research. Both the qualitative wave configuration analysis and the quantitative pressure profile analysis show that the present mathematical models and numerical methods are able to solve the problem of the interaction of the shock wave with a water column.

Figure 5.

The comparison between the present simulation results (left side) and the experimental results (right side) from Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016) of the interaction of a planar shock wave with a water column for $M_0 = 2.4$.

Figure 6.

Experimental and numerical pressure profile for $M_0 = 2.4$. The locations of sensors 2 and 3 are given in the work of Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016). Note that the diaphragm diameter of the sensor is 5.54 mm, and therefore the values obtained are averaged across the sensor's face area. Similarly, numerical simulation results are also averaged.

It is necessary to illustrate that the planar incident wave, in experiments (Sembian et al. Reference Sembian, Liverts, Tillmark and Apazidis2016), is a blast wave with gradually decreasing strength rather than a shock wave with a uniform post-shock flow. This substitution causes different pressure distributions inside the water column and then induces different cavitation behaviours. However, the main purpose of the present work is to demonstrate the propagation characteristics of waves inside the water column, which are almost the same between the blast wave with a decreasing post-shock flow and the shock wave with a uniform post-shock flow when the shock intensity is the same. Hence, the planar shock with a constant Mach number is chosen in the present study for simplification as much as possible and to improve computing efficiency. Of course, when it is necessary to deeply understand the evolution characteristics of the cavitation phenomenon inside a water column impacted by a weak shock wave, an analytical solution for blast waves, described by Bach & Lee (Reference Bach and Lee1970), can be chosen to generate a stable blast wave without much increase in computational cost.

3.1. The generation and propagation of the transmitted shock wave

The first stage begins at $t_0$ ($t^* = 0.0$) when the incident curved shock wave just impacts the left pole of the water column, figure 7(a). During the interaction of a cylindrical converged shock wave with a water column, the contact point is denoted as $P_\theta$, and the angle $\theta$ represents the angle between the line $\overline {{P_\theta }C}$ and the horizontal axis of the water column, as shown in figure 8. The radius of the cylindrical shock wave is a function of angle $\theta$ and is denoted as $R_\theta$. The angle between the tangent line of the water column and that of the converged shock wave is denoted as $\chi$

(3.1)$$\begin{gather} {R_\theta } = {R_D}\sqrt {1 + {{( {\omega - 1} )}^2} + 2( {\omega - 1} )\cos \theta }, \end{gather}$$

(3.2)$$\begin{gather}\sin \chi = \frac{{( {\omega - 1} )\sin \theta }}{{\sqrt {1 + {{( {\omega - 1} )}^2} + 2( {\omega - 1} )\cos \theta } }}. \end{gather}$$

Figure 8.

Schematic diagram of the interaction between the cylindrical converged shock wave and the water column.

At a specific contact angle $\theta$, the velocity of contact point $P_\theta$ along the water column surface can be expressed by

(3.3)\begin{equation} {V_P} = \frac{{{V_S}}}{{\sin \chi }} = {V_S}\frac{{\sqrt {1 + {{( {\omega - 1} )}^2} + 2( {\omega - 1} )\cos \theta } }}{{( {\omega - 1} )\sin \theta }}, \end{equation}

where the velocity of the cylindrical converged shock wave $V_S$ is the product of the sound speed of air ahead of the shock front $c_{g,0}$ and the Mach number of the cylindrical shock wave $M_\theta$ at $P_\theta$. Because the strength $M_0$ and the radius $R_0$ of a cylindrical converged shock wave at $t_0$ are known, the strength of the cylindrical shock wave at $P_\theta$ can be solved from (2.3).

When the contact angle is close to zero, the velocity $V_P$ is much higher than the sound speed of the water, and the transmitted shock wave is attached to the cylindrical incident shock wave and reflected shock wave at the water column surface. Since $V_P$ decreases rapidly as $\theta$ increases, as shown in figure 9, the velocity of the contact point on the column surface will catch up with the propagation speed of the transmitted shock wave, and the confined transmitted shock wave will detach from the incident shock wave. The critical time, representing the transmitted shock wave just being detached, is defined as $t_{cr}$. At the critical time, the velocity of contact point $V_P$ equals the propagation velocity of transmitted wave $V_{ts}$, and the critical contact angle $\theta _{cr}$ satisfies the following expression:

(3.4)\begin{equation} \frac{{\sqrt {1 + {{( {\omega - 1} )}^2} + 2( {\omega - 1} )\cos {\theta_{cr}}} }}{{( {\omega - 1} )\sin {\theta_{cr}}}} = \frac{{{c_{l,0}}}}{{{M_{\theta cr}}{c_{g,0}}}}, \end{equation}

where $M_{\theta cr}$ represents the strength of the incident cylindrical shock wave at $t_{cr}$. In this section ($\omega = 4.0$ and $M_0 = 2.4$), the value of the critical angle is 44.9$^\circ$.

Figure 9.

The velocity of the contact point $P$, along the column surface, varies with the contact angle $\theta$.

Once $\theta$ is larger than $\theta _{cr}$, the precursor transmitted shock wave is formed and propagates to the right pole of the water column. The time for contact point to move from the left pole of the water column to $P_\theta$ is denoted as $t_\theta$, and its expression can be written as

(3.5) \begin{equation} {t_\theta } = \int_0^\theta {\frac{{( {\omega - 1} )\sin \xi {R_D}}}{{\sqrt {1 + {{( {\omega - 1} )}^2} + 2( {\omega - 1} )\cos \xi } {M_\xi }{c_{g,0}}}}\,\text{d}\xi }. \end{equation}

Combined with (3.4), the critical instant $t_{cr}$ can be obtained when the transmitted shock wave detaches from the incident shock wave.

Based on the ray analysis method introduced in § 2.2, at each instant, an individual compression wavelet will be emitted at $P_\theta$ that will propagate with the water sound speed $c_{l,0}$. Hence, the radius of the compression wavelet emitted from the contact point $P_\theta$, at time instant $t$, can be denoted as $r(\theta,t) = (t - t_\theta )c_{l,0}$. These compression wavelets, emitted from different contact points from zero instant to time instant $t$, form a shock wave envelope, denoted as $TS_t$, figure 10(a). When $t > t_{cr}$, the propagation speed of the compression wavelet is higher than the generation speed of the new compression wavelet, and the new compression wavelet cannot catch up with the transmitted shock wave and does not make an effective contribution to the shock wave envelope, forming the precursor transmitted shock wave, as shown in figure 10(b).

Figure 10.

Schematic diagram of the generation of transmitted shock wave and the ray analysis: ($a$) the schematic diagram at critical time $t_{cr}$; ($b$) the schematic diagram at the time instant $t_1$, selected after the critical time.

Similar to the discussion in § 2.2, it can be obviously found that not all of the rays emitted from the same compression wavelet make an effective contribution to the shock wave envelope of compression wavelets. Therefore, to reveal the physical mechanism of the motion of the transmitted shock wave, it is necessary to find these special rays which have contributed to the envelope of compression wavelets. Here, the two different compression wavelets emitted from $P_\theta$ and ${{P_{\theta + \Delta \theta }}}$ are selected for detailed discussion, and these two emission points are infinitely close to each other (${\Delta \theta \to {0^ + }}$). The radii of these two compression wavelets are $r(\theta,t)$ and $r(\theta + \Delta \theta,t)$, respectively. Similar to (2.5) and (2.6), the emission angle $\alpha _\theta$, which is the angle between the vector of the emitted ray and the tangent vector of the water column at $P_\theta$, can be derived as

(3.6) \begin{align} \cos {\alpha_\theta } &= \lim_{\Delta \theta \to {0^ + }} \left[ {\frac{{( {2t - {t_\theta } - {t_{\theta + \Delta \theta }}} )}}{{2( {t - {t_\theta }} )}} \frac{{( {{t_{\theta + \Delta \theta }} - {t_\theta }} ){c_{l,0}}}}{{{R_D}\Delta \theta }}\frac{\Delta \theta/2}{\sin(\Delta \theta/2)} + \frac{{{R_D}\sin ( {{{\Delta \theta } /2}} )}}{{( {t - {t_\theta }}){c_{l,0}}}}} \right]\nonumber\\ &= \lim_{\Delta \theta \to {0^ + }} \frac{{( {{t_{\theta + \Delta \theta }} - {t_\theta }} ){c_{l,0}}}}{{{R_D}\Delta \theta }} = \frac{{{c_{l,0}}}}{{{R_D}}}{{t'}_\theta } = \frac{\kappa }{{{M_\theta }}}\frac{{( {\omega - 1} )\sin \theta }}{{\sqrt {1 + {{( {\omega - 1} )}^2} + 2( {\omega - 1} )\cos \theta } }}, \end{align}

where the dimensionless parameter $\kappa$ represents the ratio of the sound speed of water to that of air. It is worth noticing that, for the interactions of a planar shock wave (the radius ratio $\omega$ tends to the infinite and the incident shock intensity is a constant $M_0$) with a water column, the angle $\alpha _\theta$ satisfies $\cos {\alpha _\theta } = \kappa \sin \theta /{M_0}$. This theoretical result is the same as the acoustic principle used by Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021) in studying the interaction between a planar shock wave and a water column. According to (3.4), it can be concluded that $\alpha _{\theta cr}$ of the ray generated by the critical contact point $P_{\theta cr}$ equals zero. That is to verify that when $t>t_{cr}$, the transmitted shock wave presents a precursory characteristic. Only one particular ray out of the infinite number of rays generated by the same contact point, whose emission angle satisfies (3.6), effectively influences the envelope of compression wavelets. Hence, only these specific rays emitted from different contact points ($\theta >\theta _{cr}$) are considered and analysed in the subsequent analysis, and the influence of the other rays is ignored.

As the transmitted shock wave propagates forwards, the reflected rarefaction waves are observed behind the transmitted shock and a certain time, $t$ ($t^* = 0.8695$) is chosen for the following analysis, figure 11($a,b$). The position and shape of the reflected rarefaction waves are obtained from the analysis of the emitted rays. Meanwhile, it is assumed that the rays will be reflected symmetrically on the curved column surface (Wu et al. Reference Wu, Xiang and Wang2018). At a specific instant $t$, the length of the ray emitted by the contact point $P_\theta$ is $r(\theta,t)$. Moreover, if the ray is reflected from the water column surface, $r$ will represent the total length of the ray before and after reflection. Furthermore, the ray could be reflected more than once from the column surface, and the reflection times are related to the contact angle $\theta$ and the time $t$, which are elaborated as follows.

Figure 11.

The schematic diagram of the ray analysis: ($a$) the schematic diagram at $t$ ($t^* = 0.8695$); ($b$) the enlarged view of the schematic diagram at $t$. ($c$) The comparison of results between ray analysis and numerical simulation at $t^* = 0.9543$.

If the rays emitted from $P_\theta$ are not reflected at an instant $t$, the emitting angle $\alpha _\theta$ will satisfy the following:

(3.7)\begin{equation} {\alpha_\theta } \ge \arcsin \frac{{r( {\theta ,t} )}}{{2{R_D}}} = \arcsin \left[ {\left( {t - \int_0^\theta {\frac{{{{L'}_\xi }}}{{{M_\xi }{c_{g,0}}}}\,{\rm d}\xi } } \right)\frac{{{c_{l,0}}}}{{2{R_D}}}} \right] = \alpha_\theta^{(1)}, \end{equation}

where ${L'_\theta }$ represents the derivative of the equivalent propagation distance $L_\theta$ of the curved shock along the symmetrical axis of the water column for the contact angle $\theta$

(3.8)\begin{equation} {L'_\theta } = \frac{{( {\omega - 1} )\sin \theta }}{{\sqrt {1 + {{( {\omega - 1} )}^2} + 2( {\omega - 1} )\cos \theta } }}{R_D}. \end{equation}

If the rays are reflected $N$ times ($N = 1, 2, 3, \ldots$), angle $\alpha _\theta$ will satisfy

(3.9)$$\begin{gather} {\alpha_\theta } < \arcsin \frac{{r( {\theta ,t} )}}{{2N{R_D}}} = \arcsin \left[ {\left( {t - \int_0^\theta {\frac{{{{L'}_\xi }}}{{{M_\xi }{c_{g,0}}}}\,{\rm d}\xi } } \right)\frac{{{c_{l,0}}}}{{2N{R_D}}}} \right] = \alpha_\theta^{(N)}, \end{gather}$$

(3.10)$$\begin{gather}{\alpha_\theta } \ge \arcsin \frac{{r( {\theta ,t} )}}{{2( {N + 1} ){R_D}}} = \arcsin \left[ {\left( {t - \int_0^\theta {\frac{{{{L'}_\xi }}}{{{M_\xi }{c_{g,0}}}}\,{\rm d}\xi } } \right)\frac{{{c_{l,0}}}}{{2\left( {N + 1} \right){R_D}}}} \right] = \alpha_\theta^{(N + 1)}. \end{gather}$$

Thus, the emission angle $\alpha _\theta$ can be divided into different intervals according to reflection times $N$, and the intervals of angle $\theta$ can be obtained by combining with (3.6).

Due to the specific column surface, the reflected rarefaction wave has two branches: the far branch (Re-RW$_{II}$) and the near branch (Re-RW$_{I}$), as judged by the distance of the branch from the axis of the water column. Similarly, the reflected rarefaction will also reflect on the water column surface and forms the second reflected wave, it has two branches with entirely different properties: the second reflected compression wave (Re$_2$-CW) and the second reflected rarefaction wave (Re$_2$-RW). It deviates from the statement (Xiang & Wang Reference Xiang and Wang2017; Wu et al. Reference Wu, Xiang and Wang2018) that the second reflection forms the compression wave with a single characteristic, but the phenomenon obtained by the present study is consistent with the experimental result of Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016). Nonetheless, we will not elaborate further on this in the present study.

The comparison of results between ray analysis and numerical simulation at $t^* = 0.9543$ is shown in figure 11($c$), and more details about the comparison between the numerical simulation results and theoretical results are presented in Appendix C. It is obvious that the shape and position of wave structures obtained by theoretical analysis almost coincide with the distribution of wave structures in the numerical results. Hence, this verifies the reliability and accuracy of the theoretical analysis in predicting the motion characteristics of the wave structures. Moreover, due to the strong stretched effect, a negative-pressure region is found near the intersection of two branches of the first reflected rarefaction wave. Meanwhile, the pressure behind the far branch of the first reflected rarefaction waves is quickly recovered because the subsequently emitted compression wavelets catch up with these two rarefaction wave branches.

3.2. The first convergence of the first reflected rarefaction wave

The first stage ends when the transmitted shock wave reaches the right pole of the water column at the time instant $t_2$ ($t^* = 1.0$), and the second stage begins. At this moment, the transmitted shock wave is completely reflected, and two near branches of the reflected rarefaction wave on both sides of the central axis of the water column merge, forming a continuous converged rarefaction wave (Re-RW$_{IC}$). The analytical schematic is presented in figure 12($a$), which demonstrates the first reflected rarefaction wave evolution. It is observed that the intersection points of Re-RW$_{II}$ and Re-RW$_{I}$ move along the envelope of the one-time reflected rays from time $t_2$ to $t_3$, and the continuous reflected rarefaction wave (Re-RW$_{IC}$) gradually focuses inside the water column. When Re-RW$_{IC}$ completely focuses, the two branches of the far-branch rarefaction wave meet and form a continuous diverged rarefaction (Re-RW$_{IIC}$).

Figure 12.

($a$) Schematic diagram of the first reflected expansion wave propagation from $t_2$ ($t^* = 1.0$) to $t_3$ ($t^* = 1.2698$); ($b$) schematic diagram of ray analysis for the focusing of the one-time reflected rays.

Previous studies (Sembian et al. Reference Sembian, Liverts, Tillmark and Apazidis2016; Xiang & Wang Reference Xiang and Wang2017; Wu et al. Reference Wu, Xiang and Wang2018; Biasiori-Poulanges & El-Rabii Reference Biasiori-Poulanges and El-Rabii2021) have found that the complete focus of the reflected rarefaction wave decreases the water pressure rapidly, and the maximum negative pressure inside the water column reaches approximately $-$10 MPa at the complete focus instant. The pressure distributions along the centre axis of the water column, just before and after the complete focus instant, are presented in figure 13, and it can be seen that a similar focusing pressure ($-$10 MPa) is obtained in the present study. This extremely negative pressure has far exceeded the cavitation threshold pressure of $-$2.3 MPa of unpurified water reported by Sembian et al. (Reference Sembian, Liverts, Tillmark and Apazidis2016). Therefore, the region inside the water near the focus point has a high probability of capturing the cavitation phenomenon. Hence, the focus point of the first reflected rarefaction wave is called the cavitation kernel point in the present study, denoted as $P_{cav}$. Combined with figure 12($b$), $P_{cav}$ is also the left limiting position of the one-time reflected rays that intersect with the central axis of the water column.

Figure 13.

The pressure distribution along the centre axis of the water column, just before and after the complete focus instant, for the interaction between the cylindrical converged shock and the water column.

The position of $P_{cav}$ can be obtained theoretically according to the ray analysis (Obreschkow et al. Reference Obreschkow, Dorsaz, Kobel, de Bosset, Tinguely, Field and Farhat2011; Wu et al. Reference Wu, Xiang and Wang2018; Biasiori-Poulanges & El-Rabii Reference Biasiori-Poulanges and El-Rabii2021) considering the planar incident shock wave or flat rigid surface. Here, we present the case analysis considering the curved, cylindrical, incident shock wave.

Some rays emitted by the contact point whose contact angle is less than the critical contact angle $\theta _{cr}$ are presented in figure 12($b$). In order to ensure that the one-time reflected rays can intersect with the horizontal central axis, the emission angle $\alpha _\theta$ and the contact angle $\theta$ will satisfy $4{\alpha _\theta } + \theta > 180^\circ$. Selecting an arbitrary ray (the black solid lines) for analysis, the intersection points of the one-time reflected part of this ray with the column surface and the horizontal central axis are $E_\theta$ and $F_\theta$, respectively. And the distance $\lambda _\theta$ from $F_\theta$ to the column centre can be expressed as follows:

(3.11) \begin{equation} {\lambda_\theta } = \frac{{\sin {\beta_\theta }}}{\sin {\psi_\theta }}{R_D} = \frac{{{R_D}}}{{\sin \theta \sin {\alpha_\theta }( {4\cos {\alpha_\theta } - {1 /\cos \alpha_\theta}} ) + \cos \theta ( {3 - 4\cos ^2{\alpha_\theta }} )}}. \end{equation}

After substituting (3.6) into (3.11) and simplifying, the expression of distance $\lambda _\theta$ can be rewritten as follows:

(3.12)\begin{equation} {\lambda_\theta } = \frac{{{R_D}}}{{{F_1}( \theta ) - {F_2}( \theta ) + {F_3}( \theta )}}, \end{equation}

where, ${F_1}( \theta )$, ${F_2}( \theta )$ and ${F_3}( \theta )$ are shown as

(3.13)$$\begin{gather} {F_1}( \theta ) = \frac{{4\kappa ( {\omega \,{-}\, 1} ){{\sin }^2}\theta }}{{\sqrt {1 + {{( {\omega \,{-}\, 1} )}^2} + 2( {\omega \,{-}\, 1} )\cos \theta } {M_\theta }}}\sqrt {1 - \frac{{{\kappa^2}{{( {\omega - 1} )}^2} \sin^2\theta }}{{[ {1 + {{( {\omega \,{-}\, 1} )}^2} + 2( {\omega \,{-}\, 1} )\cos \theta } ]{M_\theta }^2}}}, \end{gather}$$

(3.14)$$\begin{gather} {F_2}( \theta ) = \frac{{{M_\theta }\sqrt {1 + {{( {\omega \,{-}\, 1} )}^2} + 2( {\omega \,{-}\, 1} )\cos \theta } }}{{\kappa ( {\omega - 1} )}}\sqrt {1 \,{-}\, \frac{{{\kappa^2}{{( {\omega \,{-}\, 1} )}^2}{{\sin }^2}\theta }}{{[ {1 + {{( {\omega \,{-}\, 1} )}^2} + 2( {\omega \,{-}\, 1} )\cos \theta }]{M_\theta }^2}}}, \end{gather}$$

(3.15)$$\begin{gather} {F_3}( \theta ) = 3\cos \theta - \frac{{4{\kappa^2}{{( {\omega - 1} )}^2}\sin^2\theta \cos \theta }}{{[ {1 + {{( {\omega - 1} )}^2} + 2( {\omega - 1} )\cos \theta } ]{M_\theta }^2}}. \end{gather}$$

A maximum limiting value of $\lambda _\theta$ exists. Accordingly, the position of $P_{cav}$, as well as the focus point of the one-time reflected rays, is obtained, which has the following expression:

(3.16)\begin{equation} \frac{{{\lambda_{max}}}}{{{R_D}}} = \lim_{\theta \to 0} \frac{{{\lambda_\theta }}}{{{R_D}}} = \frac{1}{{\displaystyle\lim_{\theta \to 0} {F_1}( \theta ) - \lim_{\theta \to 0} {F_2}( \theta ) + \lim_{\theta \to 0} {F_3}( \theta )}} = \frac{\kappa }{{3\kappa - {M_0}{f_s}}}, \end{equation}

where ${f_s} = {\omega /(\omega - 1)}$ is the function characterising the effect of shock wave shapes in the case of the cylindrical converged incident shock wave.

The above expression has indicated that the relative distance from the cavitation kernel point $P_{cav}$ to the column centre is determined by three dimensionless parameters: the shock Mach number $M_0$, the radius ratio $\omega$ between the cylindrical converged shock wave and the water column and the sound-speed ratio $\kappa$ between the water and the air. If the radius ratio tends to the infinite, the cylindrical converged shock wave degenerates into a planar one. Here, we define a new dimensionless parameter $n = \kappa /{M_0}$, which represents the ratio of the wave speed between the transmitted shock wave and the incident shock wave. The distance $x_f$ from the focus point $P_{cav}$ to the centre of the water column can be written as

(3.17)\begin{equation} {x_f} = \frac{n}{{3n - 1}}{R_D}, \end{equation}

which is then identical to the expression reported by Biasiori-Poulanges & El-Rabii (Reference Biasiori-Poulanges and El-Rabii2021). In other words, the theoretical analysis derived in this paper is universal for different incident shock wave shapes.

3.3. The second convergence of the second reflected waves

When Re-RW$_{IC}$ completely focuses, the continuous diverged rarefaction wave (Re-RW$_{IIC}$) propagates upstream of the water column and is reflected by the column surface, forming the second reflected compression wave (Re$_2$-CW) and the second reflected rarefaction wave (Re$_2$-RW). Due to the specific column surface, the one branch of the second reflected wave (Re$_2$-CW) firstly focuses on the central axis of the water column, $P_{tran}$, causing a high-pressure region whose maximum pressure reaches more than 60 times the initial pressure, as shown in figure 7($k$). After an extremely short time, the two symmetrical reflected rarefaction waves (Re$_2$-RW), located on the upper and lower sides of Re$_2$-CW, merge and cause a rapid decrease in the local pressure. The pressure distribution along the centre axis of the water column at six different time instants, before and after the focus of Re$_2$-CW, is shown in figure 14. Obviously, the maximum pressure in the high-pressure region increases rapidly as the focus process continues before the complete focussing of Re$_2$-CW. After that, a significant negative-pressure region appears inside the water column, and the maximum pressure in the high-pressure region decreases rapidly and then disappears.

Figure 14.

The pressure distribution along the centre axis of the water column at six different time instants for the interaction between the cylindrical converged shock and the water column.

The evolution characteristics of the second reflected wave from $t_4$ ($t^* = 2.0$) to $t_5$ ($t^* = 2.326$) based on the ray analysis are shown in figure 15($a$). Due to the different properties of the two branches of the second reflected wave, the liquid near $P_{tran}$ undergoes a violent pressure oscillation in a short time interval. The highly transient characteristics of liquid pressure have great significance in the biomedical field, which is a valuable measure to remove diseased tissue. Similar to the derivation of $P_{cav}$, the specific position of $P_{tran}$ can be obtained by the ray analysis. A schematic diagram of rays emitted from some contact points when Re$_2$-CW focuses is shown in figure 15($b$), selecting an arbitrary ray (the black lines) for further analysis. Denoting the intersection points of the two-time reflected part of this ray with the column surface and the central axis of the water column as $E_{\theta 2}$ and $F_{\theta 2}$, respectively, with the distance from $F_{\theta 2}$ to $C$ being $\delta _{\theta }$,

(3.18)\begin{equation} {\delta_\theta } = \frac{{\sin {\beta_\theta }}}{{\sin ({\rm \pi} - 5{\beta_\theta } + \theta )}}{R_D} = \frac{{\cos {\alpha_\theta }}}{{\cos ( {5{\alpha_\theta } + \theta } )}}{R_D}. \end{equation}

Figure 15.

($a$) Schematic diagram of the second reflected wave propagation from $t_4$ ($t^* = 2.0$) to $t_5$ ($t^* = 2.326$) and schematic diagram of shapes and positions of the reflected waves at time $t_5$; ($b$) schematic diagram of ray analysis for the intersection point between the reflected rays, with two-time reflection, and the central axis of the water column.

After substituting (3.6) into (3.18) and simplifying, the minimum distance of $\delta _\theta$ can be obtained as $\delta _{min}$, which has the following expression:

(3.19)\begin{align} \frac{{{\delta_{min}}}}{{{R_D}}} &= \lim_{\theta \to 0} \left[ {\frac{{\cos {\alpha_\theta }}}{{( {16\cos^4{\alpha_\theta } - 20\cos^2{\alpha_\theta } + 5} )\cos {\alpha_\theta }\cos \theta - ( {16\cos^4{\alpha_\theta } - 12\cos ^2{\alpha_\theta } + 1} )\sin {\alpha_\theta }\sin \theta }}} \right]\nonumber\\ &= \frac{\kappa }{{5\kappa - {M_0}{f_s}}}. \end{align}

Moreover, to ensure that the two-time reflected rays can intersect with the horizontal central axis before it is reflected from the column surface, the emission angle $\alpha _\theta$ and the contact angle $\theta$ will satisfy: $6{\alpha _\theta } + \theta < {360^ \circ }$. It is easy to find that, as the reflection time increases, the maximum contact angle of rays participating in the focus process gradually decreases. This means that fewer rays participate in the focusing process with the interaction process continuing, which makes the reflected wave and focus intensity weaker and weaker.

Similar to the transmitted shock wave and the first reflected rarefaction wave, the second reflected waves also reflect on the column surface, forming the third reflected compression wave, figure 15($a$), which will focus inside the water column causing a high-pressure region. As the interaction continues, the high-pressure and low-pressure areas appear alternately at different positions on the centreline. Their analysis procedure is similar to the focusing of the first reflected expansion wave and the second reflected compression wave. The following expression gives a general formula of the distance from the focusing point $x_{f,N}$ of the $N$ times reflected wave to the column centre:

(3.20)\begin{equation} {\frac{{{x_{f,N}}}}{{{R_D}}} = \frac{\kappa }{{( {2N + 1} )\kappa - {M_0}{\,f_s}}}} ,\quad {N = 1,2,3, \ldots }. \end{equation}

However, the intensity of these waves, with three times or more reflection, has been dramatically reduced compared with the intensity of the transmitted shock wave, the first reflected rarefaction wave and the second reflected wave. Therefore, we will not elaborate further in the present study.