# Wall pressure amplification of shock-collapsed multi-bubble arrays near a rigid wall

Eric Goncalves da Silva, Philippe Parnaudeau

► To cite this version:

Eric Goncalves da Silva, Philippe Parnaudeau. Wall pressure amplification of shock-collapsed multi-bubble arrays near a rigid wall. Physics of Fluids, 2025, 37 (1), pp.013319. 10.1063/5.0246108. hal-04872038

HAL Id: hal-04872038

https://hal.science/hal-04872038v1

Submitted on 7 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Wall pressure amplification of shock-collapsed multi-bubble arrays near a rigid wall

Eric Goncalves da Silva¹, a) and Philippe Parnaudeau¹

Institut Pprime, UPR 3346 CNRS, ISAE-ENSMA, 1 avenue Clément Ader,
86961 Futuroscope Chasseneuil cedex, France.

This numerical study investigates the collapse of various arrangements of gas bubbles im-
mersed in water in the vicinity of a rigid wall and impacted by a planar shock wave. Mul-
tiple bubble configurations, from 2 to 5 bubbles, are compared, focusing primarily on the
pressure loads on the wall and the potential amplification in comparison with the single-
bubble case. The three-dimensional simulations are performed using a massively parallel
compressible diffuse interface solver. The effects of the grid resolution and the mass trans-
fer term are discussed. The main characteristics of the flows are described and the dynamic
behaviors in pressure wave propagation are illustrated. A power-law is proposed for the
evolution of the maximum pressure peak on the wall as a function of the density ratio of the
bubble array. An amplification of a factor 30 is highlighted for a pyramidal arrangement.

Keywords: Multi-bubble collapse; shock wave; pressure load; compressible two-phase model; high performance computing

<footnote>a) Electronic mail: eric.goncalves@ensma.fr</footnote>
# I. INTRODUCTION

The collapse of gas bubbles in a liquid is a violent, intense, and brief phenomenon that can cause significant damage when occurring near a wall, such as the well-known case of ship propeller erosion¹, or serve therapeutic purposes, such as in kidney stone treatment². This collapse is driven by the pressure difference between the gas inside the bubble and the surrounding liquid. The pressure field distribution in the liquid determines the type of bubble collapse, such as:

• Symmetric (or spherical), when the pressure field is uniform,

• Asymmetric (or non-spherical), when the pressure field is non-uniform (gradient).

Both phenomena occur naturally and have been studied for over a century.

Lord Rayleigh³ provided one of the first comprehensive analyses and models of symmetric bubble collapse, assuming an infinite, incompressible, and inviscid fluid, while neglecting surface tension. Since then, several researchers, including Plesset⁴, Gilmore⁵, and Keller and Miksis⁶, have proposed more advanced models incorporating viscosity, surface tension, and compressibility. In symmetric collapses, the bubble contracts due to inertial forces, generating extremely high internal pressures before collapse. This overall behavior also applies in the case of asymmetric collapse. Kornfeld and Suvorov⁷ were the first to suggest the existence of asymmetric collapse, referring to "unstable surface cavities," though they did not illustrate them. Walters and Davidson⁸,⁹ later conducted theoretical and experimental studies on bubble collapse, occasionally observing non-spherical shapes. They also detailed an experimental method for generating such bubbles, characterized by a liquid jet or "liquid tongue" forming inside the bubble (ring or toroidal shape). Benjamin and Ellis¹⁰ provided a thorough theoretical explanation, attributing non-spherical collapse to the presence of a pressure gradient in the liquid. Another hypothesis involves the collapse being triggered by a planar shock wave. Later studies, such as those by Shima, Tomita, and Takahashi¹¹ and Tomita and Shima¹², focused on the shock-induced collapse of bubbles near a wall, primarily examining erosion effects. Bourne and Field¹³ were the first to comprehensively investigate the phenomena occurring during shock-induced single-bubble collapse, exploring various parameters such as the intensity of the pressure field, and the shape and size of the cavities.

Due to the complexity of simulating multiphase flows, numerical studies on shock-induced bubble collapse in fluids were conducted later. Ball et al.¹⁴ developed a numerical tool using the Free Lagrangian Method, validating their simulations with the work of Bourne and Field¹³.
Nourgaliev, Dinh, and Theofanous¹⁵ proposed an Adaptive Mesh Refinement (AMR) method coupled with a Level-Set method for this purpose. While both studies¹⁴,¹⁵ made significant strides in accurately describing the gas-liquid interface, they remained two-dimensional (2D) and lacked an in-depth analysis of the underlying physics. Johnsen and Colonius¹⁶ were among the first to use numerical simulations to analyze bubble collapse dynamics, both in free fields and near walls. They employed diffuse interface modeling with high-order spatial discretization, revisiting Shima, Tomita, and Takahashi¹¹'s work on wall erosion, though limiting their study to weak shocks. The diffuse interface method's requirement for high-resolution at the liquid-gas interface, coupled with precision schemes, resulted in high computational costs, particularly for three-dimensional (3D) simulations. However, Coralic and Colonius¹⁷ improved upon the work of Johnsen and Colonius¹⁶, achieving satisfactory results, including in 3D. The first 3D simulations of shock-induced gas (air) bubble collapse in liquid (water) were introduced by Hawker and Ventikos¹⁸, who employed a front-tracking method to model strong shock wave interactions (with a pressure jump of 1 GPa), providing detailed analysis of the collapse dynamics. Since then, other authors have conducted 3D numerical studies. Some focused on validating new numerical approaches¹⁹⁻²¹ or comparing models²², while others investigated specific aspects of bubble collapse dynamics near walls²³. The collapse of bubble arrays under shock waves has also been studied²⁴⁻²⁶.

Indeed, the academic interest in studying the collapse of an isolated bubble, with or without the presence of a wall, is obvious. In nature, however, the erosion produced by such an event involves an array or a cloud of bubbles. Preliminary experimental work²⁷ has shown that a higher density of bubbles in the fluid leads to an increase of the erosion rate. Dear and Field²⁸ propose, using a new experimental technique, a parametric study of the collapse of arrays of air cavities (2D) within a flow of a mixture of water and gelatin. These authors study the influence of the number of cavities and the type of arrangement on the pressure released by each of these arrays. A plane incident shock wave hits an array of tree cavities (arranged horizontally or vertically with respect to the incident wave), a triangular array of six cavities or an square array of nine cavities. Their main conclusion is that the collapse of such arrangements always leads to an amplification of the peaks pressure. They explain a chain collapse mechanism as "the pressure waves from the first collapsed layer collapse the next and so on".

It was not until the work of Lauer *et al.*²⁴ that the first simulations of the collapse of multiple cavities by an incident shock wave were carried out. These authors investigate the collapse of a horizontal and then a vertical array of tree air cavities placed upstream of an incident shock wave
in water. They show that the successive collapses of cavities in the context of a horizontal array lead to increasingly large successive pressure peaks, as identified in²⁸. However, this amplification is only observed when the cavities are close enough together, as in²⁷. In addition, these authors, as well as²⁸, observe negative pressures in the intercavity zones and show that (in their context) the influence of a cavitation model does not affect the temporal evolution of the maximum pressure. In agreement with¹⁸, who explain that "rotationally symmetric and fully three-dimensional calculations yielded almost indistinguishable results", Betney et al.²⁵ perform simulations using the same numerical tool and revisit²⁸ triangular arrangement by playing with the size of the bubbles. Finally, they point out that the pressure peaks increase significantly under topological conditions. The collapse dynamics of an hemispherical cloud composed of 50 air bubbles close to a wall was studied in²⁹. Authors found that the peak pressure was associated with the centremost bubble, producing a corresponding peak pressure on the wall. To demonstrate the power of their numerical tool, Wermelinger et al.³⁰ present 3D simulations of different bubble arrays without giving details about the physics. Rasthofer et al.³¹, using the same tool, provide a detailed study of the collapse of a cloud of 12,500 bubbles. Bempedelis and Ventikos²⁶, based on the work of Betney et al.²⁵, attempt to measure the energy that can be released when a bubble array collapses. They study different arrays with different bubble sizes and show that an array can produce a pressure peak up to more than 100 times greater than the intensity of the generating incident wave. In addition, Goncalves and Parnaureau³² propose a study of the collapse of two bubbles placed horizontally in front of a wall. The conclusion is that the smaller the distance between the bubbles, the greater the amplification of the pressure load on the wall, in agreement with observations²⁷.

This study, carried out exclusively with 3D simulations, focuses on the pressure in water and, in particular on the wall due to a shock-induced collapse of bubble arrays placed in front of a wall. The knowledge of the pressure load on the wall will give us an idea of the potential damage it is subject to, which is critical in many applications. To the best of our knowledge, this is the first 3D numerical study to specifically address this question through a non-symmetric collapse of a bubble array. First, the study focuses on a grid-independent sensitivity of a strong (1 GPa) shock-induced collapse of an isolated bubble immersed in water, in the free field and then close to a wall. To achieve this, metrics are being defined based on the work of Hawker and Ventikos¹⁸ and comparisons with reference works¹⁸,²¹ are proposed. The same metrics are then used to investigate the influence of mesh size on the shock-induced collapse of a multi-bubble array, including a comparison with previous works²⁵,²⁶. Based on these preliminary conclusions, a parametric study
of the shock-induced collapse of a series of multi-bubble arrangements located near a wall is provided. The paper is structured as follows. The first section presents the model and the solver used to perform the 3D simulations presented in the rest of the article. The second section concerns the mesh dependency for various kinds of collapse. The third part is the main part and is about the parametric study of the 1 GPa shock-induced collapse of multi-bubble arrays near a wall. An empirical law is proposed to estimate the maximum wall pressure as a function of the arrangement density. Finally, future investigations are discussed.

