2.1. Governing equations and constitutive models

We consider a 3-D spatial domain comprised of multiple subdomains occupied by different fluid and solid materials. Figure 1 presents a schematic drawing of the problem set-up. A spherical gas bubble with a radius of $R_{0}$ is placed in a liquid at a standoff distance $D_{0}$ from the solid material (measured from the centre of the bubble). An incident shock wave propagates towards the bubble from the opposite side, causing its collapse. Let $\varOmega _{G}$, $\varOmega _{L}$ and $\varOmega _{S}$ denote the subdomains occupied by the bubble, the ambient fluid (liquid) and the solid material, respectively. We model both the gas inside the bubble and the ambient fluid as compressible inviscid fluids. Dropping of viscosity can be justified by the high pressure of the incident shock wave in the current problem set-up (see, e.g. Brujan (Reference Brujan2019) for a different and more complex fluid environment) and a relatively small time frame considered within this study.

Figure 1.

A 3-D fluid–solid coupled model of shock-induced bubble collapse near a solid or soft material.

Therefore, in $\varOmega _{G}$ and $\varOmega _{L}$, we solve the following Euler equations in the Eulerian frame:

(2.1)\begin{equation} \frac{\partial \boldsymbol{W}(\boldsymbol{x},t)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \mathcal{F}(\boldsymbol{W}) = 0, \quad \forall \, \boldsymbol{x} \in \varOmega_{L}(t) \cup \varOmega_{G}(t), \ t>0,\end{equation}

with

(2.2a,b)\begin{equation} \boldsymbol{W}=\begin{bmatrix} \rho \\ \rho \boldsymbol{V} \\ \rho e_{t} \end{bmatrix}, \quad \mathcal{F} = \begin{bmatrix} \rho \boldsymbol{V}^T \\ \rho \boldsymbol{V} \otimes \boldsymbol{V} + p \boldsymbol{\boldsymbol{\mathsf{I}}} \\ (\rho e_{t} + p)\boldsymbol{V}^T \end{bmatrix}, \end{equation}

where $\rho$, $\boldsymbol{V}$, $p$ and $e_{t}$ denote the fluid density, velocity, pressure and the total energy per unit mass, respectively; $\boldsymbol{\boldsymbol{\mathsf{I}}}$ is the $3\times 3$ identity matrix. Equation (2.1) is closed by an equation of state (EOS) for each fluid material. We apply the perfect gas EOS to the gas inside the bubble, i.e.

(2.3)\begin{equation} p=( \gamma_{G} -1 ) \rho e ,\end{equation}

where $\gamma _{G}$ is the heat capacity ratio, and $e$ is the internal energy per unit mass. In this work, we set $\gamma _{G}=1.4$. We assume the ambient fluid is water, and apply the stiffened equation of state, i.e.

(2.4)\begin{equation} p = ( \gamma_{L}-1 ) \rho e - \gamma_{L} p_{L},\end{equation}

where $\gamma _{L}=6.59$ and $p_{L}=410$ MPa (Johnsen & Colonius Reference Johnsen and Colonius2008).

Within the subdomain of the solid material, i.e. $\varOmega _S$, we adopt the Lagrangian frame and solve the following equation of motion that enforces the balance of linear momentum:

(2.5)\begin{equation} \rho_{s} \ddot{\boldsymbol{u}}(\boldsymbol{X},t) - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} (\boldsymbol{u},\dot{\boldsymbol{u}}) = \boldsymbol{b}, \quad \forall \, \boldsymbol{X} \in \varOmega_{S}(0), \ t >0, \end{equation}

where $\boldsymbol {u}$ denotes displacement, $\rho _{s}$ the material's mass density and $\boldsymbol {\sigma }$ the Cauchy stress tensor; $\boldsymbol {b}$ denotes the body force acting in $\varOmega _{S}$, which is assumed to be zero in this study.

We adopt the generalized Maxwell model of viscoelasticity (Cheneler Reference Cheneler2016), which combines linear springs and linear dashpots. Assuming zero strain at $t=0$, the stress at any time instance $t>0$ can be written in the form of the hereditary integral,

(2.6)\begin{equation} \boldsymbol{\sigma}(t) = \int^{t}_{0} G_{R}(t-s) \frac{\textrm{d} \boldsymbol{\epsilon}(s)}{\textrm{d} s}\ \textrm{d} s, \end{equation}

where $\boldsymbol{\epsilon}$ denotes the strain tensor; $G_{R}(t)$ is the relaxation modulus. In this work, we express $G_{R}(t)$ as a Prony series, i.e.

(2.7)\begin{equation} G_{R}(t) = G_{0} \left[ 1 - \sum^{N}_{i=1} g_{i} ( 1 - \textrm{e}^{-t/\tau_{i}} ) \right], \end{equation}

where $G_{0}=G_{R}(0)$ is the instantaneous modulus, $g_{i}$ and $\tau _{i}$ are material-specific model parameters. The Prony series is appealing particularly from a numerical viewpoint, as it does not require storing of the time history of strain. Specifically, we compute $\sigma (t)$ using a recursive formula, as shown in Goh, Charalambides & Williams (Reference Goh, Charalambides and Williams2004).

For some of the materials investigated in this study (e.g. a strengthened gypsum, metals), the viscoelastic effect of the material is negligible. These materials are modelled essentially by linear elasticity, with constitutive equation

(2.8)\begin{equation} \boldsymbol{\sigma} = \frac{E\nu}{(1+\nu)(1-2\nu)} tr(\boldsymbol{\boldsymbol{\epsilon}})\boldsymbol{\boldsymbol{\mathsf{I}}} + \frac{E}{1+\nu} \boldsymbol{\epsilon}, \end{equation}

