## Abstract

Material deformation and pitting from cavitation bubble collapse is investigated using fluid and material dynamics and their interaction. In the fluid, a novel hybrid approach, which links a boundary element method and a compressible finite difference method, is used to capture non-spherical bubble dynamics and resulting liquid pressures efficiently and accurately. The bubble dynamics is intimately coupled with a finite-element structure model to enable fluid/structure interaction simulations. Bubble collapse loads the material with high impulsive pressures, which result from shock waves and bubble re-entrant jet direct impact on the material surface. The shock wave loading can be from the re-entrant jet impact on the opposite side of the bubble, the fast primary collapse of the bubble, and/or the collapse of the remaining bubble ring. This produces high stress waves, which propagate inside the material, cause deformation, and eventually failure. A permanent deformation or pit is formed when the local equivalent stresses exceed the material yield stress. The pressure loading depends on bubble dynamics parameters such as the size of the bubble at its maximum volume, the bubble standoff distance from the material wall and the pressure driving the bubble collapse. The effects of standoff and material type on the pressure loading and resulting pit formation are highlighted and the effects of bubble interaction on pressure loading and material deformation are preliminarily discussed.

## 1. Introduction

Cavitation is the explosive growth and intense collapse of bubble nuclei in a liquid when exposed to large pressure variations. It has been long known by engineers for its deleterious effects in a variety of applications involving fluid machinery, pumps, propellers, valves, sluice gates, lifting surfaces, etc. Unwanted effects of cavitation include material erosion, performance degradation, and noise and vibrations.

Over the past century, since the early works of Besant [1] and Rayleigh [2], a significant amount of studies have been dedicated to seeking out measures to mitigate, if not eliminate, cavitation. Cavitation erosion has attracted wide attention from the materials and fluids engineering, and scientific communities not only because significant damage to the machinery is frequently observed, but also because the underlying physics of the erosion process itself involves interesting and complex fluid–structure interaction dynamics, which need to be addressed.

Conversely, cavitation can be generated for useful purposes and many studies have been conducted to explore cavitation merits and enhance its intensity for efficient underwater cleaning; e.g. [3–6], ultrasonic cleaning [7–9], chemical compounds oxidation, e.g. [10,11], microorganisms disinfection, e.g. [12], algae oil extraction, e.g. [13,14], etc.

Recent developments in utilization of ultrasound cavitation for biomedical applications reveal a need for delicate control of the cavitation intensity, because the boundaries between meritorious effects (e.g. tumour removal) and deleterious effects (e.g. kill of nearby healthy cells and tissue) depend on some critical irradiation set-up conditions. For example, use of cavitation to produce cell sonoporation produces therapeutic effects only when the detrimental ‘side effects' on cells due to over dosage does not occur [15–17]. Numerous studies have also observed that cavitation activity could result in deleterious bio-effects such as haemolysis [16,18], haemorrhage [19–21], kidney failure [22] when patients are treated with ultrasound-based medical treatments such as shock wave lithotripsy or high-intensity focused ultrasound.

Even though different mechanisms can cause cavitation and different cavitation types exist, the underlying physics of materials erosion and of cell membranes poration is the same from the microscopic bubble dynamics point of view. It has been known for more than a century [1,2] for cavitation and underwater explosion bubbles that volume implosion of the bubble can generate very high-pressure pulses and shock waves. Many pioneering studies have also shown, experimentally as well as analytically, that the collapse of these bubbles near rigid boundaries results in high-speed re-entrant liquid jets, which penetrate the highly deformed bubbles to strike a nearby rigid boundary generating water hammer like impact pressures [23,24]. The resulting re-entrant jet has been found a key element in pit formation on the material surface [25–27] as well as to induce the erosion of endothelia [28] and haemorrhage [29].

Re-entrant jets can also be induced through bubble/bubble interaction. For example, re-entrant jets induced by tandem bubble interactions have recently been introduced and gained attention for its precise control on cell membrane poration [30,31]. Also, multi-bubble interactions in a bubble cloud can result in collective synergetic behaviour, which can result in very high pressures and loading on nearby structures [6,32–36].

Cavitation bubble collapse near boundaries has been extensively investigated numerically by assuming flow incompressibility and using the boundary element method (BEM) [37–41]. The BEM enables to accurately describe the re-entrant jet formation, development, advance through the bubble and impact on nearby boundaries. It can provide the jet geometric and dynamic characteristics as functions of time and space as the bubble wall velocities (including the re-entrant jet velocity) are most often small relative to the speed of sound in water. On the other hand, during the bubble explosive growth, rebound, and at jet impact, compressible effects can be non-negligible and need to be included. These dynamic stages may lead to shock wave formation and shock propagation and impact on a nearby material. Compressible flow solvers have been developed for this purpose [42–44]. However, such flow solvers often use a finite difference method which requires very fine spatial resolution and small time-step sizes especially to resolve the re-entrant jet accurately [45,46]. This makes them not very efficient for simulating the relatively long duration bubble period but excellent at describing the shock phase.

To overcome this, this work applies a novel hybrid numerical procedure involving a time domain decomposition into incompressible and compressible stages in order to capture the full bubble dynamics period as well as the shock phase occurring during re-entrant jet impact, bubble collapse and bubble rebound. This numerical procedure takes advantage of an accurate shock-capturing method and of a BEM shown to be both very efficient in modelling cavitation bubble dynamics problems and very accurate in capturing the re-entrant jet and its dynamics.

In proximity of a deforming boundary, the complex flow phenomena and the bubble dynamics itself can be significantly altered [27,47–51]. In this study, a finite-element structure code is coupled with the hybrid fluid flow solver to investigate material response to the pressure field generated by the bubble dynamics.

This paper presents first the methods used, then describes the pressure loading and material deformation resulting from a single bubble collapse near a material. The effect of bubble/bubble interaction on pressure loading and on material deformation resulting for tandem bubbles and bubble clouds are then considered. Although the actual case studies shown are for cavitation erosion on a flat material surface, the key findings and the numerical approaches apply as well to complex geometries (with appropriate corrections) and to biomedical applications such as for the study of sonoporation and breakage and removal of calculi.