## II. GOVERNING EQUATIONS AND METHODOLOGY

The hereafter simulations were performed with the in-house one-fluid solver SCB based on a 4-equation formulation³³, and massively parallelised using MPI and OpenACC, allowing the use of supercomputers with GPUs³⁴. Previous computations and comparisons with the literature on various shock-bubble interactions, as presented in²²,³²,³⁴, highlight the capability of SCB to accurately carry out the parametric analyses proposed in this study.

### A. Four-equation model

The main assumptions are that the mixture is assumed to be homogeneous and each phase is assumed to be inviscid and compressible. Within the temporal and spatial scales considered in this study, bubble collapse is primarily governed by inertial forces, with other mechanisms being negligible due to their minimal influence on the flow dynamics. Previous studies of shock-induced bubble collapse in free fields¹⁶,²⁵,²⁹ have demonstrated that relevant non-dimensional parameters, such as the Reynolds and Weber numbers, are exceedingly large, underscoring the dominance of inertial forces. More recently, when accounting for the presence of a wall, Goncalves and Parnaureau³² showed that viscous effects have little to no significant impact on the computation of pressure loads, even under quite similar physical parameters. Similarly, the Weber number (=$\frac{\rho u^2 D}{\sigma}$) exceeds 10⁶ (based on the jet velocity) indicating that surface tension effects are negligible compared to the dominant inertial force. Thus, diffusive effects and surface tension are neglected. Furthermore, mass transfer is also neglected, although specific work on the mass transfer hypothesis is discussed in sectionIV E. In the following, the subscripts *l* and *g* refer to the liquid and gaseous phases, respectively and the unindexed quantities refer to the mixture. The
void fraction and the mass fraction of the gas are $\alpha$ and $Y = \alpha\rho_g/\rho$. The expression of the model is:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0, \qquad (1) $$

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + P\mathbf{1}) = \mathbf{0}, \qquad (2) $$

$$ \frac{\partial \rho E}{\partial t} + \nabla \cdot (\rho H \mathbf{u}) = 0, \qquad (3) $$

$$ \frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = K \nabla \cdot \mathbf{u}, \quad \text{where} \quad K = \frac{\rho_l c_l^2 - \rho_g c_g^2}{\frac{\rho_l c_l^2}{1-\alpha} + \frac{\rho_g c_g^2}{\alpha}} \qquad (4) $$

with three conservation equations applied to the mixture: mass (1), momentum (2) and total energy (3), as well as a transport equation (4) for the void fraction with the source/sink term on the right side. Where $\rho = \alpha\rho_g + (1 - \alpha)\rho_l$ is the density, $\mathbf{u}$ the velocity vector, $P$ the pressure, $\mathbf{1}$ the identity tensor, $E$ the total energy and $H$ the total enthalpy. Furthermore, $\rho_l$, $\rho_g$, $c_l$ and $c_g$ denote respectively the density and sound velocity for each phase. The total energy is $E = e + \frac{1}{2}(\mathbf{u} \cdot \mathbf{u})$ with $e$ the internal energy and the total enthalpy is $H = h + \frac{1}{2}(\mathbf{u} \cdot \mathbf{u})$ with $h$ the enthalpy.

Pure phases follow the stiffened gas equation of state (EOS). From the thermal and mechanical equilibrium assumption, an expression for the mixture pressure and temperature can be deduced (see Goncalves and Parnaureau 32 for more details). In this study, the material parameters for the stiffened gas EOS are the same as those used in 26.

## B. Solver

The model (1)-(4) is discretised using the cell-centered finite volume method. The numerical flux and the non-conservative term are calculated using a HLLC scheme 35 and the second order (in space discretisation and time integration) is obtained by the MUSCL-Hancock extrapolation 36. A minmod slope limiter is also used. Finally, the time step is chosen to satisfy the Courant-Friedrichs-Lewy (CFL) criterion (fixed to 0.1 for all the simulations). Code parallelization is achieved through the hybrid use of the MPI library in combination with the OpenACC programming paradigm. Scalability, performance and benchmark validation are presented in Dubois, Goncalves, and Parnaureau 34.
### III. ASSESSMENT OF SPATIAL RESOLUTION

Before presenting the main part of this work, one must study the influence of three elements: domain size in transverse directions, mesh size and boundary conditions. The question of convergence study, in the context of inviscid flow equations, is a complex issue. Fortunately some metrics can address this. The section III A 1 is an overview of the incident shock wave induced collapse of a spherical air bubble in water in free field, using the same simulation parameters as in the literature¹⁸,²¹. In section III A 2, the non-reflecting outlet condition is replaced by a wall condition to study the mesh convergence in the presence of a solid wall.

#### A. 1 GPa shock-induced collapse of air bubble in water

According to previous studies¹⁸,²¹,²⁵, it is sufficient to treat only a quarter of the bubble, taking advantage of symmetries. The domain configuration is thus: An air bubble with a radius ($R = 0.003$ m) is placed in a domain of size $x \in [0, 6.6R]$, $y \in [0, 3R]$ and $z \in [0, 3R]$. At the start, the center of the bubble is defined as: $(x_b, y_b, z_b) = (4.6R, 0, 0)$. Initially the density, velocity and pressure of water and air are given in 5.

$$ (\rho, \mathbf{u}, p) = \begin{cases} (998 \text{ kg/m}^3, \mathbf{0} \text{ m/s}, 10^5 \text{ Pa}) & \text{Water} \\ (1.204 \text{ kg/m}^3, \mathbf{0} \text{ m/s}, 10^5 \text{ Pa}) & \text{Air} \end{cases} \qquad (5) $$

Finally, a post-shock condition is applied at position $x_{sh} = 3.1R$ with a pressure intensity $P_{sh} = 1$ GPa, corresponding to a Mach number $M_{sh} = 1.42$.

##### 1. Free-field case

Non-reflecting conditions are set in the inlet and outlet boundary conditions, while symmetric conditions are set in $y_0$, $z_0$ and reflective conditions in $y_{max}$ and $z_{max}$. Snapshots of the collapse are shown in figure 1. In the bottom plane the velocity magnitude contours are proposed. Perpendicularly to the bottom plane, Schlieren-like representations coloured by the dimensionless pressure ($P/P_{sh}$) and an iso-surface of the void are shown. The ratio $\alpha = 0.10$ is superimposed on these two maps. Simulations are performed using a uniform mesh corresponding to 200 points per bubble radius (ppbr).
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 1: Shock-induced collapse of an air bubble in water. Velocity magnitude, Schlieren-like representation colored by the dimensionless pressure ($P/P_{sh}$) and iso-surface of void ratio $\alpha = 0.10$ with ppbr = 200.</figcaption></figure>

The incident shock wave propagates with velocity $u_{sh} = 430 \ m.s^{-1}$ and hits the bubble at time $t = 0.650 \ \mu s$. One observes the deformation of the interface bubble, its acceleration (proof of the existence of the liquid jet at that moment), the reflected rarefaction wave propagating inside the liquid, while the transmitted shock propagates through the gas as an incident shock wave at time $t = 2.65 \ \mu s$. At time $t = 3.81 \ \mu s$, the liquid jet has just reached the other edge of the bubble. This is the definition of the end of the first phase¹⁸, which is the first metric used for the mesh convergence. The second metric, called jet speed at impact, corresponds to the velocity of the liquid jet at that time. In figure 1, a shock wave is emitted (commonly refered to as a water hammer shockwave) and its pressure intensity is the third metric. At time $t = 4.33 \ \mu s$, the main speed jet is then accelerated and the pressure on the central axis of the bubble increases significantly, until the liquid jet splits the bubble in two parts and bubble fragments are observed. At time $t = 5.00 \ \mu s$, the complex interaction of sheet-jetting shocks, reflected sheet jet impacts and the water-hammer leftward front produce the maximum pressure peak on the axis of symmetry in front of the bubble,
as well as the maximum velocity of the liquid jet. This is the fourth metric, referred to as the pressure of the interaction waves¹⁸. We underline that this peak pressure is of significant magnitude which exceeds the pressure of the water-hammer shock, as commented in¹⁸. All these metrics are reported in Table I.

The metrics of first phase time, jet speed at impact, and water hammer shockwave pressure indicate that a resolution of 200 ppbr is sufficient to capture the key features of collapse. However, a larger gap is observed when considering the pressure of the interaction waves. This pressure peak is associated to waves propagating and interacting in the reduced volume area of bubble lobes and therefore requires a very fine grid resolution. This indicates that the intensity of the peak is influenced by the size and spatial distribution of the bubble remnants, making it inherently challenging to determine the optimal initial resolution. Furthermore, discrepancies with Hawker and Ventikos¹⁸ may also arise from differences of EOS and numerical methods employed. Despite these challenges, the intensity of this peak remains the highest among all observed metrics, and its ratio to the water hammer peak is relatively consistent across both resolutions tested. Notably, previous two-dimensional studies¹⁸,²⁵,³⁷ have suggested that a minimum resolution of approximately 400 ppbr is required to achieve grid independence across all metrics. However, excluding the interaction wave pressure metric, our results show that most other metrics exhibit similar values at both 200 and 400 ppbr. Consequently, to balance accuracy and computational efficiency, we assume that a resolution of 200 ppbr is sufficient for capturing the overall dynamics and relative intensities of the phenomenon in this study, while keeping computational costs manageable.