where $E$ and $\nu$ are the Young's modulus and Poisson's ratio of the material, respectively; $tr(\cdot )$ denotes the trace operator.

At the fluid–solid interface, $\varGamma _{FS} = \partial \varOmega _{S}(t) \cap (\partial \varOmega _{L}(t) \cup \partial \varOmega _{G}(t))$, we enforce the continuity of normal velocity and surface traction, i.e.

(2.9)\begin{equation} \left. \begin{array}{l} (\boldsymbol{V} - \dot{\boldsymbol{u}})\boldsymbol{\cdot} \boldsymbol{n} = 0, \\ -p\boldsymbol{n} = \boldsymbol{\sigma}(\boldsymbol{u},\dot{\boldsymbol{u}})\boldsymbol{\cdot} \boldsymbol{n}, \end{array}\right. \quad \hbox{on } \varGamma_{FS},\end{equation}

where $\boldsymbol {n}$ denotes the unit normal to $\varGamma _{FS}$.

At the liquid–gas interface, $\varGamma _{FF}(t) = \partial \varOmega _{L}(t) \cap \partial \varOmega _{G}(t)$, we assume that the two fluid materials are immiscible. Also, surface tension is negligible compared to the pressure of the prescribed shock wave. Thus, the interface conditions are those describing a contact discontinuity, i.e.

(2.10)\begin{equation} \left. \begin{array}{l} (\boldsymbol{V}_{L}-\boldsymbol{V}_{G}) \boldsymbol{\cdot} \boldsymbol{n} = 0, \\ p_{L} = p_{G}, \end{array} \right. \quad \hbox{on }\varGamma_{FF}.\end{equation}

2.2. Fluid–solid coupled computational framework

We apply a recently developed fluid–solid coupled computational framework to solve the above model equations. The framework couples a finite volume, multiphase compressible flow solver with a finite element solid dynamics solver using an embedded boundary method and a staggered partitioned procedure (figure 2).

Figure 2.

Illustration of the spatial (a) and temporal (b) discretization methods applied in the computational framework.

As illustrated in figure 2(a), the fluid governing equations are semi-discretized in an augmented fluid domain, $\tilde {\varOmega } = \varOmega _L \cup \varOmega _G \cup \varOmega _S$, using an unstructured, node centred, non-interface-conforming finite volume mesh, denoted by $\tilde {\varOmega }^{h}$. For each node of the mesh (e.g. node $i$), a dual cell (i.e. control volume) $C_{i}$ is constructed. Integrating (2.1) within $C_i$ gives

(2.11)\begin{equation} \frac{\partial \boldsymbol{W}_i}{\partial t} + \frac{1}{\lVert C_i \rVert} \sum_{j\in Nei(i)} \int_{\partial C_{ij}} \mathcal{F}(\boldsymbol{W})\boldsymbol{\cdot} \boldsymbol{n}_{ij} \ \textrm{d} S = 0, \end{equation}

where $\boldsymbol{W}_i$ denotes the average of $\boldsymbol{W}$ in $C_i$, $\lVert C_i \rVert$ denotes the volume of $C_i$, $Nei(i)$ denotes the set of nodes connected to node $i$ by an edge, $\partial C_{ij} = \partial C_i \cap \partial C_j$ and $\boldsymbol {n}_{ij}$ is the unit normal to $\partial C_{ij}$. Depending on the location of nodes $i$ and $j$, four scenarios arise in the calculation of the surface integral over $\partial C_{ij}$:

1. (a) If nodes $i$ and $j$ belong to the same fluid subdomain ($\varOmega _L$ or $\varOmega _G$), the flux across $\partial C_{ij}$ is calculated using the well-known monotonic upwind scheme conservation law (MUSCL) scheme (Van Leer Reference Van Leer1979) and Roe's flux (Roe Reference Roe1981).

2. (b) If nodes $i$ and $j$ belong to different fluid subdomains (i.e. one in $\varOmega _L$, the other in $\varOmega _G$), a one-dimensional (1-D) two-fluid Riemann problem is constructed along edge $i$–$j$, i.e.

   (2.12)\begin{equation} \frac{\partial \boldsymbol{w}}{\partial \tau} + \frac{\partial \mathcal{F}(\boldsymbol{w})}{\partial \xi} = 0, \quad \text{with } \boldsymbol{w}(\xi,0)=\begin{cases} \boldsymbol{w}_{i} & \text{if } \xi \leq 0,\\ \boldsymbol{w}_{j} & \text{if } \xi > 0, \end{cases} \end{equation}
   where $\tau$ denotes the time coordinate; $\xi$ denotes the spatial coordinate along the 1-D axis aligned with $\boldsymbol {n}_{ij}$ and centred at the midpoint between nodes $i$ and $j$. The initial states $\boldsymbol {w}_{i}$ and $\boldsymbol {w}_{j}$ are projections of $\boldsymbol{W}_{i}$ and $\boldsymbol{W}_{j}$ on the $\xi$ axis. Equations (2.10) are enforced at the moving interface for $\tau >0$. The exact solution of this 1-D Riemann problem is supplied to Roe's flux function to compute the flux across $\partial C_{ij}$ (Farhat et al. Reference Farhat, Gerbeau and Rallu2012).