## 2. Numerical approach

The numerical approach applied to model material pitting in this paper is part of a general hybrid approach which was developed by the authors to simulate fluid–structure interaction (FSI) problems involving shock and bubble pressure pulses [52,53]. As illustrated in figure 1, for a highly inertial bubble such as a spark-generated bubble, an underwater explosion bubble (Undex), or a laser-generated bubble [54–57], a compressible–incompressible link is required at the beginning to handle the emitted shock wave and the flow field generated by the exploding bubble. Cavitation bubbles on the other hand, generate a small pressure peak and no shock wave during the growth phase. As a result, no initial shock phase compressible solution is required. An incompressible BEM code can then be used to simulate most of the bubble period until the end of the bubble collapse where, due to high liquid speeds or to the bubble re-entrant jet impacting on the liquid or on the structure, compressible flow effects prevail again.

In this study, the solution of the incompressible BEM code is passed to a fully compressible code capable of shock capturing to simulate re-entrant jet impact and bubble ring collapse.

### 2.1. Incompressible flow modelling

The potential flow code used in this study, 3DynaFS-Bem^{©}, is based on a BEM [58–60]. The code solves the Laplace equation for the velocity potential, *ϕ*,
2.1with the velocity vector defined as . A boundary integral method is used to solve the Laplace equation based on Green's theorem
2.2

In this expression, *Ω* is the domain of integration having elementary volume d*Ω*. The boundary surface of *Ω* is *S*, which includes the surfaces of the bubble and the nearby boundaries with elementary surface element d*S*. **n** is the local normal unit vector. is Green's function, where **x** corresponds to a fixed point in Ω and **y** is a point on the boundary surface *S*. Equation (2.1) reduces to Green's formula with being the solid angle under which x sees the domain, *Ω*:
2.3where is the solid angle. To solve (2.3) numerically, the BEM, which discretizes the surface of all objects in the computational domain into panel elements, is applied.

Equation (2.3) provides a relationship between *ϕ* and ∂*ϕ*/∂*n* at the boundary surface *S*. Thus, if either of these two variables (e.g. *ϕ*) is known everywhere on the surface, the other variable (e.g. ∂*ϕ*/∂*n*) can be obtained.

To advance the solution in time, the coordinates of the bubble and any free surface nodes, **x**, are advanced according to
2.4

The velocity potential *ϕ* on the bubble and free surface nodes is obtained through the time integration of the material derivative of *ϕ*, i.e. d*ϕ*/d*t*, which can be written as
2.5where ∂*ϕ*/∂*t* can be determined from the Bernoulli equation
2.6 is the hydrostatic pressure at infinity at *z* *=* *0*, where *z* is the vertical coordinate. *p*_{l} is the liquid pressure at the bubble surface, which balances the internal pressure and surface tension,
2.7where *p*_{v} is the vapour pressure, *σ* is the surface tension, and is the local bubble wall curvature. *p*_{g} is the gas pressure inside bubble and is assumed to follow a polytropic law with a compression constant, *k*, which relates the gas pressure to the gas volume, , and reference (often initial) value, and .
2.8

### 2.2. Compressible flow modelling

The multi-material compressible Euler equation solver used here is based on a finite difference method [46,61]. The code solves continuity and momentum equations for a compressible inviscid liquid in Cartesian coordinates. These can be written in the following format:
2.9and
where *ρ* is the fluid density, *p* is the pressure, *u*, *v* and *w* are the velocity components in the *x*, *y*, *z* directions, respectively (*z* is vertical), *g* is the acceleration of gravity, and *e*_{t} *=* *e* *+* 0.5(*u ^{2}*

*+*

*v*

^{2}*+*

*w*) is the total energy with

^{2}*e*being the internal energy. The system is closed by using the equation of state of each material, which provides for this material the pressure as a function of the specific internal energy and the density. In this study, a

*γ*-law (with

*γ*= 1.4) is used for the gas–vapour mixture [62] 2.11and the Tillotson equation is used for water [63] where

*ω*,

*A*,

*B*and

*C*are constants, and

*p*

_{0},

*e*

_{0}and

*ρ*

_{0}are the reference pressure, reference specific internal energy and reference density, respectively:

*ω*= 0.28,

*A*= 2.20 × 10

^{9}Pa,

*B*= 9.54 × 10

^{9}Pa,

*e*

_{0}= 3.54 × 10

^{5}m

^{2}s

^{−2},

*p*

_{0}= 1.0 × 10

^{5}Pa,

*ρ*

_{0}= 1000 kg m

^{−3}.

The compressible flow solver uses a high-order Godunov scheme and employs the Riemann problem to construct a local flow solution between adjacent cells. The numerical method is based on a higher-order MUSCL scheme and tracks each material. To improve efficiency, an approximate Riemann problem solution replaces the full problem. The MUSCL scheme is augmented with a mixed cell approach [64] to handle shock wave interactions with fluid or material interfaces. This approach uses a Lagrangian treatment for the cells including an interface and an Eulerian treatment for cells away from interfaces. A re-map procedure is employed to map the Lagrangian solution back to the Eulerian grid [46,61] The code has been extensively validated against experiments [45,49].

### 2.3. Compressible–incompressible link procedure

Both incompressible and compressible flow solvers are able to model the full bubble dynamics on their own. However, each method has its shortcomings when it comes to specific parts of the bubble history. The BEM-based incompressible flow solver is efficient, reduces the dimension of the problem by one (line integrals for an axisymmetric problem, and surface integrals for a three-dimensional problem) and thus allows very fine gridding and increased accuracy with reasonable computations times. It has been shown to provide re-entrant jet parameters and speed accurately [39,58,60,65]. However, it has difficulty pursuing the computations beyond surface impacts (liquid–liquid and liquid–solid).