<table><caption>TABLE I: 3D simulations of $p = 1$ GPa shock-induced collapse of an air bubble in water. Comparison of characteristic collapse metrics with Hawker and Ventikos¹⁸, Bempedelis and Ventikos²¹.</caption><thead><tr><th></th><th colspan="5">Present study</th><th>Bempedelis -3D</th><th>Hawker - 3D</th></tr></thead><tbody><tr><td>Resolutions (ppbr)</td><td>50</td><td>100</td><td>200</td><td>400</td><td>600</td><td>50</td><td>100</td></tr><tr><td>First phase time (t* = tc/R)</td><td>3.08</td><td>3.07</td><td>3.04</td><td>3.04</td><td>3.04</td><td>2.72</td><td>2.58</td></tr><tr><td>Jet speed at impact (m/s)</td><td>2815</td><td>2840</td><td>2840</td><td>2842</td><td>2847</td><td>2776</td><td>2754</td></tr><tr><td>Water hammer shockwave pressure (GPa)</td><td>4.18</td><td>4.32</td><td>4.61</td><td>4.62</td><td>4.97</td><td>5.464</td><td>4.04</td></tr><tr><td>Interaction waves pressure (GPa)</td><td>8.48</td><td>12.48</td><td>14.84</td><td>19.85</td><td>18.55</td><td></td><td>5.5</td></tr></tbody></table>

Table I shows also a comparison with previous results¹⁸,²¹ and a global agreement is observed. A discrepancy is noticeable for the dimensionless first phase time $t^* = tc/R$ where $c$ is the mixture speed of sound. It is certainly due to the difference of EOS used by authors. A similar shift
was reported in Hawker and Ventikos¹⁸ with the collapse time estimated by Nourgaliev, Dinh, and Theofanous¹⁵, whereas in Goncalves and Parnaudeau²², an agreement was observed with the same reference.

The choice of a slip boundary condition at $y_{max}$ and $z_{max}$ is driven by the primary objective of this study: to assess the loading pressure on the wall. This boundary condition effectively simulates the collapse of the bubble near a solid wall while preserving symmetry. While periodic boundary conditions could be considered as an alternative, they would significantly increase computational costs by a factor 4. Regarding the results of figure 2(a), one can observe that as soon as the distance between the center of the bubble and $y_{max}$ and $z_{max}$ is sufficient ($\ge 3R$), results are quite close. Furthermore, if the domain is too narrow ($L_y = L_z = 1.5R$), the effects are visible after the last pressure peak (e. g. when $t > 5.2 \mu s$).

<figure><img src="image_4.png" /><figcaption>FIG. 2: Time history of the maximum dimensionless pressure ($P/P_{sh}$) of the shock-induced bubble collapse.</figcaption></figure>

## 2. Collapse near a solid wall

The non-reflecting outlet condition is now replaced by a wall boundary condition. All other boundary conditions remain the same as in section III A 1. The bubble is detached from the wall and the distance between the centre of the bubble and the wall is set to $\Gamma = 2R$.

Jamaluddin et al.³⁸ proposed to simulate the shock-induced collapse of a single bubble near a wall, focusing on a comparison with the free-field case. They concluded that the presence of the wall appears to elongate the bubble during the collapse, and also increases the velocity of the liquid jet, ultimately producing a more powerful pressure peak than for the free-field case. In more
recent works, Goncalves and Parnaureau³², Turangan et al.³⁹ shared these conclusions for a va-
riety of bubble sizes and the pressure intensities of the incident shock wave. The comparison of
the evolution of the maximum pressure for cases with and without wall is plotted in Figure 2(b)
and shows that the process is exactly the same for both cases before the water-hammer rightward
front is reflected by the wall and returns to the bubble lobes (up to time t = 5.20 µs). After this
first stage, the process is very different.

Snapshots of the collapse are shown in figure 3. In the lower plane, the velocity magnitude con-
tours are illustrated. In the plane perpendicular to the bottom plane, Schlieren-like representations
coloured by the dimensionless pressure (P/P<sub>sh</sub>) are presented, and an iso-surface of the void ra-
tio α = 0.10 is superimposed on these two planes. As commented previously, the pressure peak
of the interaction waves is illustrated at time t = 5.00 µs and corresponds to the first peak (I)
in figure 4(a). Later, the water-hammer front is reflected on the wall generating pressure loads.
Several secondary shock waves are emitted inside the bubble lobes (sheet-jetting shocks) forming
multiple complex interactions with the leftward water-hammer front and with themselves. Right-
ward front waves impact and reflect on the wall generating pressure peaks during the time interval
[5.3, 5.7] µs (peaks (1) and (2) in figure 4(a)). The collision of these reflected waves, the reflected
water-hammer front, the reflected incident shock with the bubble fragments produces the largest
pressure peak in water around time t = 6.20 µs (peak (II) around 28 GPa in figure 4(a)). At time
t = 7.40 µs, the impact of a rightward wave on the wall generates a strong pressure peak (peak
(4) in figure 4(a)), which is the most intense. A third pressure peak in water (peak (III) around
17.65 GPa in figure 4(a)) is highlighted at the same instant due to the recollapse of bubble rem-
nants by a leftward wave. These waves resulting of multiple interactions are illustrated in Figure 5
where an enlargement of Schlieren-like visualizations on the cutting plane z = 0 is plotted. At time
t = 6.90 µs, these waves are located close to the bubble ring and propagate in opposite directions.
At time t = 7.40 µs, the rightward front has impacted the wall and the leftward front continue its
propagation. A few moments later, the bubble remnants attach to the wall. The maximum pres-
sure records show that the number and intensity of peaks are greater in the presence of a wall, as
already mentioned in literature³²,³⁸,³⁹.

The time history of the pressure loads on the wall (Figure 4(a)) shows the presence of four
peaks. At time t = 4.5 µs, an increase in parietal pressure is observed when the incident shock
wave hits the rigid wall, followed by the water-hammer shockwave, producing pressure loads. The
wall pressure peaks (1) and (2), between times t = 5.4 µs and t = 5.70 µs, are associated with
<figure><img src="image_5.png" /><figcaption>FIG. 3: Shock-induced collapse of an air bubble in water near a solid wall. Velocity magnitude, Schlieren-like representation colored by the dimensionless pressure ($P/P_{sh}$), dimensionless wall pressure ($P/P_{sh}$) and iso-surface of void ratio $\alpha = 0.10$ with ppbr = 200.</figcaption></figure>

the impact of secondary waves on the wall. These first two parietal pressure peaks are about 3 times weaker than the peak due to the interaction of shock waves upstream of the bubble in water (peak (I) in Figure 4(a)). The third peak appears about 0.1 µs after peak (II), with an intensity of about 14 GPa. This peak is related to the reflected waves emitted during the collapse of the bubble remnants. Finally, the strongest peak (17.65 GPa) is associated with a shock impact on the wall.

Figure 4(b) shows the evolution of the maximum wall pressure $P_{wallmax}$ normalised by the value obtained on the finest grid, for each peak as a function of the grid resolution. One can observe that, as in the free-field case, acceptable results can be obtained for a resolution of 200 ppbr, with an offset to the finest grid of less than 10%, except for the highest peak. This last pressure peak associated to complex wave interactions generated by the recollapse of bubble pieces requires a very fine grid to be well captured, that induce prohibitive computation costs for a parametric study.
<figure><img src="image_6.png" /><figcaption>FIG. 4: Time history of the maximum dimensionless pressure ($P/P_{sh}$) of the shock-induced collapse of an air bubble in water near a solid wall.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 5: Shock-induced collapse of an air bubble in water near a solid wall. Enlargement of Schlieren-like visualizations with ppbr = 200.</figcaption></figure>

## B. Shock-induced collapse of a triangular air bubble array in water

Following a 2D study²⁵, Bempedelis and Ventikos²⁶ investigated the energy released by the collapse of bubble arrays in 3D. One of the aims of this study was to determine whether 3D instabilities exist. To answer this question, they performed simulations without any symmetry assumptions. They concluded that there was no 3D instability but this hypothesis had consequences
for their grid resolution. In addition, attention was given to bubble arrays with different bubble sizes, where the finest one had a marginal resolution (e. g. 30 ppbr). These authors first considered two arrays of three bubbles of different sizes and found that in the most sparse case there was almost no increase in collapse intensity compared to the collapse of a single bubble. In this subsection, this configuration is re-examined to determine the influence of the mesh. And, in order to limit the computational cost of the simulations, the symmetries are maintained on the basis of their conclusions.

The setup is as follows: $R_b$ (= 0.0012 m) is the radius of the finest bubble and $R_B$ is the radius of the two large bubbles, with $R_B = 2.5R_b$. The size of the domain in each direction is $x \in [0, 6.6R_B]$, $y \in [0, 4.3R_B]$ and $z \in [0, 4.3R_B]$. Thanks to the symmetries, only one half of a large bubble and a quarter of the finest bubble are computed: the first half (the larger-one) is placed at $(x_{B1}, y_{B1}, z_{B1}) = (4.16R_B, 1.5R_B, 0)$ and the quarter of the finest bubble is placed at $(x_b, y_b, z_b) = (5.16R_B, 0, 0)$, resulting in an interbubble distance equal to $1R_B$. A 2D scheme is provided in figure 6. The incident shock wave is placed at $x_{sh} = 2.66R_B$, with the same intensity and

<figure><img src="image_8.png" /><figcaption>FIG. 6: 2D-Schematic of the triangular array of three bubbles. The dashed lines are part of the domain that is ignored in the computation due to symmetry conditions.</figcaption></figure>

Mach number as in the previous cases. According to the previous sections, the resolution should be set to $ppbr_b = 200$ for the finest central bubble, but the cost of such simulation would be too high (greater than 15 billion cells on a regular mesh) for current computing resources. Therefore, in the
remainder of this section, the influence of the grid on the collapse dynamics of the smallest bubble will be investigated. Four grids were used for this purpose, with the resolution $ppbr_b = (30, 60, 80$ and $100)$ for the finest bubble, and $ppbr_b = (75, 150, 200$ and $250)$ for the largest bubble. The main dynamics of the collapse of the three-bubble array will be described below.