3. (c) If one of the two nodes (e.g. node $i$) belongs to a fluid subdomain ($\varOmega _L$ or $\varOmega _G$), while the other node belongs to the solid subdomain ($\varOmega _S$), a 1-D fluid–solid Riemann problem with a moving wall boundary is constructed, i.e.

   (2.13)\begin{gather} \frac{\partial {\boldsymbol{w}}}{\partial \tau} + \frac{\partial {\mathcal{F}}({\boldsymbol{w}})}{\partial \xi} = 0, \quad \tau>0,\ \xi < v_S \tau, \end{gather}
   (2.14)\begin{gather}{\boldsymbol{w}}(\xi,0) = \boldsymbol{w}_i, \quad \xi<0, \end{gather}
   (2.15)\begin{gather}v(v_S\tau,\tau) = v_S, \quad \tau>0, \end{gather}
   where $\xi$ is the spatial coordinate along the normal direction of the fluid–solid interface and centred at the midpoint between nodes $i$ and $j$. The initial state $\boldsymbol {w}_i$ is the reconstructed fluid state at the interface; $v$ denotes the fluid velocity at the moving wall boundary, and $v_S$ the normal velocity of the wetted surface of the solid structure (i.e. the fluid–solid interface) at its intersection with fluid edge $i$–$j$, which is computed by the solid dynamics solver. Similar to the previous case, the exact solution of this 1-D Riemann problem is supplied to Roe's flux function to compute the flux across $\partial C_{ij}$ (Wang et al. Reference Wang, Rallu, Gerbeau and Farhat2011).

4. (d) If nodes $i$ and $j$ both belong to the solid subdomain ($\varOmega _S$), the flux across $\partial C_{ij}$ is set to $0$.

The above method is referred to as FIVER (Farhat et al. Reference Farhat, Gerbeau and Rallu2012). It is ideally suited for solving the present problem, as the problem involves both multi-fluid (i.e. liquid–gas) and fluid–solid interfaces.

FIVER requires tracking of the liquid–gas and fluid–solid interfaces in the unstructured, non-interface-conforming mesh $\tilde {\varOmega }^{h}$. We track the evolution of the liquid–gas interface by solving the level-set equation,

(2.16)\begin{equation} \frac{\partial \phi(\boldsymbol{x},t)}{\partial t} + \boldsymbol{V} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi = 0, \quad \forall \, \boldsymbol{x}\in \varOmega_L \cup \varOmega_G \cup \varOmega_S, \end{equation}

where $\phi (\boldsymbol {x},t)$ represents the signed shortest distance from $\boldsymbol {x}$ to the interface. In this way, the large deformation and topological changes (e.g. splitting and merging) of the bubble surface are naturally accommodated. To track the fluid–solid interface, we apply a collision-based computational geometry algorithm as presented in Wang et al. (Reference Wang, Grétarsson, Main and Farhat2012, Reference Wang, Lea and Farhat2015).

In this work, (2.16) is recast in the form

(2.17)\begin{equation} \frac{\partial \phi}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\phi \boldsymbol{V}) = \phi\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{V}, \end{equation}

and solved on the same fluid mesh using a finite volume method. The convection term, $\boldsymbol {\nabla } \boldsymbol {\cdot } (\phi \boldsymbol{V})$, is discretized using the MUSCL scheme with the upwinding flux. The term $\phi \boldsymbol {\nabla } \boldsymbol {\cdot} \boldsymbol{V}$ is added as a source term. Additional numerical details can be found in Main et al. (Reference Main, Zeng, Avery and Farhat2017) and Main (Reference Main2014).

A standard Galerkin finite element method is used to semi-discretize the weak form of (2.5), which yields

(2.18)\begin{equation} {\boldsymbol{\mathsf{M}}}\frac{\partial^2 \boldsymbol{u}^h}{\partial t^2} + \boldsymbol{f}^{int}\left(\boldsymbol{u}^h, \frac{\partial\boldsymbol{u}^h}{\partial t}\right) = \boldsymbol{f}^{ext},\end{equation}

where ${\boldsymbol{\mathsf{M}}}$ denotes the mass matrix, $\boldsymbol {u}^h$ denotes the discrete displacement vector; $\boldsymbol {f}^{int}$ and $\boldsymbol {f}^{ext}$ denote the discrete internal force and external force vector, respectively.

We use a staggered fluid–solid time integrator presented in Farhat et al. (Reference Farhat, Rallu, Wang and Belytschko2010) to advance the fluid and solid subsystems. The fluid subsystem is time integrated using the explicit fourth-order Runge–Kutta scheme, while the solid subsystem is time integrated using the second-order central difference scheme. Notably, the fluid and solid time steps are offset by half a step (figure 2b). This is a designed feature to achieve second-order accuracy in time while maintaining optimal numerical stability.

3.1. Wave propagation across a plane fluid–solid interface

We first apply the computational framework described in § 2.2 to solve a 3-D model problem featuring the propagation of an impulsive pressure wave across a planar fluid–solid interface. The objective is to verify the computational framework for predicting wave transmission and deflection at a fluid–solid interface, which is an important feature of the bubble–material interaction problem under investigation. Figure 3(a) presents the set-up of the problem. We consider a cylindrical fluid–solid domain with a diameter of $12$ mm and a height of $20$ mm. In the fluid subdomain, we impose a spherical pressure wave generated from a monopole source, similar to the shock wave emitted from the rapid collapse of a bubble. Specifically, we specify a Ricker wavelet, which is widely used as a model of the seismic waves from a point excitation. Table 1 shows the geometric and material information involved in the simulation.

Figure 3.

(a) Problem set-up (the dynamic pressure of the fluid at $t = 0$ is visualized). (b) Snapshots of numerical solution at four different time instances, showing the dynamic pressure of the fluid and the maximum principal stress in the solid.

Table 1.

Parameters of the verification problem.

$c_{L}$: longitudinal wave speed; $c_{T}$: transverse wave speed.

We prescribe the Ricker wavelet as the initial condition of the fluid. Specifically, the initial pressure at a distance $r$ from the source is given by

(3.1)\begin{equation} p(r) = p_{0} + \frac{Q(-r/c_{0})}{r}, \end{equation}