On the other hand, the compressible flow solvers used here (Gemini developed by the US Navy at NSWCIH [46] and 3DynaFS-Comp [61]) are most adequate to model shock wave emission and propagation, liquid–liquid, and liquid–solid impacts. The method requires, however, very fine grids and very small time steps to resolve shock wavefronts. This makes it appropriate to model time portions of the bubble dynamics. Concerning bubble–liquid interface and the re-entrant jet dynamics, the procedure is diffusive as the interface is not directly modelled, and re-entrant jet characteristics are usually less accurate than obtained with the BEM approach.

Hence, we employ here a general novel approach, which combines the advantages of both methods and consists in executing the following steps:

(1) Set-up the initial flow field using the Eulerian compressible flow solver and run the simulation until the initial shock fronts go by and the remnant flow field can be assumed to be incompressible with bubble dynamics independent of initial treatments.

(2) Transfer at that instant to the Lagrangian BEM potential flow solver, all the flow field variables needed for the solution to proceed: geometry, bubble pressure and moving boundary normal velocities, .

(3) Solve for bubble growth and collapse using fine BEM grids to obtain a good description of the re-entrant jet until the point where the jet is very close to the opposite side of the bubble.

(4) Transfer the solution back to the compressible flow solver with the required flow variables. To do so, compute using the Green equation all flow field quantities on the Eulerian grid.

(5) Continue advancing the solution with the compressible code to obtain pressures due to jet impact and to the collapse of the remnant bubble ring.

### 2.4. Structure dynamics modelling

To model the dynamics of the material, the finite-element model Dyna3D is used. Dyna3D is a nonlinear explicit structure dynamics code developed by the Laurence Livermore National Laboratory [66]. Here, it computes the material deformation when the loading is provided by the fluid solution. Dyna3D uses a lumped mass formulation for efficiency. This produces a diagonal mass matrix **M**, to express the momentum equation as
2.13where **F**_{ext} represents the applied external forces, and **F**_{int} the internal forces. The acceleration, for each element, is obtained through an explicit temporal central difference method. Additional details on the general formulation can be found in [66].

### 2.5. Material models

In this study, two metal alloys: Aluminium Al7075 and Stainless Steel A2205, two rubbers: a Neoprene synthetic rubber and a Polyurethane, and a Versalink-based polyurea are considered. The metals are modelled using elastic-plastic models with linear slopes (moduli); one for the initial elastic regime and the second, a tangent modulus, for the plastic regime. The material parameters required for the current elastic-plastic model are specified as follows:

— Aluminium Al7075: Young modulus 71.7 GPa, tangent modulus 680 MPa, yield stress 503 MPa and density 2.81 g cm

^{−3}.— Stainless Steel A2205: Young modulus 190 GPa, tangent modulus 705 MPa, yield stress 515 MPa and density 7.88 g cm

^{−3}.

Concerning the rubbers, a Blatz-Ko's hyper-elastic model is used. This model is appropriate for materials undergoing moderately large strains and is based on the implementation in [67]. The material motion satisfies the formal equation
2.14where *ρ _{s}* is the material density,

**x**is the position vector of the material and

**is the Cauchy stress tensor.**

*σ*If we denote **F** the deformation gradient tensor, then the Cauchy stress tensor is computed by
2.15where *J* is the determinant of **F**, and **S** is the second Piola–Kirchoff stress tensor and is computed by
2.16

In the above equation, *G* is the shear modules, *V* is the relative volume to original which is related to excess compression and density, *C _{ij}* is the right Cauchy–Green

*strain*tensor, and

*ν*is Poisson's ratio and is set to 0.463. In this study, the shear module and density are specified as 99.7 MPa and 1.18 g cm

^{−3}for the Neoprene synthetic rubber (Rubber no. 1) and 34.2 MPa and 1.20 g cm

^{−3}for Polyurethane (Rubber no. 2).

The Versalink-based polyurea was modelled as a viscoelastic material with the following time-dependent values for the shear modulus, *G*:
2.17with
2.18

While the time dependence of *G* is considered in this model, the bulk modulus, *K*, was assumed to be constant. The material parameters in (2.18) are derived from simplification of fourth-order Prony series in [68] taking the second term as the major contributor.

### 2.6. Fluid–structure interaction coupling

Fluid–structure interaction effects are captured in the simulations by coupling the fluid codes (Gemini or 3DynaFS-Comp^{©}) with Dyna3D through a fluid–structure coupler interface (FSC). The coupling is achieved through the following steps:

— The fluid code solves the flow field and deduces the pressures at the structure surface using the positions and normal velocities of the wetted body nodes

**.**— In response, the structure code computes material stresses, strains and deformations and velocities of the wetted interface in response to this loading.

— The new coordinates and velocities of the structure surface nodes become the new boundary conditions for the fluid code at the next time step.

Additional details on the procedure can be found in [41,45,51,59]. This FSI coupling procedure has only a first-order time accuracy. A predictor–corrector approach is also implemented in the coupling to iterate and improve the solution but was not used here. This is because the numerical error due the time lag is negligible thanks to the very small time steps used. These are controlled by the steep pressure waves, which have a time scale that is two orders of magnitude shorter than the time response of the material. This method has been shown in Undex studies to correlate very well with experiments [69].

## 3. Single cavitation bubble collapse near wall

### 3.1. Problem description

We consider an initially spherical bubble of radius 50 µm, located at a distance of *X* = 1.5 mm from a flat material surface and subject it to a time-varying pressure field as represented in figure 2 and expressed as follows:
2.19

This imposed pressure variation is different from that used in many classical studies on bubble collapse near a wall and where a bubble with a maximum radius is suddenly subjected to a pressure higher than the internal pressure such as in [71–74]. Here, bubble growth is included (this allows one to include standoff distances smaller than the bubble maximum radius and covers a large range of applications), and the time-varying pressure field represents for example the pressure encountered by a bubble nucleus captured in the shear layer of a cavitating jet. In the considered expression, the bubble nucleus which is in equilibrium at atmospheric pressure is entrained in a jet low pressure region, travels to the wall and enters in stagnation region of the cavitating jet impacting the wall. The bubble transport in the flow is not included here. The stagnation pressure at the wall, which is a portion of the nozzle pressure can exceed 10 MPa as 10^{4} psi jets (70 MPa) are commonly used in practice. Similarly, a bubble nucleus travelling near a propeller blade can encounter a pressure, which drops to a value close to the vapour pressure and then, in the pressure recovery region, can be impacted by a shock wave emitted from the collapse of a nearby cloud collapse, which can exceed 10 MPa.