The incident shock wave reaches the largest bubbles at the same time as in the free-field case for a single bubble. The shock wave carries along the two bubbles, resulting in the formation of a main transverse liquid jet. The water-hammer shockwave is emitted with a similar pressure intensity as in the case of a single bubble, in agreement with previous results$^{21}$. The first pressure peak corresponds to the water-hammer shockwave of the two largest ones (peak (1) in Figure 7). As mentioned previously$^{18,24,25}$, the generation of negative pressure (green area in figure 8) is mainly due to the density of the arrangement. This negative pressure is the evidence for the cavitation phenomenon, which is not considered in this simulation (e. g. without mass transfer term). The rarefaction waves of the first largest bubbles merge and creates a low-pressure ahead of the central finest bubble, as observed by Lauer et al.$^{24}$. As this negative pressure builds up in front of the

<figure><img src="image_9.png" /><figcaption>FIG. 7: The effect of the grid resolution for the shock-induced collapse of a triangular array of three air bubbles in water: Evolution of the maximum pressure.</figcaption></figure>

<table><thead><tr><th>ppbr<sub>B</sub></th><th>75</th><th>150</th><th>200</th><th>250</th></tr><tr><th>ppbr<sub>b</sub></th><th>30</th><th>60</th><th>80</th><th>100</th></tr></thead><tbody><tr><th>Peak 1</th><td>4.28</td><td>4.48</td><td>4.54</td><td>4.58</td></tr><tr><th>Peak 2</th><td>8.22</td><td>10.08</td><td>10.41</td><td>10.29</td></tr><tr><th>Peak 3</th><td>7.16</td><td>10.07</td><td>11.36</td><td>11.89</td></tr></tbody></table>

small bubble, the incident shock wave continues to advance, causing it to collapse on its axis of symmetry. As the water-hammer shockwaves of the large bubbles are emitted ($t = 3.90$ µs), the finest bubble is deformed until the emission of its water hammer shockwave ($t = 4.80$ µs). The water-hammer shockwave of the finest bubble has the same pressure intensity as that of the large bubbles, but it is not observable on the maximum pressure timeline because the largest bubbles are simultaneously subjected to a very high intensity phenomenon on their respective lobes, called
transmitted expelled gas shocks¹⁸. These shocks hit the finest bubble just after the water-hammer shockwave of the latter is released, which produces a pressure peak not along its axis of symmetry, but on its sides at time $t = 5.00$ µs. This causes an acceleration of the deformation of the finest bubble. Two pressure peaks (one for each bubble) are observed on their central axis at time $t = 5.15$ µs (peak (2) in Figure 7) due to the interaction of different shock waves emitted in the previous instants. This peak can be clearly identified as the interaction waves pressure peak of the single bubble appearing on its central symmetry axis. Finally, at time $t = 5.25$ µs, on the central axis of the finest bubble, the largest pressure peak in water is observed (peak (3) in Figure 7).

<figure><img src="image_10.png" /><figcaption>FIG. 8: Shock-induced collapse of a triangular array of three air bubbles in water. Schlieren-like representation colored by the dimensionless pressure ($P/P_{sh}$) and iso-surface of void ratio $\alpha = 0.10$ with ppbr = 200. Green: Negative pressure.</figcaption></figure>

Peak (1) has a pressure intensity close to that of an isolated bubble and is slightly affected by the resolution. Turning now to the influence of the mesh for peak (2), two aspects should be highlighted. Firstly, the intensity of the peak is significantly lower than that of the isolated bubble (about 30%). The origin of this significant drop is related to the merge of the rarefaction waves, as explain above. Secondly, its intensity remains independent of the resolution if $ppbr_B > 150$. Finally, looking at peak (3), which corresponds to the pressure peak on the central axis of the smallest
bubble, figure 7 shows two opposite conclusions depending on the grid resolution. Indeed, at a marginal resolution ($ppbr_B = 75$), peak (3) is not the strongest, and worse, it is much weaker than that observed for an isolated bubble (about 40 %). However, as the resolution increases, this peak becomes the most intense, similar to the case of an isolated bubble. If the resolution of the smallest bubble is sufficient, one can observe an enhancement of the collapse intensity for the smallest bubble compared to what it is for an isolated bubble. On the other hand, if the resolution is marginal, one could conclude that the collapse intensity is attenuated. In agreement with Bempedelis and Ventikos²⁶, one can notice that this arrangement does not cause an amplification of the collapse, but rather to a chain reaction, which, depending on the grid resolution, may be understood as an attenuation of the collapse intensity of the smallest bubble. In conclusion, it appears that the shock-induced collapse of this triangular bubble array is quite sensitive to the grid resolution.

On the basis of this section, the study will focus on the energy released by the collapse of multi-bubble arrays and the corresponding wall pressure. To streamline the subsequent parametric study and reduce costs, the analysis will concentrate on two parameters: the number of bubbles and their spatial arrangement, while keeping the bubble sizes constant.

#### IV. SHOCK-INDUCED COLLAPSE OF ARRAYS OF AIR BUBBLES IN WATER NEAR A WALL

This section investigates the effect of the number and density of bubbles on the intensity of the energy released during their collapse. The primary objective is to understand the extent to which this collapse intensity can generate a pressure peak on the wall, potentially leading to damage. To achieve this, various bubble arrangements are analyzed, and distinct parameters are defined for each case. The first parameter is $Γ$, which represents the distance between the solid wall and the center of the nearest bubble.

The $Γ_x$, $Γ_y$ and $Γ_z$ parameters are defined as the distance between the centres of each bubble in the x, y and z directions, respectively. All these parameters are dimensionless by the radius of the air bubble (R). The post-shock condition is applied at position $x_{sh} = 2.6R$ with a pressure intensity $P_{sh} = 1$ GPa, corresponding to a Mach number of $M_{sh} = 1.42$ in water.
<table><caption>TABLE II: Bubble arrangement and domain parameters</caption><thead><tr><th>Case</th><th>L<sub>x</sub>/R</th><th>L<sub>y</sub>/R</th><th>L<sub>z</sub>/R</th><th>Γ/R</th><th>Γ<sub>x</sub>/R</th><th>Γ<sub>y</sub>/R</th><th>Γ<sub>z</sub>/R</th></tr></thead><tbody><tr><td>1B</td><td>6.6</td><td>3</td><td>3</td><td>2</td><td></td><td></td><td></td></tr><tr><td>2B-X-1</td><td>7.6</td><td>3</td><td>3</td><td>2</td><td>3</td><td></td><td></td></tr><tr><td>2B-X-2</td><td>8.6</td><td>3</td><td>3</td><td>2</td><td>4</td><td></td><td></td></tr><tr><td>2B-X-3</td><td>9.6</td><td>3</td><td>3</td><td>2</td><td>5</td><td></td><td></td></tr><tr><td>2B-Y-1</td><td>6.6</td><td>4.6</td><td>4.6</td><td>2</td><td></td><td>3</td><td></td></tr><tr><td>2B-Y-2</td><td>6.6</td><td>5</td><td>5</td><td>2</td><td></td><td>4</td><td></td></tr><tr><td>2B-Y-3</td><td>6.6</td><td>4.6</td><td>4.6</td><td>2</td><td></td><td>5</td><td></td></tr><tr><td>3B-XY-1</td><td>7.6</td><td>4.6</td><td>4.6</td><td>2</td><td>3</td><td>3</td><td></td></tr><tr><td>3B-XY-2</td><td>7.6</td><td>5</td><td>5</td><td>2</td><td>3</td><td>4</td><td></td></tr><tr><td>3B-XY-1i</td><td>7.6</td><td>4.6</td><td>4.6</td><td>2</td><td>3</td><td>3</td><td></td></tr><tr><td>3B-XY-2i</td><td>7.6</td><td>5</td><td>5</td><td>2</td><td>3</td><td>4</td><td></td></tr><tr><td>5B-XYZ-1</td><td>8.6</td><td>4.6</td><td>4.6</td><td>2</td><td>4</td><td>3</td><td>3</td></tr><tr><td>5B-XYZ-2</td><td>8.6</td><td>5</td><td>5</td><td>2</td><td>4</td><td>4</td><td>4</td></tr></tbody></table>

The full set of bubble array configurations and the dimensions of the computational domains are summarised in table II. All other boundary conditions remain the same as in section III A 2 and, following the conclusion of the section III, the resolution is set to ppbr = 200.

In a previous study³², the influence of the distance between twin bubbles in a horizontal array of bubbles positioned perpendicular to the direction of the incident shock wave (Γₓ) and in front of a wall was investigated. Considering only the case of detached bubbles, the higher pressure peak is observed for Γₓ = 3R. This previous study used a shock wave of lower intensity (Psh = 120 MPa). Therefore, the influence of the intensity of the incident shock wave is considered here by playing with different values of Γₓ: 2B-X-[1,2,3] cases. Then, an arrangement of two bubbles in a vertical column was performed with the 2B-Y-[1,2,3] cases. Next, two types of triangular arrays are studied, which can be understood as 2D arrays because they consist of three bubbles placed in the same plane (Figure 15). The first case corresponds to an isolated air bubble close to the incident shock wave and the following bubbles are superimposed, these are the 3B-XY-[1,2] cases. 3B-XY-[1i,2i] is a second type of triangular arrangement where the isolated bubble is close to the wall. Finally, pyramidal array consisting of five bubbles with an isolated bubble placed in front of the four following bubbles and close to the incident shock wave, are studied (cases 5B-XYZ-[1,2]).
### A. Horizontal array of two bubbles