with

(3.2)\begin{equation} Q(s) = [ 1-2{\rm \pi}^{2}f_{0}^{2}(s+s_{1})^2 ] \exp({-{\rm \pi}^2f_{0}^2 (s+s_{1})^2}), \end{equation}

where $p_{0}$ denotes the hydrostatic pressure, $c_{0}$ the speed of sound in the fluid, $f_{0}$ the frequency of Ricker wavelet and $s_{1}$ is a constant that controls the initial position of the wavelet. The values of these parameters are given in table 1. The initial velocity of the fluid is then prescribed according to the acoustic theory. The initial pressure wave has a peak magnitude of $913$ Pa and a width of ${\sim }2$ mm in the radial directions (figure 3a). Its central peak is located at ${\sim }1.31$ mm above the material surface.

The simulation is carried out using an unstructured, finite volume fluid mesh with $21.9$ million tetrahedral elements and a finite element solid mesh with $6.7$ million tetrahedral elements. In the most refined region, the characteristic element size is approximately $0.01$ mm for the fluid and $0.04$ mm for the solid. A constant time step size of $5.0 \times 10^{-4}\ \mathrm {\mu }$s is used in both solvers.

Figure 3(b) presents the numerical solution at four time instances, showing the fluid pressure field and the maximum principal stress in the solid. The simulation captures the transient dynamics of wave propagation across the fluid–solid interface, including the reflection of the pressure wave and the elastodynamic response of the solid. More specifically, figure 4 shows different types of waves that occur at a specific time instance, including the incident and the reflected waves in the fluid, the longitudinal (or $P$) and transverse (or $S$) waves in the solid, and the head waves generated by the advancing of the transmitted $P$ wave along the interface.

Figure 4.

Velocity field at $t = 3\ \mathrm {\mu }$s. The incident and reflected waves in the fluid, the transmitted waves in the solid and the head wave generated at the interface are clearly captured.

In this problem, the maximum dynamic pressure in the fluid is two orders of magnitude smaller than the hydrostatic pressure. Therefore, the incident wave can be considered as a small disturbance and modelled as an acoustic wave propagating in a homogeneous fluid medium. Given the assumption that the solid material is isotropic and linear elastic, this problem can thus be modelled adequately by coupling the linear acoustic wave equation with the equation of motion for a linear elastic solid. The exact solution of this problem can be obtained using the Cagniard–de Hoop method (de Hoop & Van der Hijden Reference de Hoop and Van der Hijden1984), therefore providing a reliable reference for verifying our numerical solver. We use an open-source code, Gar6more3D (INRIA 2013), to compute the exact solution. Figure 5 presents a comparison of the numerical and the analytical solutions at three sensors located in the fluid and in the solid. A close agreement is obtained. For example, the peak dynamic pressure recorded at sensor ($-$1.4, 1.05, 0.35) mm differs by only $0.66\,\%$ between the numerical and the analytical solutions.

Figure 5.

Comparison of the numerical solution obtained using the CFD (computational fluid dynamics) – CSD (computational solid dynamics) coupled solver with the analytical solution. (a) Time history of pressure at sensor ($-$1.4, 1.05, 0.35) mm in the fluid. (b–d) Time history of velocity at sensors (4, 0, 0) mm and (1.4, 1.05, 0.35) mm in the solid.

3.2. Collapse of a spherical bubble in an infinite medium

To demonstrate the numerical model's capability of capturing single bubble dynamics, we simulate the collapse of a spherical bubble in an infinite liquid medium. The simulation is based on the experiment of Kröninger et al. (Reference Kröninger, Köhler, Kurz and Lauterborn2010) in which the dynamics of a laser-generated spherical bubble in water is captured using high-speed photography. We focus on simulating the initial collapse of the bubble from its maximum radius $R_{0}=0.747$ mm. We assume that when the bubble is at its maximum size, velocity is zero both inside and outside of the bubble, and the thermodynamic variables have constant values within each subdomain. Specifically, the initial density and pressure of liquid water are set to $1000$ kg m$^{-3}$ and $0.101$ MPa, respectively; and the stiffened equation of state described in § 2.1 is applied. The composition and state of the gas inside the bubble was not measured in the experiment. In the simulation, we set the gas material to be air, modelled using the perfect gas equation of state ((2.3)). In previous computational studies, the initial state of the bubble is often estimated via theoretical models or data fitting (e.g. Müller et al. Reference Müller, Bachmann, Kröninger, Kurz and Helluy2009; Koch et al. Reference Koch, Lechner, Reuter, Köhler, Mettin and Lauterborn2016) due to the lack of experimental data. Although the time history of bubble radius in those studies shows good agreement with the experimental data, the initial pressure and temperature applied therein are often unrealistically small (e.g. $p_{0} = 4.579\ \text {Pa}$ in Müller et al. Reference Müller, Bachmann, Kröninger, Kurz and Helluy2009). Here, we apply the initial density calculated in Müller et al. (Reference Müller, Bachmann, Kröninger, Kurz and Helluy2009) (i.e. $\rho _{0}=0.957\times 10^{-3}$ kg m$^{-3}$) but use a more realistic initial pressure, $p_{0} = 100$ Pa. Leveraging the spherical symmetry in this problem, the computational domain is set to be a tall and slender tetrahedron whose apex coincides with the centre of the bubble, with symmetric boundary condition applied at its three side surfaces.

Five simulations are carried out using meshes of different resolutions. Specifically, the characteristic element size ($\varDelta h$) within and around the region of the bubble is refined from $60$ to $3.75\ \mathrm {\mu }$m.