The duration of the pressure drop is selected to correspond to practical distances of travel of the bubble in the above two cases. The exact value is selected to achieve a normalized bubble standoff from the wall of 0.75 in order to compare between cases. The nucleus size is chosen arbitrarily to be 50 µm. A complete study requires varying this size over the range of bubble nuclei distribution in water. However, we expect the results to depend more on *R*_{max}. Figure 2*b* shows an approximate time history of the bubble radius in response to the imposed pressure field obtained by solving the Rayleigh–Plesset equation ignoring for now the presence of the wall. It is seen that the bubble starts an explosive grow as soon as the pressure suddenly drops to a low pressure below the critical pressure. The bubble continues to grow asymptotically until the pressure rises back to the bubble ‘collapse driving pressure’, *P*_{d}, here of 10^{7} Pa.

### 3.2. Non-spherical bubble dynamics and re-entrant jet development

To simulate the bubble dynamics near the wall, the incompressible BEM solver is first applied with a total of 400 nodes and 800 panels used to discretize the bubble surface. This corresponds to a grid density, which provides grid independent solution [41]. Figure 3 shows the variations of the bubble outer contours as time advances. As the bubble grows between *t* = 0 and *t* ∼ 2.4 ms, it behaves almost spherically on its portion away from the wall, while the side close to the material flattens and expands in the direction parallel to the wall actually never touching the wall as a layer of liquid remains between the bubble and the wall. Such a behaviour has been confirmed experimentally by Chahine *et al*. [41]. Note that at maximum bubble volume, the non-dimensional standoff, is less than one (here ), where *R*_{max} is the maximum equivalent bubble radius deduced from its volume.

The bubble will continue expanding following pressure reversal due to the inertia of the outward flow of the liquid. Due to the asymmetry of the flow, the pressures at the bubble interface on the side away from the wall are much higher than those near the material; thus, the collapse proceeds with the far side moving towards the material wall. The resulting acceleration of the liquid flow perpendicular to the bubble-free surface develops a Rayleigh–Taylor instability [75–78] at the axis of symmetry, which results into a re-entrant jet that penetrates the bubble and moves much faster than the rest of the bubble surface to impact the opposite side of the bubble and the material boundary.

In the present approach, the simulation of the bubble dynamics is switched from the BEM to the compressible flow solver right before the jet touches the opposite side of the bubble (dubbed touchdown). Ideally, the time of the ‘link’ between incompressible and compressible approaches should be at the time the bubble becomes multi-connected. However, to avoid increased errors/fluctuations in the BEM solution when the distance between a jet panel and the opposite bubble side panels continues to decrease as the jet advances, the ‘link’ time is selected to be when the distance between the jet front and the opposite bubble surface becomes less than or equal to 1.5 times the local panel size. This results in an underestimate of *t*_{link} by less than 1%.

Figure 4*a* shows the pressure contours and velocity vectors at the selected time, *t*_{link}, when the compressible flow solver starts its computations and FSI effects become non-negligible [79]. Figure 4*b* shows the corresponding velocity vectors and velocity magnitude contour levels. Note that, for this bubble collapse condition, prior to jet impact on the opposite side of the bubble, the liquid velocities near the tip of the jet exceeds 1400 m s^{−1}. The maximum liquid velocity of the jet exceeds the sound speed after the ‘link’ time when the computation was continued with the compressible code. The peak value reaches about 1600 m s^{−1} at the time when the jet touches down the opposite side of the bubble. This value is much higher than values reported in the literature for bubble collapse under atmospheric conditions, such as in ‘shallow’ underwater explosions or for laser or spark-generated bubbles where the reported values are of the order of 100 m s^{−1}. This is because the re-entrant jet velocity is proportional to the square root of the collapse driving pressure, *P*_{d}. Detailed discussions of the effect of *P*_{d} on the re-entrant jet velocity can be found in [25,26].

### 3.3. Re-entrant jet impact and bubble ring collapse

The compressible flow solver is used to continue the simulation after re-entrant jet touchdown. An axisymmetric domain with a total of 220 × 1470 grid points in a 1 m × 1 m domain was used with stretched grids concentrated in the immediate region surrounding the bubble. The grids were distributed such that there was a uniform fine mesh with a size of 10 µm in the area of interest where the interaction between bubble and plate is important as shown in figure 5. The axisymmetric computational domain ended at the large distance of 1 m (20 000 *R*_{0}) in the far field (radial direction and away from the wall) and at the plate material wall. A reflection boundary condition was imposed on the axis of symmetry, i.e. all physical variables such as density, pressure, velocities and energy are reflected from the axis, while transmission non-reflective boundary conditions (i.e. the flow variables are extrapolated along the characteristic wave direction) were imposed at the far-field boundaries.

Figure 6 shows the bubble shapes and corresponding pressure contours computed at six time instances after *t*_{link} = 2.435 ms. It is seen that at *t*−*t*_{link} = 0.05 µs, the jet has completely penetrated the bubble and touched the opposite side. The liquid–liquid impact event generates a localized high-pressure region which then expands quasi-spherically to reach the material liquid interface at *t*−*t*_{link} = 0*.*2 µs. The volume of the bubble ring remaining after the jet touchdown shrinks and reaches a minimum at *t*−*t*_{link} = 0.7 µs. The collapse of the bubble ring generates another high-pressure wave, which then propagates towards the axis of the cylindrical domain and reaches the wall at *t*−*t*_{link} = 0.9 µs.