The shock-induced collapse of a horizontal array of two bubbles has already been considered<sup>24,25,28,32,40</sup>. For the sake of clarity, the following analysis is proposed in two parts. The first one (Section IV A 1) begins when the incident shock wave reaches the first bubble and ends when the first shock waves are reflected from the wall. After the second part (Section IV A 2) begins.

#### 1. Before reflection

As already mentioned in Dear and Field<sup>28</sup> and confirmed by Lauer *et al.*<sup>24</sup>, the water-hammer shockwave of the first bubble always hits the second bubble first (see Figure 9(a)) and plays the role of a shock wave with a curved incidence. The water-hammer shockwave pressure peak followed by the interaction waves pressure peak, reached during the collapse of the first bubble, are similar to those of a single bubble collapse (Section III A 2), with a pressure intensity of 4.61 GPa and about 14.90 GPa, respectively (see Figure 10(a)).

The pressure in the interbubble region diminishes, revealing a negative pressure state. As noted by Dear and Field<sup>28</sup> and Lauer *et al.*<sup>24</sup>, cavitation occurs primarily where the bubbles are closest together (case 2B-X-1). When the water-hammer shockwave impacts the second bubble, it generates a rarefaction wave. Additionally, in the case 2B-X-1, the first bubble is drawn into the second one, as illustrated in Figure 9(b), leading to a pronounced deformation of the remaining portion of the first bubble. This phenomenon arises mainly from the proximity of the two bubbles, which creates a low-pressure region in the interbubble zone, as highlighted by Betney *et al.*<sup>25</sup>.

Examining the jet speed at impact during the collapse of the second bubble, it becomes evident that this speed decreases with increasing distance between the bubbles (e. g. $\Gamma_x$). Specifically, in the case 2B-X-1, the jet speed upon impact of the second bubble is comparable to that of the first bubble, monitored to 2850 m/s. In contrast, for the configuration 2B-X-3, the speed drops to 2650 m/s, as illustrated in Figure 10(a). For the shortest interbubble distance ($\Gamma_x = 3R$), the pressure exerted by the water-hammer shockwave (6.10 GPa) and the pressure from interaction waves (16.05 GPa) on the second bubble are significantly higher than those experienced by the first bubble, exceeding by approximately 30% and 8%, respectively. An increase in the interbubble distance means equal or smaller pressure peak for the second bubble collapse than for the first one (see Figure 10(a)). Thus, the pressure intensity of the water-hammer shockwave for the second
<figure><img src="image_11.png" /></figure>

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0246108

<figure><img src="image_12.png" /><figcaption>FIG. 9: Shock-induced collapse of a horizontal array of two air bubbles in water near a wall. Velocity magnitude, Schlieren-like representation colored by the dimensionless pressure ($P/P_{sh}$), dimensionless wall pressure ($P/P_{sh}$) and iso-surface of void ratio $\alpha = 0.10$.</figcaption></figure>

<footnote n="402">402 bubbles are 4.61 and 4.55 GPa, and the interaction waves pressure are 14.42 and 13.24 GPa for the cases 2B-X-2 and 2B-X-3, respectively. This dependence on ($\Gamma_x$) aligns with previous studies$^{25,32}$.</footnote>

<footnote n="404">404 Furthermore, the location of the significant pressure peak during the collapse of the second bubble</footnote>
<figure><img src="image_13.png" /><figcaption>FIG. 10: Time history of the maximum dimensionless pressure ($P/P_{sh}$) of a horizontal array of two bubbles in water near a wall.</figcaption></figure>

consistently occurs along the central axis of the interbubble region (see Figure 9(c)).

One can see a reflected wave on the upper boundary condition, which can interact with the bubble ring (cases 2B-X-2 and 2B-X-3). As discussed in section III A 1 about the size of the computational domain, the main conclusions are not modified if the domain is extended, but this leads to a significant increase in computation costs.

## 2. After reflection

After the first pressure peak on the symmetry axis has disappeared, the first shock waves are reflected on the wall and collide with the remnants of the second bubble, creating a new pressure peak on the central axis between the two bubbles at time $t=t = 10.16 \ \mu s$ (see Figure 11(a)). The distance between the bubbles significantly changes the intensity of the events after this second pressure peak. For $\Gamma_x = 3R$, a contraction of the first bubble remnants is observed and their recollapse leads to a strong peak pressure (83.01 GPa) in water, at time $t = 10.81 \ \mu s$ (see Figures 11(b) and 10(b)). Finally, the largest parietal pressure peak (48.73 GPa) at time $t = 11.48 \ \mu s$ (see Figures 11(c) and 10(b)) happens, on the central axis, when a shock wave resulting from numerous interactions impacts the wall. In fact, multiple waves propagate in all directions between the two bubbles, creating complex interactions. The recollapse of the second bubble generates leftward fronts interacting with reflected waves which impact the first bubble ring around time $t = 10.55 \ \mu s$.
<figure><img src="image_14.png" /></figure>

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0246108

<figure><img src="image_15.png" /><figcaption>FIG. 11: Shock-induced collapse of a horizontal array of two bubbles in water near a wall. Velocity magnitude, Schlieren-like representation colored by the dimensionless pressure ($P/P_{sh}$), dimensionless wall pressure ($P/P_{sh}$) and iso-surface of void ratio $\alpha = 0.10$.</figcaption></figure>

<footnote n="422">422 $\mu$s producing shock waves in the central torus of the bubble ring, as illustrated in Figure 12 at</footnote>

<footnote n="423">423 time $t = 10.89 \mu$s. Later, the rightward wave front from this collapse propagates towards both the</footnote>

<footnote n="424">424 second bubble ring and the wall (see Figure 12 at time $t = 11.14 \mu$s). By time $t = 11.48 \mu$s the</footnote>
<figure><img src="image_16.png" /><figcaption>FIG. 12: Shock-induced collapse of a horizontal array of two bubbles in water near a wall, case 2B-X-1.
Velocity magnitude, Schlieren-like representation colored by the dimensionless pressure ($P/P_{sh}$),
dimensionless wall pressure ($P/P_{sh}$) and iso-surface of void ratio $\alpha = 0.10$.</figcaption></figure>

This pattern changes as the distance between the bubbles increases ($\Gamma_x = 4R$). At time $t = 11.50$ µs the pressure peak on the central axis is weaker (21.20 GPa) which is divided by around a factor 4 compared to the previous case. This can be explained by the fact that the low-pressure in the interbubble area is less important than in the previous case. Finally, at time $t = 12.60$ µs, the maximum parietal pressure peak (62.75 GPa) is observed on the central axis (see Figures 11(c) and 10(b)) due to the impact of a strong shock wave, e. g. the same scenario as for $\Gamma_x = 3R$.

The impact of shock waves on the wall is always associated with the maximum parietal pressure peak, but this peak is not necessarily associated with the maximum pressure peak in the water. Finally, for $\Gamma_x = 5R$, the distance between the bubbles eliminates the amplification phenomena described above and tends towards a damping of the pressure peak on the wall (22.50 GPa) which is quite similar to the case of the isolated bubble.

The shock-induced collapse of a horizontal array of two bubbles near a wall can be summarised as follows:

* The collapse of the second bubble is always initiated by the water-hammer shockwave of the first bubble.

* Rarefaction waves emitted during the collapse of the second bubble create a low-pressure
region in the interbubble area and negative pressure could occur if the two bubbles were close enough together.

* The maximum pressure peak on the wall is always greater in comparison with the single bubble case.

* The maximum pressure peak in water does not match with the maximum pressure peak on the wall.

* The maximum pressure peak on the wall results from a shockwave generated by a cascade of wave interactions between the bubbles and the wall.

## B. Vertical array of two air bubbles

<figure><img src="image_17.png" /><figcaption>FIG. 13: Shock-induced collapse of a vertical array of two bubbles in water near a wall. Dimensionless longitudinal velocity component $u/u_{sh}$, Schlieren-like representation colored by the dimensionless pressure ($P/P_{sh}$), dimensionless wall pressure ($P/P_{sh}$) and iso-surface of void ratio $\alpha = 0.10$.</figcaption></figure>

Ball et al.¹⁴ carried out an experimental-numerical comparison on the interaction of a shock
wave with bubbles in the context of a vertical arrangement. Lauer et al.²⁴ investigated the discrepancy observed between experimental and numerical measurements by Ball et al.¹⁴ with numerical simulations of the shock-induced collapse of a vertical array of bubbles. Unfortunately, the two studies differ in the intensity of the incident shock wave ($P_{sh} = 1.9$ GPa for Ball et al.¹⁴, $P_{sh} = 260$ MPa for Lauer et al.²⁴). The intensity of the incident shock wave in the present study ($P_{sh} = 1.0$ GPa) is closer to that of the study of Ball et al., while maintaining a vertical array of two bubbles ($\Gamma_y \in [3, 4, 5]$R) close to the study of Lauer et al.²⁴ ($\Gamma_y = 4R$). Thanks to the symmetry conditions, the computational domain is limited to half of one bubble.