Figure 6 presents the time evolution of bubble radius obtained using the five meshes. Evidently, as the mesh gets refined, the numerical result converges gradually to the experimental data for the period of bubble collapse. Figure 7 presents the velocity magnitude and the pressure field around bubble at four time instances for the converged solution. Driven by pressure difference, the bubble interface accelerates toward its interior and remains spherical during the collapse. When the bubble reaches its minimum volume, a strong shock wave is emitted.

Figure 6.

Comparison between numerical results and experimental data for the time history of bubble radius. Five sets of numerical result obtained with different mesh resolutions ($\varDelta h$) are presented (the radius is calculated based on the volume of the gas domain). The experimental data are extracted from Kröninger et al. (Reference Kröninger, Köhler, Kurz and Lauterborn2010).

Figure 7.

Converged solution at four time instances during collapse. Upper figure: pressure field in the tetrahedron computational domain. Lower figure: pressure and velocity magnitude in a 2-D plane across the bubble centre (the white curve represents the bubble surface). The full solution field is reconstructed by post-processing (reflection) the solution field on a plane across the apex of the tetrahedron domain.

It is noteworthy that for the time period of the bubble's rebound, the converged numerical result does not closely match the experimental data. This is likely due to the fact that the numerical model does not account for some physical processes (e.g. liquid/gas phase change) that become important when the bubble is at its minimum volume, as well as the inconsistent modelling of the bubble content (c.f. Koch et al. Reference Koch, Lechner, Reuter, Köhler, Mettin and Lauterborn2016). However, since the main focus of this work is the initial collapse of the bubble, the agreement with the experimental data presented above provides support for the validity of the numerical model.

4.1. Set-up of numerical experiment

In this section, we consider shock-induced bubble collapse near three representative materials, namely a strengthened gypsum (BegoStone 15 : 3, Esch et al. Reference Esch, Simmons, Sankin, Cocks, Preminger and Zhong2010; Simmons et al. Reference Simmons, Cocks, Zhong and Preminger2010), an elastomer (polyurea P1000, Amirkhizi et al. Reference Amirkhizi, Isaacs, McGee and Nemat-Nasser2006; Chahine & Hsiao Reference Chahine and Hsiao2015) and a foam rubber (Rubatex R8702 styrene-butadiene rubber (SBR) foam, Deigan Reference Deigan2007). Table 2 summarizes the main properties of these materials. In particular, the relative acoustic impedance of the material compared to that of the ambient fluid (i.e. $Z/Z_{0}$) varies from $5.2$ in the case of BegoStone to $1.1$ in the case of polyurea and $0.2$ in the case of the foam. Both polyurea and the SBR foam are modelled as a viscoelastic material. For the BegoStone, viscoelastic effects are negligible within the small time frame considered in this study, therefore it is modelled as a purely elastic material. As a reference, we also simulate the same bubble collapsing near a rigid wall, which can be considered as a material with infinite acoustic impedance.

Table 2.

Material properties of three representative materials.

Figure 8 shows the set-up of the numerical experiment. Initially, a spherical air bubble with a radius of $R_{0}=0.05$ mm is placed in the liquid, at a distance of $D_{0}=2R_{0}$ from the material surface (measured from the centre of the bubble). The incident shock wave that induces the bubble's collapse is adopted from Johnsen & Colonius (Reference Johnsen and Colonius2008). It has a compressive front with a positive peak of $35$ MPa, followed by a tensile phase with a negative peak of $10.1$ MPa (figure 8). This type of waveform can be generated in water through electrohydraulic or electromagnetic mechanisms (e.g. Church Reference Church1989; Fovargue et al. Reference Fovargue, Mitran, Smith, Sankin, Simmons and Zhong2013). We prescribe the incident shock wave as the initial condition of the fluid governing equations, following the method described in Cao et al. (Reference Cao, Zhang, Liao, Zhong and Wang2019). Away from the incident shock wave, the ambient flow velocity and hydrostatic pressure are set to $v_0 =0$ mm s$^{-1}$ and $p_{0}=0.101$ MPa, respectively. The 3-D model described in § 2.1 is applied. Notably, this problem can also be simulated using a 2-D axisymmetric model, which may be computationally less expensive. Nonetheless, the 3-D model can be easily extended in future to investigate more complicated cases without rotational symmetry, such as problems involving multi-bubble interaction and misalignment between the shock wave and the bubble. For computational efficiency, we consider a $45^{\circ }$ slice of the cylindrical domain with symmetry boundary conditions applied to the two cut planes, instead of modelling the entire 3-D domain. The computational domains for the solid and the fluid are discretized using $2.3$ million and $19.6$ million tetrahedron elements, respectively. The meshes are refined within a region close to the bubble and the fluid–solid interface, where the characteristic element size is $\varDelta h_{f} = 1.5\ \mathrm {\mu }\text {m}$ for the fluid mesh and $\varDelta h_{s} = 3\ \mathrm {\mu }\text {m}$ for the solid mesh. For comparison, the compressive front of the incident shock wave is approximately $700\ \mathrm {\mu }\text {m}$ in the liquid.

Figure 8.

Set-up of numerical experiment. A shock wave with a magnitude of $35$ MPa is applied to trigger the collapse of an air bubble with radius $R_{0}=0.05$ mm located at $D_{0}=2R_{0}$ from the material surface. Computations are performed on a $45^{\circ }$ slice of the cylindrical fluid–solid domain.

5.1. Acoustic impedance as a key parameter

The results presented in § 4 show that the acoustic impedance of the solid material has significant impact on the bubble and fluid dynamics. A complexity, however, is that acoustic impedance is a function of the material's density and elastic properties. Specifically (see, e.g. Brekhovskikh & Godin Reference Brekhovskikh and Godin2012),