To display better the pressure field dynamics during the bubble collapse impulsive loads period, a zoom of the pressure contours in the region between the re-entrant jet impact on the opposite side of the bubble and the rigid wall is shown in figure 7. Eight time steps are selected between the jet impact and the shock wave reflection from the wall. The re-entrant jet impact in figure 7*a* forms a strong quasi-ellipsoidal shock wave, which later propagates outwards in all directions losing intensity from reflections into the bubble-free surface but remaining strong in the axial region of the bubble (figure 7*b*,*c*). The shock front advances towards the solid wall at the fluid sound speed and impacts in between frames (*d*) and (*e*). It then reflects as a reinforced shock wave between frames (*e*) and (*f*). One can also observe the generation of a small bubble ring (dark blue spots on either side of the axis) in frames (*e*) and (*f*) due to the generated local high tensile stresses. This bubble ring continues to grow and becomes a clear two-phase interface as seen in figure 6*b*–*e*, while the main larger bubble ring continues to shrink. This is associated with the vortical flow induced in the re-entrant jet shear layer after the jet pierces the bubble.

The corresponding pressure distribution along the axis of symmetry at the last four time steps is shown in figure 8. It can be seen that the incoming pressure of maximum amplitude 700 MPa (red curve) is almost doubled following reflection at the rigid wall to about 1.3 GPa (green curve). The magnitude of this high-pressure loading level has been extracted or deduced from both numerical studies and from experimental measurements [26,80,81]. Figure 9 shows the time history of the bubble equivalent radius and the pressure at the axis of symmetry on the wall surface (*z* = 0). One can observe the generation by the bubble collapse of two distinct pressure peaks, one resulting from the re-entrant jet impact at the wall and the other occurring right after the remainder ring bubble reaches its minimum size. In addition to these two peaks, many other pressure peaks are observed due to pressure or shock waves bouncing back and forth between the target wall, the bubble surface, and any other daughter bubbles in the near wall flow field.

### 3.4. Dynamics response of material due to pressure loadings

The pressures generated during bubble collapse and rebound are at least two orders of magnitude higher than those generated during the bubble growth period when the bubble dynamics is that due to a driving pressure as in figure 2 [25]. Also, material deformation during bubble dynamics up to re-entrant jet impact has very little influence on the re-entrant jet impact and bubble ring collapse phases. Consequently, FSI simulations in this study are only carried out after the incompressible–compressible link has occurred and when the compressible code is used to simulate the flow field.

For the structure computations, a circular plate with a radius of 1 m and a thickness of 0.01 m is discretized using rectangular brick elements. As shown in figure 10, a stretched grid with 220 elements in the radial direction and 446 elements in the axial direction are used to discretize the plate. The elements are distributed such that a uniform fine 10 µm mesh size exists near the centre of the plate where the high-pressure loading occurs. The mesh sizes were tested to establish convergence and grid independence of the solution. The motion of the nodes at the plate bottom was restricted in all directions. The nodes along the vertical axis were only allowed to move in the vertical direction.

Figure 11 shows a time sequence of the contours of the von Mises equivalent stresses in the material for one of the two metallic alloys considered here, Al7075. It is seen that high stresses appear at the plate centre near the surface when the re-entrant jet impact pressure reaches the wall first at *t*−*t*_{link} = 0*.*2 µs. The high stress wave is observed to propagate and move radially away from the impact location. As the first high stress wave starts to attenuate, another high stress is observed initiating from the top centre of the plate at *t*−*t*_{link} = 0.9 µs, the time at which the high-pressure wave generated by the collapse of the remaining bubble ring reaches the wall (figure 6*f*).

All high stresses due to the bubble dynamics eventually attenuate. However, residual stresses remain in the material below the surface due to the plastic deformations of some regions of the plate. In the conditions of figure 11, these have their highest value occurring at a depth of 0.2 mm below the surface as shown in figure 11*f*. As the material is modelled as elastic-plastic, permanent deformation should occur wherever the local equivalent stresses exceed the material yield point. With the yield stress of Al7075 being 503 MPa, all regions that have seen red stress contour levels as shown in figure 11, experience permanent deformation due to loads from either re-entrant jet impact or bubble ring collapse.

To quantitatively examine the material response to the pressure loading, the time histories of the liquid pressure and the vertical displacement of the material surface at the centre of the Al7075 plate/liquid interface are shown together in figure 12. The material starts to get compressed as the high-pressure loading due to the re-entrant jet impact reaches it, and the plate surface centre point starts to move into the material direction at *t*−*t*_{link} = 0.45 µs. The maximum deformation occurs when the highest pressure loading peak due to the bubble ring collapse reaches the centre of the plate at time *t*−*t*_{link} = 1.15 µs. Once the pressure loading due to the full bubble dynamics has virtually vanished at *t*−*t*_{link} = 4 µs, the surface elevation continues to oscillate due to stress waves propagating back and forth through the metal alloy thickness and lack of damping in the model. Finally, a permanent deformation in the form of a pit remains as a result of the high-pressure loading causing local stresses that exceed the Al7075 elastic limit. The vertical displacement of the monitored location eventually converges to a non-zero value of about 9 µm.

The radial extent of the permanent deformation is shown in figure 13. The permanent deformation generated on the plate surface shows a profile, which is qualitatively similar to those observed in previous experimental studies such as in [80,82,83].

### 3.5. Effect of normalized standoff distance

It is known from previous studies on an explosion bubble near a rigid wall [41,65,79,84] that the impact pressure due to the re-entrant jet attains a maximum at non-dimensional distances between half and three quarters of the maximum bubble radius, i.e. In those studies, the jet hits the wall almost at the same time when re-entrant jet touchdown occurs even though a small liquid film always exists between bubble and wall. However, for this study, as shown in figure 6, the jet touchdown occurs at a small distance away from the wall. The shock generated by the jet touchdown needs to travel a distance and thus attenuates before reaching the wall. To show the effect of a direct jet impact, a smaller standoff at illustrated in figure 14 with the pressure contours at six selected time. One can see clearly that the re-entrant jet directly impacts the wall when the jet touchdown occurs at *t* − *t*_{link} = 0.05 µs. Similar to the case, the volume of the bubble ring remaining after the jet touchdown continues to shrink and reaches a minimum at about *t* − *t*_{link} = 0.92 µs. However, the high-pressure shock wave generated by the bubble ring collapse originates much closer to the wall as compared with the case. This results in a higher concentrated pressure loading at the plate centre when the shock wave propagates towards and reaches the axis at *t* − *t*_{link} = 1.1 µs.