In this part, the evolution of the collapse is mostly identical for the three cases, so only the two most compact cases are shown (see Figure 13), where the differences are mostly significant. As Lauer et al.²⁴ mentioned, there is a significant pressure drop between the two bubbles due to merging of the rarefaction waves and this is particularly true when the two bubbles are the closest together: e. g. case 2B-Y-1. In this arrangement, this low-pressure region slows down the collapse of the bubbles and reduces their pressure intensity. Indeed, the second pressure peak found on the symmetry axis in front of the bubbles, corresponding to interaction waves (in figures 13 and 14, at time $t = 7.40$ µs), is only 12.00 GPa, while in the other two cases it is very close to that observed for an isolated bubble (e. g. 14.61 GPa for cases 2B-Y-2 and 2B-Y-3). Parietal pressure also shows this decrease. Thus, the maximum pressure peak on the wall is about 14.97 GPa for the case 2B-Y-1 at time $t = 8.80$ µs (see Figure 13), while it is 15.60 GPa for cases 2B-Y-2 and 2B-Y-3 at around time $t = 9.90$ µs. These parietal pressure peaks are observed when the shock waves emitted by the collapse of each bubble overlap, and are weaker (about 10%) than a single bubble collapse near a wall case. The most obvious interaction effect observed during the collapse of such a vertical arrangement is the fact that at the shortest distance between the bubbles ($\Gamma_y = 3R$) there is a tilt towards the symmetry plane, as previously observed in²⁴. This suction phenomenon is almost non-existent as the distance between the bubbles increases ($\Gamma_y = 4R$ or $5R$). This is due to the interaction of the two initial rarefaction waves and mainly explains the attenuation of the pressure peak of the interaction waves.

The following statements summarise of a shock-induced collapse of a vertical array of two bubbles near a wall:

* The interaction is only sensitive when the bubbles are close together.

* The rarefaction waves emitted during collapse create a low-pressure region, decreasing the
<figure><img src="image_18.png" /><figcaption>FIG. 14: Time history of the maximum dimensionless pressure ($P/P_{sh}$) of a shock-induced collapse of a vertical array of two air bubbles in water near a wall.</figcaption></figure>

pressure peak on the wall compared to the single bubble case.

• For the present 3D simulations, no negative pressure appears in the low-pressure region. However, a 2D simulation reaches the opposite conclusion, in agreement with previous work²⁴,²⁸.

• The maximum pressure peak on the wall will always be less intense than in the case of a single bubble.

### C. Triangular arrays of three bubbles

This section is inspired by the earlier work of Dear and Field²⁸, which initiated the numerical studies mentioned in section III B. The present study is only interested in triangular arrangements of three bubbles of the same size, placed in the direction perpendicular to a planar incident shock wave. Based on the conclusions of sections IV B and IV A, the choice to consider only the most intense collapse for each previous case is studied and merged into several triangular arrangements (see Table II). The cases 3B-XY-1 and 3B-XY-2 show arrangements in which an incident shock wave hits the isolated bubble first (see Figure 15(a)), while the 3B-XY-1i and 3B-XY-2i arrangements propose the opposite (see Figure 15(b)).

Finally, two simulations are proposed without taking into account any symmetry, in order to evaluate this hypothesis. To summarise, in this part three influences are studied: the distance between the bubbles, the order in which the incident wave hits the bubbles, and the symmetry
<figure><img src="image_19.png" /><figcaption>FIG. 15: 2D-Schematic view of the triangular arrangements of three air bubbles. The dashed lines are part of the domain ignored in the computation due to symmetry conditions.</figcaption></figure>

conditions.

Figure 16 only shows snapshots for cases 3B-XY-1 and 3B-XY-1i, because these two cases highlight the main difference between the two types of triangular arrangements, although the differences remain visible in Figure 17.

As explained in Section IV A, for horizontal arrangements, the first bubble collapse for cases 3B-XY-1 and 3B-XY-2 is similar to the collapse of an isolated bubble (see Figure 15(a), left). For these cases, the water-hammer shockwave and the interaction waves pressure, reached during the collapse of the first bubble, have intensities of 4.61 GPa and about 14.90 GPa, respectively (see Figure 17). Then, the planar incident shock wave hits the second bubbles first, followed by the water-hammer shockwave generated by the collapse of the first bubble, implying the presence of two rarefaction waves (see Figure 16(a), left). On the one hand, the two bubbles are deformed by the joint action of two shock waves and thus experience a completely asymmetric deformation. On the other hand, the liquid jets that take place inside the bubbles do not follow a main direction in the plane orthogonal to the incident shock wave and have an intensity of about 3300 m/s at impact, which is about 20% higher than for a single bubble. At about time $t = 7.90 \ \mu s$ (see Figure 16(b), left), the second bubbles release their water-hammer shockwave and a secondary pressure peak occurs at their sides, which can be observed with the double peak in Figure 17, at about time $t = 7.90 \ \mu s$ for the case 3B-XY-1 and time $t = 8.20 \ \mu s$ for the case 3B-XY-2. Similarly to the horizontal array, the water-hammer shockwave pressure is a bit stronger for the second bubble
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0246108

<figure><img src="image_20.png" /><figcaption>FIG. 16: Shock-induced collapse of two triangular arrays of three bubbles in water near a wall. Dimensionless longitudinal velocity component $u/u_{sh}$, Schlieren-like representation colored by the dimensionless pressure ($P/P_{sh}$), dimensionless wall pressure ($P/P_{sh}$) and iso-surface of void ratio $\alpha = 0.10$. Green: Negative pressure.</figcaption></figure>
than for the first one (about 20%), and it is explained by the rarefaction waves emitted during the collapse of the second bubbles. A few moments later (about time $t = 8.60 \ \mu s$), the higher pressure peak occurs upstream with more or less the same intensity as for the first bubble.

<figure><img src="image_21.png" /><figcaption>FIG. 17: Time history of the maximum dimensionless pressure ($P/P_{sh}$) of a shock-induced collapse of a triangular arrays of three bubbles in water near a wall.</figcaption></figure>

The shock waves are then reflected off the wall and a series of recollapses of bubble fragments occur, resulting in high intensity pressure peaks in water (26.30 and 29.22 GPa for cases 3B-XY-1 and 3B-XY-2, respectively). In the context of this triangular arrangement, the maximum wall pressure peak (24.41 and 10.54 GPa for cases 3B-XY-1 and 3B-XY-2, respectively) occurs when the shock waves emitted during the collapse of each bubble overlap. There is a tiny amplification of the intensity of the peak pressure on the wall compared to the collapse of a single bubble, in the most compact case and no negative pressure appears in the low-pressure region. When the planar incident shock wave reaches the two bubbles first (cases 3B-XY-1i and 3B-XY-2i, Figure 15(b)), at about time $t = 3.80 \ \mu s$, the water-hammer shockwave is emitted with the same intensity for each bubble (4.61 GPa). The second pressure peak (time $t = 4.80 \ \mu s$), in front of the bubbles on their axis of symmetry, is less violent (12.40 GPa) in the case where the two bubbles are vertically close to each other (case 3B-XY-1i). This phenomenon has already been observed for the vertical arrangements: The merging of the rarefaction waves creates an area of low-pressure, which slows down the collapse of the bubbles and reduces their intensity. For the case 3B-XY-2i, when $\Gamma_y$ is greater, the second pressure peak is close to the case of an isolated bubble. Thus, similar to what is explained in section III B, a region of negative pressure is observed ahead of the
last bubble (see Figure 16,(a), right). Such a region is not observed in the context of a triangular arrangement where the isolated bubble is placed in front of the other two due to the less compact nature of this arrangement. For the collapse of the last bubble, in the case where the bubbles are closest together in the x-direction (case 3B-XY-1i), the intensity of the water-hammer shockwave pressure observed at time $t = 7.80$ µs is higher than for the first bubbles (about 20%), while the second pressure peak (time $t = 8.60$ µs) is less intense than for the first bubbles (about 30%) and the intensity of the peak pressure on the wall is only a half that of the single bubble collapse case. For the case 3B-XY-2i, the collapse of the last bubble is quite different. The intensity of the water-hammer shockwave pressure (time $t = 7.90$ µs) is similar to that of the first bubbles, while the second pressure peak (time $t = 8.55$ µs) is less intense (about 15%) and once again the peak pressure on the wall is much less intense than for the single bubble collapse.

Finally, one can see in figure (16(c)) that there is a strong asymmetry taking place during the collapse of the bubbles and for each of them. Then the question of the relevance of the symmetry hypothesis in the simulations is raised. In order to verify this assumption, two simulations (2-XY-1 and 2-BX-1i) have been performed without taking any symmetry into account. The grid of these calculations is about 6 billion cells. Comparisons of the time evolution of the maximum pressure in water and on the wall show no significant discrepancy (see Figure 18), and this again confirms that the use of symmetry conditions can be safely kept. This hypothesis drastically reduces the cost of the simulation while preserving the quality of the results, which is in agreement with results of Betney *et al.*²⁵, who did not observe any three-dimensional instabilities.

<figure><img src="image_22.png" /><figcaption>FIG. 18: Time history of the maximum dimensionless pressure ($P/P_{sh}$) of the shock-induced collapse of triangular arrays of three bubbles in water near awall.</figcaption></figure>
In conclusion, the shock-induced collapse of triangular arrangements, can be summarised as
follows:

* A region of negative pressure occurs only if the rarefaction waves of the two superimposed bubbles are close enough, and if these two bubbles are the first to be hit by the incident shock wave.

* There is no increase in the intensity of the collapse, but a chain reaction such as in section III B and according to previous observations²⁶.

* No significant amplification of the peak pressure is observed compared to the collapse of a single bubble, although the intensity of the water-hammer shockwave pressure for the second bubble for the smallest $\Gamma_x$ is significantly higher than for a single bubble.

* Despite significant asymmetric deformations during collapse, no three-dimensional instabilities are identified.

D. Pyramidal array of five bubbles

Given the lack of amplification of the peak pressure on the wall observed in the triangular arrangements compared to the collapse of a single bubble and the fact that the horizontal arrangement results in much higher pressure peaks, it can be concluded that the amplification of the peak pressure remains quite modest in the case of collapse in a plane arrangement. Bempedelis and Ventikos²⁶ proposed a pyramidal arrangement identical to the triangular one, with the addition of two transverse bubbles and observed an increase in the pressure peak on water. Based on this observation, simulations of cases 5B-XYZ-1 and 5B-XYZ-2 investigate such arrays and may be understood as transverse extensions of cases 3B-XY-1 and 3B-XY-2.