(5.1)\begin{equation} Z = \sqrt{\frac{\rho_s E(1-\nu)}{(1+\nu)(1-2\nu)}}.\end{equation}

Therefore, it is impossible to vary acoustic impedance while keeping all the other properties the same, or to fix acoustic impedance while varying only one of the other properties. Indeed, because the materials considered in § 4 are selected from real-world materials of practical importance, they differ not only in $Z$, but also in $\rho _s$, $E$ and $\nu$.

To gain more insight into the role of acoustic impedance, we consider three numerical tests that involve artificial materials with prescribed properties. In each test, we fix the value of $Z$, and vary two of the three parameters: $\rho _s$, $E$ and $\nu$.

1. (i) Test 1: shock-induced bubble collapse near a series of materials that have the same values of $Z$ and $\rho _s$, but different values of $E$ and $\nu$.

2. (ii) Test 2: shock-induced bubble collapse near a series of materials that have the same values of $Z$ and $\nu$, but different values of $\rho _s$ and $E$.

3. (iii) Test 3: shock-induced bubble collapse near a series of materials that have the same values of $Z$ and $E$, but different values of $\rho _s$ and $\nu$.

In Test 1, we consider four materials as shown in table 4, including BegoStone. Their Young's modulus and Poisson's ratio vary by a factor of $16$ and $1.8$, respectively. For each material, we perform the numerical analysis described in § 4.1, using the same computational model.

Table 4.

Test 1 material properties.

Table 5 summarizes the results obtained from Test 1, including the bubble's collapse time ($t_{c}$), the minimum bubble volume ($V_{min}/V_{0}$), the maximum pressure at the sensor shown in figure 13 and the maximum values of maximum principal and shear stresses at the sensor shown in figure 16. The table shows that when $Z$ and $\rho _s$ are fixed, the bubble and fluid dynamics are relatively insensitive to variation in $E$ or $\nu$. The small changes can be explained by the fact that while the material's density and longitudinal wave speed are fixed, the transverse wave speed does vary among the four materials. For the same reason, it is not surprising that the stress field inside the material clearly varies among the four cases.

Table 5.

Test 1 results.

MPS and MSS: the temporal maximum of maximum principal stress and maximum shear stress at the sensor shown in figure 16. For MPS, the negative (i.e. compressive) peak is shown.

In the same way, we perform Test 2 for the materials shown in table 6, in which $\rho _s$ varies from $1500$ to $10\ 000\ \text {kg}\ \text {m}^{-3}$, and $E$ from $20.2$ to $134.4\ \text {GPa}$. The result shows that $t_c$ varies by less than $1\,\%$, from $0.298$ to $0.300\ \mathrm {\mu }\text {s}$; $V_{min}/V_{0}$ varies by around $3\,\%$, from $0.00275$ to $0.00285$. The maximum sensor pressure varies by around $5\,\%$, from $280$ to $296$ MPa. In comparison, variation in the material's stress field is much larger. For example, at the aforementioned sensor location (figure 16), the peak value of maximum shear stress varies from $121.4$ to $200.1\ \text {MPa}$.

Table 6.

Test 2 material properties.

In Test 3, we consider the four materials shown in table 7. The material's acoustic impedance and Young's modulus are fixed, and their density and Poisson's ratio vary by a factor of $7.5$ and $1.2$, respectively. Again, the fluid results are relatively insensitive to the material variations. For example, the peak pressure at the aforementioned sensor location varies by ${\sim }5\,\%$, from $119$ to $125\ \text {MPa}$.

Table 7.

Test 3 material properties.

Combining the results of the above numerical tests and those presented in § 4, it is evident that under the simulated conditions, the solid material's acoustic impedance ($Z$) is an important parameter that controls the material's reciprocal effect to the two-phase fluid flow. The same bubble collapsing near materials with different values of $Z$ can produce clearly different bubble and fluid dynamics, with $Z=Z_0$ being a transition point. On the other hand, the same bubble collapsing near materials that have the same value of $Z$, yet differ in other elastic properties, are found to produce similar bubble and fluid dynamics (yet dissimilar stress field inside the material).

5.2. Selection of materials with different values of acoustic impedance

To characterize the effect of acoustic impedance on bubble and fluid dynamics, we select five additional materials with acoustic impedance varying in a broader range compared to those considered in § 4. The materials are selected from several broad categories, as shown in figure 17. Table 8 summarizes the properties of these materials. In particular, their acoustic impedance varies from $0.063$ to $17.28$ MPa$\cdot$s m$^{-1}$. For each material, we conduct the same numerical analysis described in § 4.1.

Figure 17.