Figure 15 shows the pressure versus time monitored at the plate centre for different standoff distances. It is seen that under the present conditions, the pressure loading due to the jet impact is much higher for because the jet directly impacts on the wall when it penetrates the other side of the bubble. The higher pressure loading due to the direct jet impact and later more concentrated pressure loading due to the ring collapse is expected to results in a different pit shape on material surface. As shown in figure 16, the pit radius is smaller with than with , while the pit depth is larger with than with because of the higher magnitude and concentrated pressure loadings.

### 3.6. Effect of material compliance

As seen in figures 9, 12 and 15, two main pressure peaks can be attributed to the jet impact and bubble ring collapse. As shown below, the magnitudes of these pressure peaks depend on the level of deformation of the material when the pressure loading on the material is a result of full interaction between the collapsing bubble and the responding material. Figure 17 shows a comparison of the time histories of the collapse impulsive load between a non-deformable wall (where FSI was not exercised and a rigid boundary condition was used) and an Al7075 deformable plate (FSI used) at the moment when the initial shock wave reaches the wall. The figure shows that the pressure peaks felt at the plate centre are notably smaller when the solid boundary deforms and absorbs part of the energy.

To further study the effect of material compliance on the pressure loading, we consider compliant materials with the properties described in §2.5. Due to the low compressive and shear strength of the considered compliant materials, a smaller collapse driving pressure, *P*_{d} = 0.1 MPa, is chosen as these materials cannot handle the large *P*_{d} value used with the metallic alloys. Figure 18*a* shows a comparison of the time histories of the pressures resulting from the bubble dynamics for the same base case shown in the previous section. The figure shows for illustration the results on a rigid fixed plate, on A2205, Al7075, Rubber no.2 and Versalink.

Figure 18*b* shows a zoom on the pressure peaks due to the bubble ring collapse. This figure clearly shows that the pressures felt by the metallic alloys and the rigid plate are very close. On the other hand, due to absorption of energy through deformation, lower magnitude of the impact pressures are felt by the softer materials. This was already observed experimentally by Chahine & Kalumuck [85]. Also, a delay in the peak timing indicates that a soft material tends to elongate the bubble period as compared with the rigid material as shown in figure 19.

## 4. Second bubble effect on bubble collapse near wall

As in the previous sections of this paper, studies of the interaction of inertial bubbles with nearby boundaries have mostly concentrated on isolated bubbles. However, bubbles are seldom isolated and their interactions with each other affect their dynamics and the way by which they load a nearby boundary. Also, there are configurations, such as for controlled cell sonoporation [30,31], or for underwater explosion applications [86,87], where the interactive dynamics of two or more bubbles is purposely sought.

In this section, we consider for illustration a simple example of two identical inertial bubbles initially located as illustrated in figure 20. The two bubbles are initially spherical and at the same normalized standoff distance, from the plate. Their centres are separated from each other by a distance, *D*. The initial bubble radius for each of the bubbles is *R*_{0} = 50 µm and the initial gas pressure (*P*_{g0} = 1500 GPa) is specified such that the bubble will grow to a maximum radius *R*_{max} of 2 mm if it was isolated in a free field when the ambient pressure (bubble collapse pressure) is *P*_{d} = 100 MPa, for example for a bubble generated by energy deposition in a the high stagnation pressure of a cavitating jet. The very high initial pressure was obtained by applying the Rayleigh–Plesset equation to obtain an initial bubble radius of 50 µm and a maximum value of 2 mm. This would correspond to a laser-generated bubble where actually the value of *R*_{0} and *P*_{g0} do not really correspond to a physical configuration but are selected to represent mathematically the energy deposit in the laser beam, which generates a 2-mm bubble.

Here, we consider the effects of the normalized bubble spacing, on the pressure loading and the pit shapes for the same Figure 21 presents for three different values of the normalized bubble spacing, 1.5 and 2.5, the bubble shapes at a time when the re-entrant jet (appearing as a darkened area of the grid in figure 21) is about to impact the opposite side of the bubble. The jet angle is obtained from the velocity vector at the fastest node on the bubble re-entrant jet surface. While a single bubble collapsing near a rigid wall forms a re-entrant jet perpendicular to the wall, in the two bubbles problem the re-entrant jet direction is controlled by both the presence of the second bubble and the wall. Therefore, the jet forms a shallower angle, between 20° and 70° for the various bubbles spacing.

To compare quantitatively the effect of the spacing, , on the re-entrant jet behaviour, we introduce as a jet characteristic the momentum averaged jet velocity vector, **V**_{mom}, defined as the average of the liquid velocity integrated in the whole re-entrant jet volume (region of the bubble surface with negative curvature) as follows:
4.1where **V** is the liquid velocity at any point inside the jet, and is the jet volume.

The jet speed versus time for the three values of can be seen in figure 22. We can observe in the figure that the presence of other nearby dynamically behaving bubbles does not necessarily enhance the cavitation dynamic effects and can negatively affect the strength of each bubble re-entrant jet impact at the wall. These effects depend strongly on the bubble and wall spacing as well as on the timing of the two bubbles and deserve a full study.

Figure 22 shows that the single bubble re-entrant jet is the shortest lived and the fastest, reaching sound speed under the present conditions. The other cases are subsonic, with the condition resulting in the lowest jet speed with almost one-third of the isolated bubble jet speed.

Figure 23 analyses further the effect of the bubble spacing on the jet velocity at touchdown. It can be seen from the figure that the jet speed of the tandem bubbles at touchdown approaches that of a single bubble when increases. **V**_{mom} is only 5.5% less than the single bubble value for and 2.5% for This is consistent with previous work on the limit of influence of a rigid wall (where the image plays the role of the second bubble), which corresponds to a standoff from the wall of about 3 [88].