The incident wave first hits the isolated bubble, which collapses and emits a shock wave with a pressure intensity of 4.61 GPa (water-hammer shockwave, see Figure 19 at time $t = 5.50 \ \mu s$). This is followed by a pressure peak (interaction waves peak) behind the bubble, on its symmetry axis, with an intensity of 14.90 GPa (see Figure 20). Then, the incident wave reaches the other bubbles which begin to deform. A few moments later, the water-hammer shockwave of the first bubble reaches them (see Figure 19 at time $t = 6.35 \ \mu s$). As in the case of the triangular arrangement, this means that the fast jet forming inside the bubbles follows a different direction to that of the
<figure><img src="image_23.png" /><figcaption>FIG. 19: Shock-induced collapse of a pyramidal array of five bubbles in water near a wall. Case 5B-XYZ-1. Dimensionless longitudinal velocity component $u/u_{sh}$, Schlieren-like representation colored by the dimensionless pressure ($P/P_{sh}$), dimensionless wall pressure ($P/P_{sh}$) and iso-surface of void ratio $\alpha = 0.10$. Green: Negative pressure.</figcaption></figure>

incident wave, and the deformation of the bubbles is asymmetric. The presence of a significant area of negative pressure between the bubbles and its persistence throughout the collapse process of the pyramidal array explains why the pressure intensity of the water-hammer shockwave emitted during the collapse of the last bubbles (see Figure 19 at time $t = 8.90 \ \mu s$) is quite similar in intensity to that of the first bubble and why the intensity of the pressure peak behind the bubbles is lower compared to the single bubble collapse (about 30% and 20% for cases 5B-XYZ-1 and 5B-XYZ-2, respectively).

Figure 19 shows that pyramidal configurations greatly amplify the intensity of shock waves emitted during bubble collapse, and figure 20(b) illustrates this result with an amplification (up to almost 30 times) of the maximum loading pressure on the wall compared to the case of a single bubble.

After the incident shock wave has reflected off the wall and the bubbles have collapsed, the negative pressure area is still present on the symmetry axis. Finally, at about time $t = 11.20 \ \mu s$
<figure><img src="image_24.png" /><figcaption>FIG. 20: Time history of the maximum dimensionless pressure ($P/P_{sh}$) of a shock-induced collapse of a pyramidal array of five bubbles in water near a wall.</figcaption></figure>

and $t = 11.50$ µs for cases 5B-XYZ-1 and 5B-XYZ-2, respectively, the negative pressure region reaches the wall and just after a very strong pressure peak is applied on the wall. The intensity of the event is 30 times greater than the isolated bubble collapse, and this configuration is the one that releases the most energy, as already mentioned in²⁶. To summarise the key points for the shock-induced collapse of the pyramidal arrangement:

* Negative pressure is displayed in the interbubble space.

* Compared to the single bubble collapse, a strong amplification of the collapse intensity and a higher potential wall damage, are detected.

* The maximum pressure peak on the wall will always correspond to the maximum pressure peak in water.

E. Influence of the mass transfer term

As discussed in sections III B, IV A, IV C and IV D, negative pressure leading to cavitation, often occurs during bubble collapse, particularly in regions where bubble density is higher. Forehand *et al.* ⁴¹ demonstrated the likelihood of cavitation inside a water droplet during interaction with a normal shock wave. These results emphasize the need to revisit the mass transfer hypothesis and assess its impact. An additional term is introduced as a source term in the transport equation for
the void ratio (4) following the approach outlined by Goncalves and Charriere³³, Goncalves⁴². Thus, the void ratio equation (4) becomes:

$$ \frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = K \nabla \cdot \mathbf{u} + \frac{1}{\rho_l} \dot{m}, \quad \text{where} \quad \frac{1}{\rho_l} = \frac{\frac{c_g^2}{\alpha} + \frac{c_l^2}{1-\alpha}}{\frac{\rho_l c_l^2}{1-\alpha} + \frac{\rho_g c_g^2}{\alpha}} \quad (6) $$

where $\rho_l$ the interfacial density, $\dot{m}$ is the mass transfer between phases. Goncalves⁴² suggested that the mass transfer could be considered as proportional to the velocity divergence and proposed the following formulation:

$$ \dot{m} = \frac{\rho_l \rho_g}{\rho_l - \rho_g} \left( 1 - \frac{c^2}{c_{wallis}^2} \right) \nabla \cdot \mathbf{u} \quad (7) $$

The liquid density $\rho_l$ is assumed to be in its equilibrium state, e. g. : $\rho_l = \rho_l^{sat}(T_{ref})$, and the gas (vapor) density $\rho_g$ follows the stiffened gas EOS and varies with the temperature. To determine the effect on each type of arrangement, cases 2B-X-2, 3B-XY-1i and 5B-XYZ-2 were repeated with mass transfer taken into account.

Figure (21) illustrates the evolution of maximum pressure in the fluid and at the wall for each case. Since the primary objective of this study is to assess the amplification of wall pressure resulting from the collapse of a bubble array compared to a single isolated bubble, only the evolution of the maximum pressure is considered to evaluate the impact of mass transfer. It is worth noting that different areas of focus could lead to alternative conclusions.

Overall, the differences between cases with and without mass transfer are negligible (Figure 21). However, the conclusions may vary slightly depending on the considered arrangement.

For the horizontal arrangement cases (Figure 21, case 2B-X-2) and where the isolated bubble is positioned in front of four others (Figure 21, case 5B-XYZ-2), no significant differences are observed. A slight increase in wall pressure is noted when mass transfer is excluded, although this difference is minimal, less than 0.5%.

However, when the isolated bubble is placed downstream, just in front of the wall (Figure 21, case 3B-XY-1i), the differences in pressure peak intensity become more pronounced. Specifically, notable discrepancies emerge in the evolution of the maximum pressure, both in the fluid and at the wall. At time $t = 9.60$ µs, when the reflected wave strikes the remnants of the bubble closest to the wall, a significantly higher pressure peak (around 30% greater) occurs if mass transfer is included. Subsequently, at time $t = 9.86$ µs, a second peak appears resulting from a bubble recollapse, with lower intensity when mass transfer is accounted for, in similar proportions. Ultimately, while the
<figure><img src="image_25.png" /><figcaption>FIG. 21: Time history of the maximum dimensionless pressure ($P/P_{sh}$) of a shock-induced collapse of various arrays of vapor bubbles in water near a wall, with or without mass transfer term.</figcaption></figure>

maximum pressure in the fluid and on the wall remains nearly the same with or without mass transfer, the most intense events occur at different times.

## F. Brief quantitative analysis

Therefore, only some of the studied arrangements resulted in an increase in shock wave intensity at the wall compared to the reference case with a single bubble. While interactions between bubbles can be observed when they are close but not merged, these interactions do not systematically lead to intensify pressure loads on the wall.

To facilitate the analysis of the results, it is necessary to define the bubble density in an array. Let's assume that each bubble arrangement is enclosed within a sphere. The center of this sphere is the center of the arrangement, and its radius is defined as the distance to the farthest bubble from the center.

For simplicity, arrangements with two or three bubbles can be considered in a plane, where only the area of the enclosing circles is relevant. For arrangements involving five bubbles, the volume of the sphere must be considered.

The density of bubbles in an array is calculated as the ratio of the area (for two- and three-
bubble cases) or volume (for five-bubble cases) of the enclosing circle or sphere to the area or volume of the individual bubbles, multiplied by the number of bubbles in the array.

Thus, for the arrangements 2B-X-1 (or 2B-Y-1), 2B-X-2 (or 2B-Y-2), and 2B-X-3 (or 2B-Y-3), the densities are 4/5, 2/3, and 4/74, respectively. For cases 3B-XY-1, 3B-XY-2, 3B-XY-1i, and 3B-XY-2i, the densities are 12/25 and 6/15, respectively. Finally, for cases 5B-XYZ-1 and 5B-XYZ-2, the densities are 4/15 and 5/27, respectively.

For simplicity, only the cases where the maximum pressure at the wall is higher than that of a single bubble are considered. In all other cases, the pressure remains comparable to the reference case.

Table 22 provides a comparison of the maximum pressure peaks observed during single bubble collapse near a wall (reference case) with those in other configurations, namely bubble arrays, where the pressure peaks surpass the reference case.

<table><thead><tr><th></th><th>Maximum collapse<br>pressure peak (GPa)</th><th>Maximum wall<br>pressure peak (GPa)</th></tr></thead><tbody><tr><td>1B (Reference case)</td><td>28.08</td><td>17.65</td></tr><tr><td>2B-X-1</td><td>82.70</td><td>48.72</td></tr><tr><td>2B-X-2</td><td>62.23</td><td>62.75</td></tr><tr><td>2B-X-3</td><td>22.00</td><td>22.50</td></tr><tr><td>5B-XYZ-1</td><td>130.00</td><td>130.00</td></tr><tr><td>5B-XYZ-2</td><td>477.00</td><td>477.00</td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>FIG. 22: Maximum pressure generated during bubble collapse near a wall for various configurations, with proposed law describing the pressure peak intensity on the wall.</figcaption></figure>

The preliminary study on mesh influence highlights the importance of fine resolution around each bubble, especially for accurately capturing the cavitation phenomenon, which consistently occurs when the arrangement amplifies shock wave intensity. Despite the limitations of the proposed parametric study, a power law seems to be applicable:

$$ \frac{P_{wall,max}}{P_{sh}} = \left( \frac{p_{\infty}}{P_{sh}} \right)^k d^{-\frac{k}{2}} \qquad (8) $$