Selected materials for parametric study (the Young's modulus–density diagram is adapted from Ashby Reference Ashby2010).

Table 8.

Properties of five additional materials with different acoustic impedance.

In the following, we combine the new materials in table 8 with those considered in § 4 (table 2), and characterize the effect of acoustic impedance on three fluid quantities of interest, namely the bubble's collapse time, the minimum bubble volume and the maximum hydrodynamic pressure on the material surface.

5.3. Effect of acoustic impedance on bubble collapse time

Figure 18 shows in black dots the bubble's collapse time, $t_c$, for the eight materials considered. As the material's acoustic impedance increases, $t_c$ decreases monotonically. This trend is consistent with what we found in § 4 using the three representative materials, and has been explained there. The dashed line in figure 18 marks $Z/Z_0=1$. For materials with $Z/Z_{0}>1$, the variation in $t_c$ is relatively small. For example, as $Z/Z_{0}$ decreases from $10.8$ (aluminium) to $1.1$ (polyurea), $t_c$ increases by $10\,\%$. When $Z/Z_{0}<1$, $t_c$ becomes more sensitive to $Z$. It increases from $0.36\ \mathrm {\mu }\text {s}$ in the case of the rubber ($Z/Z_{0}=0.66$) to $0.465\ \mathrm {\mu }\text {s}$ in the case of the SBR foam ($Z/Z_{0}=0.2$).

Figure 18.

Effect of material's acoustic impedance on the bubble's collapse time. Dashed line: $Z=Z_0$.

The trend shown in figure 18 can be captured using a simple model. We recall that in the case of free-field Rayleigh collapse modelled by the Rayleigh–Plesset equation, after neglecting the gas content inside the bubble, the collapse time can be derived as (Brennen Reference Brennen2014)

(5.2)\begin{equation} t_{c}^{R} = 0.915 \sqrt{\frac{\rho_{l}}{p_{\infty}-p_{v}}} R_{0},\end{equation}

in which the superscript R denotes Rayleigh collapse; $\rho _l$ denotes the density of the liquid, $p_\infty$ the far field pressure, $p_v$ the vapor pressure and $R_0$ the initial bubble radius. Further, for the non-spherical Rayleigh collapse of a bubble near a rigid wall, Rattray (Reference Rattray1951) proposed that the bubble's collapse time can be approximated by

(5.3)\begin{equation} t_{c}^{R,w} \approx t_{c}^{R} \left( 1+0.205 \frac{R_{0}}{D_{0}} \right) = 0.915 \sqrt{\frac{\rho_{l}}{p_{\infty}-p_{v}}} R_{0} \left( 1+0.205 \frac{R_{0}}{D_{0}} \right), \end{equation}

where the superscript w on $t_c$ denotes wall, and $D_0$ is the distance between the centre of the bubble to the wall (cf. Johnsen & Colonius Reference Johnsen and Colonius2009; Koch et al. Reference Koch, Lechner, Reuter, Köhler, Mettin and Lauterborn2016).

In the present study, the bubble's collapse is induced by a shock wave, which is different from a uniform pressure $p_\infty$ in Raleigh collapse. Moreover, in § 4 we have shown that the bubble's collapse time depends sensitively on both the incident shock and its reflection against the material surface, while the amplitude of the reflected wave is a function of $Z/Z_0$. Therefore, we replace $p_\infty$ in (5.3) by

(5.4)\begin{equation} p_s = p_i + c_r p_r = p_i \left(1 + c_r \frac{Z/Z_0 - 1}{Z/Z_0 + 1}\right), \end{equation}

where $p_i$ is a characteristic pressure of the incident shock, $p_r$ the pressure of its reflection ((4.1)), and $c_r$ a non-dimensional coefficient that accounts for the time lag between the impact of the incident shock wave and that of the reflected wave. In other words, we propose that for the current study, the bubble's collapse time can be described by

(5.5)\begin{equation} t_{c} = 0.915 \sqrt{\frac{\rho_{l}}{p_{i} \left( 1 + c_{r}\dfrac{Z/Z_0 - 1}{Z/Z_{0} + 1} \right) - p_{v}}} R_{0} \left( 1+0.205\dfrac{R_{0}}{D_{0}} \right).\end{equation}

The solid line in figure 18 shows the prediction of the above model, with $p_i = 20\ \text {MPa}$ (approximately half of the peak pressure of the incident shock) and $c_r=0.7$. It agrees well with the result of the numerical analysis.

5.4. Effect of acoustic impedance on minimum bubble volume

In figure 19, the black dots show the bubble's minimum volume normalized by its initial volume, i.e. $V_{min}/V_{0}$, for the eight materials considered. As the material's acoustic impedance increases, $V_{min}/V_0$ decreases monotonically. The trend is nonlinear and again, $Z/Z_0=1$ roughly marks a turning point. The minimum bubble volume is more sensitive to the nearby material's acoustic impedance when $Z/Z_{0} < 1$. For example, $V_{min}/V_{0}$ in the case of polypropylene foam ($Z/Z_{0} = 0.04$) is more than $9$ times higher than that in the case of rubber ($Z/Z_{0} = 0.66$). In comparison, the difference between the case of Ultracal-30 stone ($Z/Z_{0} = 3.4$) and the case of aluminium ($Z/Z_{0} = 10.8$) is less than $20\,\%$.

Figure 19.

Effect of material's acoustic impedance on the bubble's minimal volume. Dashed line: $Z=Z_0$.

The effect of acoustic impedance on minimum bubble volume can also be captured using a phenomenological model derived from the Rayleigh–Plesset equation. Again, we recall that in the case of free-field Rayleigh collapse modelled by the Rayleigh–Plesset equation, the minimum bubble volume can be derived analytically (Brennen Reference Brennen2014). After neglecting surface tension, we have

(5.6)\begin{equation} R^{R}_{min} = R_{0} \left[ \frac{p_{G_0}}{(\gamma_{G}-1) \left( p_{\infty}-p_{v} \right)} \right]^{1/[3(\gamma_{G}-1)]},\end{equation}

where $R_{min}^R$ denotes the minimum radius of the bubble, which is assumed to remain spherical all the time. Again, the superscript $R$ denotes Rayleigh collapse; $p_{G_0}$ denotes the partial pressure of gas inside the bubble; $p_v$ denotes the vapor pressure; and $\gamma _G$ the specific heat ratio of the internal gas.

To account for non-spherical collapse, we generalize the definition of bubble radius ($R$) as

(5.7)\begin{equation} R = \left( \frac{V}{4 {\rm \pi}/3} \right)^{1/3}. \end{equation}

Then, we generalize the above analytical result from Rayleigh–Plesset equation in the same way described in § 5.3, which gives

(5.8)\begin{equation} R_{min} = R_{0} c_{s} \left[ \frac{p_{0}}{(\gamma_{G}-1)\left[ p_{i} \left( 1+c_{r}\displaystyle \frac{Z-Z_{0}}{Z+Z_{0}} \right) - p_{v} \right]} \right]^{1/[3(\gamma_{G}-1)]}.\end{equation}

The solid line in figure 19 shows the minimum bubble volume predicted by (5.8) with $p_i=20\ \text {MPa}$, $c_{r}=0.7$ and $c_{s}=7.0$. The prediction agrees well with the result of numerical analysis.

5.5. Effect on the maximum pressure on material surface

Figure 20 shows in black dots the maximum hydrodynamic pressure on the material surface, denoted by $p_{{w,max}}$, for the eight materials considered. Again, a monotonic relationship is observed. Specifically, as $Z/Z_0$ varies from $0.04$ (in the case of the polypropylene foam) to $10.8$ (in the case of aluminium 7075), $p_{{w,max}}$ increases by $26$ times, from $11.3$ to $291.6$ MPa. The increase of maximum pressure is nonlinear with respect to $Z/Z_0$, and the slope reduces as the acoustic impedance becomes larger (i.e. concave). It can be expected that as $Z/Z_{0}$ goes to infinity, $p_{{max}}$ would approach the value obtained from the case with a rigid wall (cf. figure 13).

Figure 20.

Effect of material's acoustic impedance on the maximum hydrodynamic pressure on material surface. Dashed line: $Z=Z_0$.

For the eight materials considered (with $Z/Z_{0} \in [0.04,10.8]$), $p_{{w,max}}$ is approximately a linear function of $\log (Z/Z_{0})$. Using least-squares fitting, the relationship is given by

(5.9)\begin{equation} p_{{w,max}}^{(1)} = \rho c_{0}^2 \left[ 0.07 \log_{10}\left( \frac{Z}{Z_{0}} \right) + 0.0416 \right], \end{equation}

where $c_0$ denotes the speed of sound in the fluid. This fitting function is plotted in figure 20. Despite the simplicity, a linear function like (5.9) cannot capture the correct asymptotic behaviour as $Z/Z_0$ goes to infinity. Accounting for the case of $Z/Z_0=\infty$ (i.e. a rigid wall), we find that $p_{{w,max}}$ can also be fitted using a bounded growth function,

(5.10)\begin{equation} p_{{w,max}}^{(2)} = \rho c_{0}^2 \frac{0.0591(Z/Z_0)}{1+0.4434(Z/Z_0)},\end{equation}

which is also shown in figure 20.

5.6. Effect of bubble's stand-off distance

While the material's acoustic impedance dictates the reflection of shock wave at the material surface, the bubble's stand-off distance determines the time when the reflected shock waves reach and interact with the bubble, and hence, is an important parameter that affects the subsequent bubble dynamics. To study how the stand-off distance $D_{0}$ influences the material's reciprocal effects discussed above, we re-run the three representative cases presented in § 4 using different $D_{0}$ between $R_{0}$ and $6R_{0}$.

Table 9 compares the results obtained with a stand-off distance of $D_{0} = 4 R_{0}$ and the results for $D_{0} = 2R_{0}$ from table 3. For both stand-off distances, similar trends in bubble's collapse time and the peak surface pressure with respect to $Z/Z_{0}$ can be observed. However, the differences between three materials in the case of $D_{0} = 4R_{0}$ becomes smaller, indicating that the effect of nearby materials becomes weaker with increased stand-off distance. This is because it takes longer (i.e. approximately twice as the case of $D_{0} = 2R_{0}$) for the reflected wave to reach and interact with the collapsing bubble, which thus weakens the material's effect on the collapse process. Another interesting observation is that because of the reduced effect from the nearby material, a shock wave is emitted in the case of the SBR foam with $D_{0} = 4R_{0}$.

Table 9.

Comparison of the dynamics of bubble collapse with two stand-off distances ($D_{0} = 2 R_{0}$ and $D_{0} = 4 R_{0}$).

Furthermore, there exists a limit value of standoff distance, above which the bubble would reach its minimum volume before the reflection of the incident shock arrives, therefore its initial collapse is not affected by the solid material. Our numerical test suggests that for the current setting, the limit stand-off distance is approximately $D^{*}_{0} = 5.5R_{0}$.

We also find that (5.5) provides a good estimate of this limit stand-off distance. Specifically, the time for the shock wave to travel between the bubble and the material surface (back and forth) is given by

(5.11)\begin{equation} t_{wp} = \frac{2D_{0}}{c} = 2D_{0} \sqrt{\frac{\rho_{l}}{\gamma_{L} (p+p_{L})}} \approx 2D_{0} \sqrt{\frac{\rho_{l}}{\gamma_{L} p_{L}}},\end{equation}

where $D_{0}$ denotes the bubble's stand-off distance, $\gamma _{L}$ and $p_{L}$ are constant coefficients of the stiffened equation of state. At the limit stand-off distance, denoted by $D^{*}_{0}$, $t_{wp}$ is equal to the bubble's collapse time, $t_{c}$. Combining (5.5) with $c_{r} = c^{*}_{r} = 0$ and (5.11) gives

(5.12)\begin{equation} 2D^{*}_{0} \sqrt{\frac{\rho_{l}}{\gamma_{L} p_{L}}} = 0.915 \sqrt{\frac{\rho_{l}}{p_{i}-p_{v}}} R_{0} \left( 1+0.205 \frac{R_{0}}{D^{*}_{0}} \right) \end{equation}

which is a quadratic equation of $R_{0}/D^{*}_{0}$. Solving this equation gives $R_{0}/D^{*}_{0} = 0.1809$, or equivalently, $D^{*}_{0}=5.53 R_{0}$, which matches well the result of our numerical test.