Figures 24, 25 and 26 show the pressure contours at different instances during the bubble collapse near the wall for the three normalized standoffs respectively. Each time sequence is shown with its own pressure scale in order to be able to highlight the important high-pressure regions locations and the displacement of these ‘hot spots' with time.

When the two bubbles are very close, (figure 24), the re-entrant jets are so much deviated from the direction of the wall that they impact on each other (and not on the bottom wall) in the plane of symmetry of the problem. A high impact pressure is then generated away from the bottom wall. Figure 24*a* shows the initiation of the shock when the two re-entrant jets impact each other. The pressure then reaches in figure 24*b* its highest peak when the top halves of the two bubbles complete their collapse. A high-pressure wave is then observed to propagate radially away from the highest pressure spot with a pressure front moving towards the wall. The pressure amplitude is however significantly attenuated before it reaches the bottom wall in figure 24*c*.

For each of the two bubbles is equidistant from the rigid wall and the plane of symmetry. This equal attraction results in two very wide re-entrant jets both moving towards the centre of the plate. Figure 25*a* shows the complex shape of each bubble before touchdown on both the plane of symmetry and the bottom plate as seen in figure 25*b*. This generates high-pressure spots at the three impact locations: two on the bottom plate and one in the plane of symmetry. As the re-entrant jet front continues moving towards the plate centre, the residuals of bubbles near the centre eventually collapse completely and result in a highly focused pressure loading at the centre of the plate.

In the other extreme case, (figure 26), the interaction between the bubbles is much weaker, the re-entrant jets remain close to perpendicular to the wall, and two shock centres are generated first when each jet impacts the wall. These shocks then move along the wall surface towards the plane of symmetry. The pressure loading on the wall is then much higher than for the case.

A quantitative comparison of the time history of the pressure at the plate centre is shown in figure 27. It is seen that results here in a highest pressure loading among the three tandem bubble cases at the centre of plate. However, when compared with the single bubble case, the pressure is only about 88% of the isolated bubble case. Figure 28 shows a comparison of the shape of the final pit in Al 7075 for the three bubble spacing cases. It is seen that the *D* = 1.5 *R*_{max} results in the deepest pit among the three tandem bubble cases because of the highest concentrated loading at the centre of plate. However, the case results in a wider pit because of the two separate high-pressure impacts from each of the bubbles. The single bubble case is however the most erosive for the present conditions. Figure 29 shows a comparison of the volume of the final pit in Al 7075 for the three bubble spacing cases. It is seen that the total volume of the pit, as for the depth, increases as the spacing is increased and it surpasses the sum of the pit volumes from the two isolated single bubbles for case (figure 29).

Figure 30 investigates another tandem bubble configuration to illustrate conditions where the two bubble sizes and the two standoffs are not the same. The first bubble is selected to have an initial radius *R*_{0} = 142 µm, and initial gas pressure *P _{g}*

_{0}= 115 GPa, and is located at a normalized standoff, , while the second bubble has

*R*

_{0}= 50 µm,

*P*

_{g}_{0}= 1500 GPa and The initial pressures are again selected such that each bubble would grow, if isolated and in an infinite medium, to a maximum radius

*R*

_{max}= 2 mm when the ambient pressure is

*P*

_{d}= 100 MPa The normalized spacing between the two bubbles is . As the second bubble, located in the right-hand side of the figure, collapses first, the incompressible–compressible link in this case was initiated right before the re-entrant jet of the second bubble touches down, as shown in figure 30.

Figure 31 shows pressure contours at different instances during the collapse of the two bubbles. It is seen that the right-hand side bubble collapses first at *t* = 15.55 µs and results in a high-pressure loading of about 1700 MPa on the material surface as seen in figure 31*c*. The left-hand side bubble collapses next at *t* *=* 17.92 µs and results in a pressure loading of about 1500 MPa as seen in figure 31*f*.

Figure 32 compares the resulting permanent deformation outlines in the plane of symmetry of the problem for Al 7075. The figure compares the tandem bubble results with each of the two bubbles when alone. It is seen that each of the two isolated bubbles results in about the same pit depth with the second bubble producing a slightly wider pit. Both single bubble pits are, however, much larger than when the tandem bubbles act together. The tandem bubbles generate two pits with the right bubble, which collapses first and generates a higher pressure loading, resulting in a wider and deeper pit on the material than the left bubble. The sum of the volumes of the two pits is about the same as that due to two of the isolated single bubbles. A systematic study on the various combinations of tandem bubbles sizes and position would need to be conducted in the future to further understand the various effects of bubble/bubble interaction (constructive versus destructive) on the permanent deformation.

## 5. Bubble cloud effects

In this section, we discuss, as an illustration of multi-bubble dynamics, the behaviour of an inertial cloud of bubbles interacting with a soft material, Polyurea. The selected three-dimensional computational domain is 100 × 100 × 50 mm and possesses two planes of symmetry (the *XOY* and *YOZ* planes) as shown in figure 33. The bubble ‘cloud’ behaviour is triggered by the dynamic of a bubble of initial radius 3.15 mm and initial internal pressure 25 atm placed in the domain which has an ambient pressure of 1 atm. The driving bubble is located on the *Z*-axis 20 mm away from the material surface whose initial interface with the liquid is located in the plane *XOY*. The material is a 6 mm thick block of polyurea (Versalink). The bubble cloud is composed for ease of analysis of bubbles with the same radius initially at equilibrium with the 1 atm ambient pressure.

Figure 34 illustrates such a cloud configuration in the one-fourth volume (because of double symmetry). The bubbles in the cloud are located within the projected square [−2 mm < *x* < 2 mm, −2 mm < *y* < 2 mm] and are equally spaced in the vertical direction within a 1 mm height [*z* < 1 mm]. The bubbles in the cloud are distributed in four layers: the bottom layer (closest to the Versalink) has 10 × 10 bubbles, the second layer has 8 × 10 bubbles, the third layer has 10 × 10 bubbles and the top layer has 8 × 10 bubbles in the quarter domain. Bubble centres in each layer are staggered by half bubble spacing from the adjacent layers. The total number of bubbles is 360 in the full domain. The finest mesh size of 25 µm is used near these bubbles. A total of 4.7 million cells are used to model the problem. In the far field, no reflection boundary conditions are applied. The simulations in this section were obtained from the compressible code only without resorting to the compressible–incompressible link.