where *d* is the density ratio and *k* is equal to 4.4.
Arrangements of two or three bubbles provide valuable cases for understanding the dynamics of multiple bubble collapses near a wall. Notably, wall damage is consistently greater compared to the single-bubble scenario, with maximum wall pressure potentially reaching up to 500 times that of the incident wave that triggered the collapse.

## V. FUTURE INVESTIGATIONS

To further advance this parametric study, several aspects could be explored in the future. Given the high intensity of the incident shock wave in this context, it would be beneficial to reduce it to focus on medical issues. This adjustment will impact the cost of the parametric study, as the slower propagation speed of the incident wave will require longer simulation times.

To study wall damage, a comprehensive analysis must be performed using a fluid-structure interaction (FSI) simulation with strong coupling. Initially, the cost of this type of simulation will necessitate the use of 2D simulations.

Although the collapse of a bubble cloud has been studied but not near a wall, and considering mesh influence studies here, it seems that using our in-house code and current computational resources for such analysis is unfortunately not feasible in the near future. Nevertheless, given the current results, understanding the complex interactions between collapses in a bubble cloud near a wall is essential. An approach based on modeling with empirical laws, like the one proposed here, can help to predict the maximum intensity of pressure peaks under specific limiting assumptions. However, only comprehensive simulations can provide complete insights.

## ACKNOWLEDGMENTS

This research was supported by the French National Research Agency ANR (project 18-CE46-009). Computations have been performed on the HPC resources of GENCI under allocations A0112A10981, A0132A10981 and during a Grand Challenge project on the GPU-extension of the Jean Zay supercomputer.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request. Some videos are available. Please email the authors for access.
## 709 AUTHOR CONTRIBUTIONS

All of the authors were involved in the preparation of the manuscript and have read and approved the final manuscript version.

## 712 REFERENCES

1. S. S. Cook, "Erosion by water-hammer," Proc. R. Soc. Lond. A 119 783 (1928).

2. T. Kodama and Y. Tomita, "Cavitation bubble behavior and bubble-shock wave interaction near a gelatin surface as a study of in vivo bubble dynamics," Appl. Physics. B 70 1 (2000).

3. Lord Rayleigh, "On the pressure developed in a liquid during the collapse of a spherical cavity," Phil. Mag. 34 (1917).

4. M. S. Plesset, "The dynamics of cavitation bubbles," J. Appl. Mech. 16 (1949).

5. F. R. Gilmore, "The growth or collapse of a spherical bubble in a viscous compressible liquid," California Institute of Technology, Technical Report No 26-4 (1953).

6. J. Keller and M. Miksis, "Bubble oscillations of large-amplitude," J. Acoust. Soc. Am 68 (80).

7. M. Kornfeld and L. Suvorov, "On the destructive action of cavitation," J. Appl. Phys. 15 (1944).

8. J. K. Walters and J. F. Davidson, "The initial motion of a gas bubble formed in an inviscid liquid. Part 1. The two-dimensional bubble." J. Fluid Mech. 12 (1962).

9. J. K. Walters and J. F. Davidson, "The initial motion of a gas bubble formed in an inviscid liquid. Part 2. The three-dimensional bubble and the toroidal bubble." J. Fluid Mech. 17 (1963).

10. T. B. Benjamin and A. T. Ellis, "The collapse of cavitation bubbles and the pressures thereby produced against solid boundaries," Phil. Trans. R. Soc. Lond. A 260 (1966).

11. A. Shima, Y. Tomita, and K. Takahashi, "The collapse of a gas bubble near a solid wall by a shock wave and the induced impulsive pressure," Institution of Mechanical Engineers, Proceeding 198 (1984).

12. Y. Tomita and A. Shima, "Mechanisms of impulsive pressure generation and damage pit formation by bubble collapse," J. Fluid Mech. 169 (1986).

13. N. K. Bourne and J. E. Field, "Shock-induced collapse of single cavities in liquids," J. Fluid Mech. 244, 225-240 (1992).

14. G. J. Ball, B. P. Howell, T. G. Leighton, and M. J. Schofield, "Shock-induced collapse of a cylindrical air cavity in water: a free-Lagrange simulation," Shock Waves 10 (2000).
¹⁵R. R. Nourgaliev, T. N. Dinh, and T. G. Theofanous, "Adaptive characteristics-based matching for compressible multifluid dynamics," J. Comput. Phys. **213** (2006).

738

739

¹⁶E. Johnsen and T. Colonius, "Numerical simulations of non-spherical bubble collapse," J. Fluid Mech. **629** (2009).

740

741

¹⁷V. Coralic and T. Colonius, "Finite-volume weno scheme for viscous compressible multicomponent flows," J. Comput. Phys. **274** (2014).

742

743

¹⁸N. Hawker and Y. Ventikos, "Interaction of a strong shockwave with a gas bubble in a liquid medium: A numerical study," J. Fluid Mech. **701**, 59-97 (2012).

744

745

¹⁹S. A. Beig and E. Johnsen, "Maintaining interface equilibrium conditions in compressible multiphase flows using interface capturing," J. Comput. Phys. **302** (2015).

746

747

²⁰Y.-L. Yoo and H.-G. Sung, "Modeling for non isothermal cavitation using 4 equation models," Int. J. Heat Mass Transf. **127** (2018).

748

749

²¹N. Bempedelis and Y. Ventikos, "A simplified approach for simulations of multidimensional compressible multicomponent flows: the grid-aligned ghost fluid method," J. Comput. Phys. **405** (2020).

750

751

752

²²E. Goncalves and P. Parnaureau, "Comparison of multiphase models for computing shock-induced bubble collapse," International Journal of Numerical Methods for Heat and Fluid Flow **22**, 3845-3877 (2020).

753

754

755

²³E. Lauer, X. Y. Hu, S. Hickel, and N. A. Adams, "Numerical modelling and investigation of symmetric and asymmetric cavitation bubble dynamics," Comput. & Fluids **69** (2012).

756

757

²⁴E. Lauer, X. Y. Hu, S. Hickel, and N. A. Adams, "Numerical investigation of collapsing cavity arrays," Phys. Fluids **24** (2012).

758

759

²⁵M. R. Betney, B. Tully, N. A. Hawker, and Y. Ventikos, "Computational modelling of the interaction of shock waves with multiple gas-filled bubbles in a liquid," Phys. Fluids **27** (2015).

760

761

²⁶N. Bempedelis and Y. Ventikos, "Energy focusing in shock-collapsed bubble arrays," J. Fluid Mech. **900** (2020).

762

763

²⁷A. Shima, Y. Tomita, and T. Ohno, "Collapse of multiple gas bubbles by a shock wave and induced impulsive pressure," J. Appl. Phys. **56** (1984).

764

765

²⁸J. P. Dear and J. E. Field, "A study of the collapse of arrays of cavities," J. Fluid Mech. **190**, 409-425 (1988).

766

767

²⁹A. Tiwari, C. Pantano, and B. Freund, "Growth-and-collapse dynamics of small bubble clusters near a wall," J. Fluid Mech. **775**, 1-23 (2015).

768

769
³⁰F. Wermelinger, U. Rasthofer, P. Hadjidoukas, and P. Koumoutsakos, “Petascale simulations of compressible flows with interfaces,” J. Computational Science **26** (2018).

³¹U. Rasthofer, F. Wermelinger, P. Karnakov, J. Sukys, and P. Koumoutsakos, “Computational study of the collapse of a cloud with 12500 gas bubbles in a liquid,” Phys. Rev. Fluids **4** (2019).

³²E. Goncalves and P. Parnaudeau, “Numerical study of pressure loads generated by a shock-induced bubble collapse,” Phys. Fluids **33** (2021).

³³E. Goncalves and B. Charriere, “Modeling for isothermal cavitation with a four-equation model,” International Journal of Multiphase Flow **59**, 54–72 (2014).

³⁴R. Dubois, E. Goncalves, and P. Parnaudeau, “High performance computing of stiff bubble collapse on cpu-gpu heterogeneous platform,” Comput. Math. Appl. **99**, 246–256 (2021).

³⁵E. Toro, M. Spruce, and W. Speares, “Restoration of the contact surface in the hll-riemann solver,” Shock Waves **4**, 25–34 (1994).

³⁶B. van Leer, “On the relation between the upwind-differencing schemes of godunov,” Journal on Scientific and Statistical Computing **5**, 1–20 (1984).

³⁷J. H. J. Niederhaus, J. A. Greenough, J. G. Oakley, D. Ranjan, M. H. Anderson, and R. Bonazza, “A computational parameter study for the three-dimensional shock–bubble interaction,” J. Fluid Mech. **594**, 85–124 (2008).

³⁸A. R. Jamaluddin, J. Ball, C. K. Turangan, and T. G. Leighton, “The collapse of single bubble and approximation of the far-field acoustic emissions for cavitation induced by shock wave lithotripsy,” J. Fluid Mech., 1–37 (2001).

³⁹C. K. Turangan, J. Ball, A. R. Jamaluddin, and T. G. Leighton, “Numerical studies of cavitation erosion on an elastic–plastic material caused by shock-induced bubble collapse,” Proc. R. Soc. A **473** (2017).

⁴⁰N. Apazidis, “Numerical investigation of shock induced bubble collapse in water,” Phys. Fluids **28** (2016).

⁴¹R. W. Forehand, K. C. Nguyen, C. J. Anderson, R. Shannon, S. M. Grace, and M. P. Kinzel, “A numerical assessment of shock–droplet interaction modeling including cavitation,” Phys. Fluids **35** (2023).

⁴²E. Goncalves, “Numerical study of expansion tube problems: Toward the simulation of cavitation,” Computers & Fluids **72**, 1–19 (2013).