Figure 35 shows the pressure field near the bubble cloud 10, 19 and 21 µs after the initiation of the driving bubble growth and as the high-pressure wave reaches the bubble cloud. The cloud is seen to first act as a shield, preventing the high pressure generated by the source bubble from reaching the wall. This is seen in figure 35*a*,*b* where the high pressure (red region) impacts the wall outside the bubble cloud while the pressure under the bubble cloud remains close to the initial 1 atm. As first observed in [89], the bubbles on the outer edge of the cloud begin to collapse first, and then the phenomenon propagates inward to the centre of the bubble cloud. The bubbles in the cloud then collapse strongly (figure 35*c*) and then go through multiple collapse–rebound cycles, which generate multiple pressure pulses.

Figure 36 shows the time history of the radii of a few selected bubbles in the cloud. Overall, the bubbles in the cloud practically collapse and rebound almost simultaneously. The first collapse is the strongest as this brings the bubbles to their smallest volume during their history (figure 35*c*). Figure 36 also shows the internal pressure of the driving bubble, and the radii of the bubbles in the cloud when they all have initial radii of 150 µm. Even though the outside contour of the cloud has not changed from the previous case, the cloud period is increased–in the same way as the individual bubble radii following the Rayleigh collapse time—and is almost doubled. Here again all bubbles collapse and rebound almost in unison resulting in a very high pressure at collapse. Figure 37 shows the corresponding pressure time history at the origin of coordinate, i.e. centre of the material surface. The first pressure peak from the collapse of the 150 µm bubbles is more than four times higher than that of the 100 µm bubbles.

Several computations with a different number of bubbles in the cloud are illustrated in figures 38 and 39. Figure 38 shows the time histories of the pressure at the centre on the material obtained with different number of bubbles in the cloud both near a rigid wall and when the wall is responsive and made of Versalink. It can be observed that the period of the bubble cloud collapse and the pressure peak both increase as the number of bubbles in the cloud increases. This is the case for both material responses: one with a rigid wall and the other with a compliant Versalink wall. The bubble cloud works as an energy accumulator; the energy from the incoming pressure wave is absorbed by the bubble cloud, and then re-emitted as the bubble cloud collapses at a later time. A larger cloud with a larger number of bubbles accumulates more energy and then generates later a larger pressure peak at the cloud collapse.

Figure 39 shows four curves selected from figure 38. The pressure peaks with the rigid wall are observed to be two to three times higher than those with the compliant wall. As the high pressure reaches the compliant wall, the wall responds by compressing first and absorbing energy through deformation. This fluid–structure interaction and energy absorption by deformation reduces the pressure peaks.

## 6. Conclusion

In this contribution, the material pitting due to cavitation bubble collapse is studied by modelling the dynamics of growing and collapsing bubbles near responding and deforming materials. The pressure loading on the material surface during the bubble collapse is shown to be due to the re-entrant jet impact and to the collapse of the remaining bubble ring. The high-pressure loading results in high stress waves, which propagate radially from the loading location into the material and cause material deformation. Permanent deformation in the shape of a pit is formed when the local equivalent stresses exceed the material yield stress.

The impulsive pressure loading due to the bubble collapse is highly dependent on the initial standoff distance between the bubble and the nearby boundary. This standoff distance affects the jet characteristics in a non-monotonic fashion. Higher jet velocities occur at the larger standoff distances. However, a higher jet velocity does not necessarily result in a higher impact pressure, as the impact pressure also depends on the distance the jet front has to travel to impact the wall after touching down the opposite side of the bubble. The energy transferred to the wall is maximum at a normalized standoff distance close to . A more concentrated pressure loading on the material surface is obtained for smaller standoffs where both the jet touches downs and the bubble ring collapses very close to the wall. Such concentrated pressure loadings result in deeper but narrower pits. As a result, the shape of the pit, i.e. the ratio of pit radius and depth does not vary monotonically with standoff.

The magnitude of the pressure peak felt by the material depends on the response and amount of deformation of the solid. Fluid–structure interaction simulations show that the load on the material is damped with material deformation. The load reduction with wall response increases when the solid boundary deformation increases due to increased energy absorption. Impact pressures for metallic alloys are very close to those on a rigid plate while compliant materials deform and absorb energy. This results in lower magnitude of the impact pressures for the coatings and delays in peak occurrence due to lengthening of the bubble period.

Interaction between bubbles significantly influences the pressure loading on the material surface and the resulting pit shape. This interaction requires extensive study as both enhancement and negative interference of the interaction on the resulting damage can be seen. A classification of these effects deserves further investigations. Similarly, the work needs to be extended to the study of actual damage and material loss, which requires fracture and damage models as opposed to the present relatively simple elastic-plastic model used here. Such models would involve various criteria including tension and shear failure such as for hardened and brittle materials and heating effects such as for coatings.

## Competing interests

We declare we have no competing interests.

## Funding

This work was conducted under support from Dynaflow, Inc. internal IR&D and partial support from the Office of Naval Research under contract N00014-12-M-0238, monitored by Dr Ki-Han Kim.

## Acknowledgements

We thank Dr Kim for his support. We thank Mr Gregory Harris from the Naval Surface Warfare Center, Indian Head, for allowing us to access the Gemini code and giving us the opportunity to contribute to the coupling within Dysmas of our code 3DynaFS and the Dyna3D Structure code. We are also grateful to many colleagues at Dynaflow, who have contributed to several aspects of this study, most particularly, Dr Jin-Keun Choi and Dr Anil Kapahi for their contributions.

## Footnotes

One contribution of 13 to a theme issue ‘Amazing (cavitation) bubbles: great potentials and challenges’.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.