Cavitation over solid surfaces: microbubble collapse, shock waves, and elastic response

We discuss the interaction of the strongly nonlinear fluid motion induced by the collapse of a vapor microbubble over a planar surface and the elastic dynamics of the underlying solid. The fluid is described using an extension of the Navier-Stokes equations endowed with distributed capillary stresses in the context of a diffuse interface approach. The collapse of the bubble is triggered by overpressure in the liquid and leads to an intense jet that pierces the bubble, changing the bubble topology from spheroidal to toroidal, and impinges the solid wall inducing an intense and strongly localized load. Moreover, at bubble collapse, a compression wave is launched into the liquid surrounding the bubble. By propagating along the solid surface, the compression wave combined with the liquid jet excites the dynamics of the elastic solid, producing a complex system of waves, including, longitudinal, transversal, and Rayleigh waves, propagating in the solid. It is conjectured that the intense deformation of the solid induced by the strongly localized liquid jet may lead to the plastic deformation of the solid producing the surface pitting observed in many applications subject to cavitation-induced material damage.


Introduction
The subject of cavitation is rather wide and may be taken to range from nucleation of a vapor phase in a liquid, due to pressure decrease below a threshold, to the growth and coalescence of bubbles, and their dynamical coupling with the fluid motion [1]. The overall process is strongly influenced by the environment where it takes place, whether in free space or in presence of boundaries of different nature (e.g. rigid, deformable, soft, elastic, plastic). The confinement effect may play a big role as shown by a number of instances where cavitation occurs in narrow cavities found, e.g., in botany for the fern sporangium [2,3] or for plant xylems [4,5], in dentistry for root canal Abstract We discuss the interaction of the strongly nonlinear fluid motion induced by the collapse of a vapor microbubble over a planar surface and the elastic dynamics of the underlying solid. The fluid is described using an extension of the Navier-Stokes equations endowed with distributed capillary stresses in the context of a diffuse interface approach. The collapse of the bubble is triggered by overpressure in the liquid and leads to an intense jet that pierces the bubble, changing the bubble topology from spheroidal to toroidal, and impinges the solid wall inducing an intense and strongly localized load. Moreover, at bubble collapse, a compression wave is launched into the liquid surrounding the bubble. By propagating treatment [6], in ophthalmic surgery [7] or for cavitation assisted drug delivery [8].
Among the different facets of the problem, bubble collapse is probably one of the most spectacular and will be the subject of the present paper focused on the strong effects associated with bubble implosion close to a planar, elastic wall. This configuration is of paramount importance for hydraulic turbomachinery and ship propellers, for example, where cavitation bubble implosion on the blades is the main cause of erosion [9,10], and in bearings, where pitting is a common and disruptive phenomenon [11][12][13]. Given the great practical and economical relevance, the issue raised considerable attention and sophisticated experimental set-ups were conceived to analyze and understand this complex phenomenon [14][15][16][17].
Before entering the details of this specific problem, it is probably worth recalling what is the simplest, yet extremely rich description of bubble dynamics in free space. We owe Rayleigh [18] the derivation of the fundamental equation describing the dynamics of a spherically symmetric bubble immersed in a liquid which, in its modern form, became known as the Rayleigh-Plesset (RP) equation [19]. This model equation is so simple and illuminating that its derivation may be worth being briefly recalled before digging deeper into the subject of the paper. Assuming spherical symmetry, the only non-vanishing, radial velocity component u r of the unbound (incompressible, constant mass density 0 ) liquid enclosing a bubble of radius R can be expressed in terms of the potential (r) = ∫ ∞ r u r dr . Since mass conservation requires u r =ṘR 2 ∕r 2 , the velocity potential follows at once as (r, t) =ṘR 2 ∕r , allowing to integrate the Navier-Stokes equation from the bubble interface (where u r =Ṙ ) to infinite, where the pressure p(r, t) is taken to be p ∞ (t) . This procedure leads to the RP equation, In the above equation, the viscous term of the Navier-Stokes equations, ∇ 2 vanishes altogether on account of radial symmetry = −∇ combined with incompressibility ∇ 2 = 0 while the pressure at the bubble interface is prescribed by the force balance T rr (R) + 2 ∕R = P i , where T rr (r, t) = p(r, t) − 2 u r ∕ r is the force per unit area exerted by the liquid on the bubble which, combined with the surface tension , balances the pressure P i of the internal gaseous phase. Note that only through this boundary condition does the liquid viscosity play a (dissipative) role. The behavior of the bubble also depends on the thermodynamics of the gaseous phase. Two simple models are usually considered, either the pressure inside the bubble is constant, given by the saturated vapor pressure for a vapor bubble, or the gas trapped inside is assumed to undergo a polytropic transformation with exponent , P i ∕P 0 = R∕R 0 3 . Figure 1 illustrates the RP evolution of the radius of a spherical bubble, subject to a stepwise compression. For small compressions, the bubble radius oscillates around a new equilibrium with a slightly reduced radius. The small sinusoidal oscillations obey the linearized RP equation around the new equilibrium condition. However, as the compression becomes more intense, while the oscillation period decreases as a consequence of the increased stiffness, significant non-linear effects show up in the sequence of compressions and re-expansions, with more gentle crests and sharper through. According to the model, the velocity of the liquid around the bubble oscillates following the bubble radius, see the expression of u r (r, t) provided above. At increasing compression, the minimum bubble radius gets smaller, the gas inside gets compressed and the temperature increases. The thin line with symbols provides Ṙ for a mild compression (panel (a), yellow symbols) and a strong compression (panel (b), green symbols), respectively. Apparently, as the compression intensifies, Ṙ becomes less symmetrical in time until, panel (b), something more and more similar to a sudden jump from negative (end of the bubble compression phase) to positive velocity (beginning of the subsequent expansion phase) develops. This is called the rebound: the bubble cyclically collapses to a minimum radius and suddenly expands back. In these conditions liquid compressibility is no longer negligible and, as discussed below, shock waves are launched in the liquid, contributing to the energy loss of the bubble together with viscous dissipation and heat conduction from the bubble, getting hotter due to compression, to the surrounding liquid.
During the collapse, pressure and temperature inside the bubble may attain extremely large values, leading to the emission of light, the so-called sonoluminescence [20], where acoustics is implied in the name since bubble collapses and rebounds are usually driven by ultrasound irradiation.
The collapse phenomenology becomes more complex if the spherical symmetry is broken by a planar solid wall. In this case, a strong jet is generated by the collapsing bubble which together with the shock waves propagating in the liquid and impinging on the solid surface excites elastic oscillations in the material.

Dynamics of a collapsing bubble
In this context, the paper aims at describing the highly non-linear dynamics of a micron-sized (vapor) bubble collapsing on top of an elastic wall, accounting for the propagation of elastic waves in the underlying solid. The overall phenomenology is illustrated in Fig. 2, which reports a snapshot taken from one of the numerical simulations to be described below. A micro-bubble is initially in unstable equilibrium near a planar boundary between the fluid domain (the upper half-space) and the elastic domain. As for the prototype problem with the RP equation, the bubble collapse is triggered by the sudden compression of the liquid. The compression induces the axisymmetric shrinkage of the bubble, much like for the RP equation. Here, however, the impermeable planar wall restrains the receding motion of the bubble portion proximal to the wall and consequently intensifies the motion of the distal region (to a first approximation, the compression produces a volume change with the distal part compensating for the reduced mobility of the proximal part). The overall process amounts to the formation of a jetting motion toward the wall that perforates the bubble, inducing its topological reshaping into a toroid. As the bubble keeps shrinking the temperature of the vapor increases due to compression, until reaching supercritical conditions. After becoming incondensable, the supercritical gaseous phase keeps on compressing, until the converging motion toward the bubble gets suddenly stopped, inducing the emission of strong compression waves in the liquid.
The top part of Fig. 2 is a visualization of the fluid field. The grey contours are contour levels of the fluid density gradient. The gradient intensity concentrates at the bubble interface and in the (almost) spherical to the liquid. In panel a, the different curves correspond to ΔP∕P 0 = 10 −1 , 1, 10 , blue, yellow, and purple, respectively. In panel b ΔP∕P 0 = 10 1 , 10 2 , 5 × 10 2 , for the purple, green, and yellow curve, respectively. Due to compression, the natural oscillation frequency increases. In the two panels, the thin line with symbols represents the (normalized) radial bubble velocity to be read on the right axis which corresponds to the case reported with the same color compression wave traveling outward. The impact of the liquid jet and the moving pressure disturbance on the fluid/solid interface excites the response of the solid, as described in the next section. The left part of the upper plot reports the fluid velocity vectors. Concerning the fluid phase, the simulation illustrated in Fig. 2 is based on a diffuse interface model [21][22][23] of the fluid which is able to account for the liquid/vapor phase transition, including thermal effects associated with the latent heat [24,25].
with ( , t) the density, ( , t) the velocity vector, E( , t) the total energy density of the system -comprising internal and kinetic energies. The stress tensor ( , t) and the energy flux e ( , t) are expressed as [25] where the capillary coefficient ( ) is a function of the temperature , and ̃ are the first and second viscosity coefficients, respectively, and k the thermal conductivity. In Eq. (3) the usual viscous stress tensor is augmented by distributed capillary effects concentrated in the narrow interfacial layer and in Eq. (4) the energy flux is generalized by adding a capillary contribution to the Fourier law. System (2) is closed using the van der Waals equation of state (EoS) for the pressure p 0 , which, in reduced variables, i.e. normalizing pressure, temperature, and density with the corresponding critical values, p c = 22 MPa, c = 647 K, c = 322 kg∕m 3 , reads and numerically solved in the cylindrical domain Ω = [0, R max ] × [0, H] by the finite element method [26] assuming axisymmetry after its rearrangement in a weak form. With regard to the fluid phase, the solid wall W is assumed rigid, adiabatic, and neutrally wettable [27], implying the boundary conditions = 0 , Given a cylindrical coordinate system (r, , z) with the axis orthogonal to the solid wall placed at z = 0 , the fluid occupies the upper half-space ( z > 0 ), and the system is initialized with a bubble of radius R = 1 m and center at Z c , as shown in the sketch in panel (a) of Fig. 3. Viscosity, thermal conductivity, and surface tension of (liquid) water are assumed, and the reference The unstructured computational grid is adaptevely refined during the simulations, starting from an initial maximum resolution at the liquid/vapor interface of Δr = Δz = R max ∕2 15 and an initial minimum resolution in the unperturbed liquid of Δr = Δz = R max ∕2 6 (being R max = 20R ). To appropriately resolve the dynamics of the bubble, a portion of the domain Ω res = [0, R max ∕2] × [0, H res ] is forced to maintain a minimum grid resolution of Δr = Δz = R max ∕2 11 . The time step Δt ranges between 10 −5 and 10 −3 , depending on the collapse phase and initial ΔP∕P 0 .
Initially, the bubble is in unstable equilibrium with the pressure in the liquid, according to the Young-Laplace law p v ( v , ) − p l ( l , ) = 2 ∕R , where p v∕l are the vapor/liquid pressures and where x n is the radial coordinate from the bubble center, the surface tension of the liquid/vapor interface [28,29], with = 0.5 c , p l = 0.0211978 p c , v = 0.0217468 c and l = 2.45824 c . Note that, Eq. (6) approximates the surface tension of the bubble with that of a planar interface, see [28,30] for the expression accounting for the interface curvature.
Initially, the liquid is in metastable conditions, i.e., in the pressure/specific volume diagram, the representative point lies in the region below the binodal line and above the spinodal, see panel (b) of Fig. 3. At the time t = 0 the pressure of the liquid environment is suddenly increased to p l + Δp , triggering the successive collapse dynamics of the bubble.
An example of the bubble evolution is provided in Fig. 4, where the bubble configurations at different instants is identified as the density isolevel (r, z, t) = 1 . During the collapse, the bubble shrinks, a liquid jet pierces the bubble and the topology of the bubble changes from spheroidal to toroidal. This Given the state of the liquid, the vapor will be at the same temperature and a density such that the pressure jump across the bubble interface satisfies the Young-Laplace law (note that the pressure in the stable vapor is larger than in the metastable liquid) sequence of events is generic and can be compared with experimental data, typically obtained for millimeter bubbles, [15,16]. For given thermodynamic conditions, the initial bubble distance from the wall strongly influences the non-sphericity of the collapse, which from perfectly spherical far from the wall, becomes more and more asymmetric [25,31]. As anticipated, a more detailed view of a snapshot along the evolution is provided in Fig. 2 corresponding to the time t = 22 , compare with Fig. 4. There, a strong compression wave irradiated by the bubble is visible in the contour of density gradient magnitude.
The combined action of the fast jet impinging the wall and the propagation of a fast compression (shock) wave traveling along the surface excite the dynamics of the elastic substrate. The normal load distribution applied to the solid surface is shown in Fig. 5. The instantaneous local loads are significantly strong and change in time at a fast rate, given the propagation speed of the compression wave traveling on top of the solid and the impulsive nature of the load generated by the liquid jet. Discussing the ensuing dynamics of the elastic solid is the subject of the foregoing section.

Elastic waves in the substrate
A snapshot of the system of waves propagating in the solid occupying the lower half-space z < 0 is shown on the bottom part of Fig. 2, which in its left part displays the norm of the elastic displacement, | (r, z, t)| and in the right part the norm of the deviatoric part of the elastic stress T vM = √ 3∕2� e − Tr � e � ∕3� -the socalled Von Mises stress. T vM is substantially localized Envelopes of the maximum normal stress transmitted to the solid surface by the traveling pressure disturbance induced by the collapsing bubble. The envelopes are strongly influenced by the distance between the bubble and the solid. As the distance decreases, the impact gets more concentrated, with increased maximum pressure. In the inset, the envelope corresponds to Z c ∕R = 0.2 together with three instantaneous configurations of the time dependent load distribution (red curves). The peak load propagates toward larger r. All load profiles are taken at constant Δp∕p l = 28.5 in the region where the jet is impinging the elastic wall. This information is important for the understanding of the damage potential of the imploding bubble. Indeed, certain elastic materials, e.g. metals, undergo a plastic transition when the local Von Mises stress exceeds a threshold, leading to permanent deformation of the surface, known as pitting, see the Discussion & Conclusions section below and [31] for an in-depth analysis of the plastic response under the dynamic loading of a collapsing bubble. The right part of the figure highlights a complex system of wavefronts propagating in the elastic medium. Two features are prominent, a (hemispherical) wave adjoining the compression wave propagating in the liquid phase at the fluid/solid interface, and the trace of a cone, appearing as the oblique line, with vertex at the common trace on the fluid/solid interface of the spherical elastic wave in the solid and the compression wave in the liquid. In order to better understand the wave features of the solution, we have to resort to the classical equations of (linear, isotropic) elasticity, where (r, z, t) is the elastic displacement and the elastic stress is e = 1 ∇ ⋅ + 2 ∇ + ∇ T , with 1 and 2 the first and second Lamé coefficients -related to the shear and Young moduli by the relations G = 2 and E = 2 3 1 + 2 2 ∕( 1 + 2 )and S the mass density of the elastic solid, assumed constant in the linearized momentum equation. In fact, the solid density is S + Δ S , where the density variation is a field Δ S (r, z, t) which changes with the local compression/expansion of the material, obeying the (linearized) mass conservation equation, Δ S ∕ t + S ∇ ⋅̇ = 0 . At the time t = 0 the elastic material is undeformed, i.e. the elastic displacement is zero, (r, z, 0) = 0.
Concerning the boundary conditions, in principle, the problem should be treated as a free boundary problem, where both the displacement, hence the velocity, and the contact force per unit surface (the tension) are required to be continuous across the fluid/solid interface W = ∶ z = z I (r, t) , where z = z I (r, t) is the Monge representation of the surface, Clearly, the interface configuration is recovered by solving the kinematic equation where = =̇ is the velocity of a geometric point of the interface and = � 1, z I ∕ r � ∕ √ 1 + ( z I ∕ r) 2 is its unit normal, oriented, e.g., toward the fluid domain. In fact, for a stiff solid, the interface displacement is so small that the boundary condition can be safely linearized by assuming an unperturbed configuration of the domain, I = 0 , that amounts to assuming the wall as being rigid, as implicitly assumed when dealing with the fluid problem in Sect. 2. Moreover, after solving the fluid problem one realizes that the tangential stress exerted by the collapsing bubble on the interface is much smaller than the normal one, which basically amounts to the contribution of the pressure field.
In these conditions, the elastic problem is posed in the undeformed lower half-space, z < 0 , and the boundary condition reduces to requiring on the plane z = 0 , where p w (r, t) is the pressure exerted by the fluid on the solid surface, Fig. 5. The (small) deformation of the fluid/solid interface is then recovered after the fact as I = (r, 0, t).
After recasting the problem in a weak form, Eq. (7), with the above initial and boundary conditions, is solved with the finite element method [26] in the cylindrical computational domain Ω e = [0, R max ] × [−H e , 0] , with R max = H e = 10 R , discretized with a uniform grid 2048 × 2048 . The time step adopted, Δt ≤ 0.9Δr∕ max (c L , v D ) , allows the resolution of the fastest timescale in the system. At variance with the fluid problem, where viscous dissipation dampens the acoustic waves, the elastic system sustains the undamped propagation of waves. In order to avoid spurious reflections at the boundary of the truncated computational domain, Ω e , Rayleigh damping layers, where the equations are modified by including the artificial term −̇ to the righthand side, are added in correspondence to the side, r = R max , and bottom, z = −H e artificial boundaries, [32]. The thickness and the damping coefficient were tuned to dissipate the waves before they are reflected back in the physically relevant region.
As anticipated, the snapshot of Fig. 2 shows a complex pattern of elastic waves which are excited by the load traveling on the fluid/solid interface. It is well known that elastic materials support the propagation of two fundamentally different wave systems: one system of longitudinal waves and two systems of transverse waves. In the case of our axisymmetric configuration, the transverse system with oscillations in the plane orthogonal to an axial section is inhibited by the symmetry. The various systems of waves are better isolated by decomposing the displacement field in terms of a scalar, , and a vector, , potential (Helmholtz decomposition), where, for an axisymmetric system, ∇ =̂ r ∕ r +̂ ∕r ∕ +̂ z ∕ z , with ̂ r , ̂ and ̂ z the radial, circumferential and axial unit vector, respectively. The vector potential is given by = (r, z, t)̂ , and the scalar potential is independent of the angular coordinate, = (r, z, t) . The equations for the scalar and vector potential are readily obtained after realizing that ∇ ⋅ = ∇ 2 and ∇ × = ∇ × (∇ × ) , where a judicious choice of the gauge allows to assume ∇ ⋅ = 0 . Indeed, taking the divergence and, respectively, the curl of Eq. (7) leads to It is immediate to realize that free waves described by the scalar potential are longitudinal, whereas those described by the vector potential are transversal (given the above equations, corresponding wave equations are obtained for the potentials in free space). In Eq. (12) c L = √ ( 1 + 2 2 )∕ S is the speed of free longitudinal (or compression) waves, while c T = √ 2 ∕ S is the speed of free transversal waves (sometimes called shear waves).
In a domain with a boundary, the two systems of waves are coupled via the boundary conditions, Eq. (10). If the applied pressure were spatially and temporally constant (e.g., p w = 0 ), the coupling of the longitudinal and the only transversal wave existing in an axisymmetric system gives rise to two other additional kinds of waves. One is the Rayleigh wave, which is, in fact, a surface wave confined to the surface layer. The propagation speed c R of the Rayleigh wave obeys the dispersion relationship, [33], which, for ductile metals, can be approximated as c R = (0.87 + 1.12 )∕(1 + ) c T , with the Poisson's ratio, [34]. The other additional kind of wave is the so-called von Schmidt, or head, wave, which can be roughly understood as the Mach cone of transversal waves generated by the interaction with the constant pressure interface of the faster longitudinal wave, which exceeds the transversal waves limit velocity, [35]. Different wave components excited by the load are disentangled in Fig. 6 which shows in the upper part a sequence of snapshots depicting the configuration of the ∇ ⋅ field. The bottom part provides the corresponding configurations of the ∇ × ⋅̂ field. In each panel, the vector field represents the displacement . For the case reported in the figure, the material parameters of the solid are selected to have c T < c a < c L , where c a = √ ( p∕ )� s , with the derivative taken at constant entropy, is the speed of sound in the liquid. In the top panels of Fig. 6 the front of the fastest longitudinal wave envelops the region where the solid is deformed. The bottom panels identify the front of the transversal system. It starts at the fluid/solid interface and appears as the section of a cone with the axial plane that progressively turns into a spherical front. The latter corresponds to the free propagation of a transversal wave. The former is generated by the intersection of the leading edge of the compression wave propagating in the liquid with speed c a with the fluid/solid interface. This can be described as a load discontinuity radially advancing along the interface with a velocity faster than the transversal speed of sound in the solid, c T . In the gas-dynamic jargon, this would correspond to a Mach cone, enveloping the region that can be reached by the transverse waves that were excited in the past instants by the moving front of the load. A careful look at the displacement field shows that across the cone the elastic displacement changes abruptly.
The Rayleigh and the von Schmidt waves are hardly seen in Fig. 6 but are clearly spotted in Fig. 7. In panel a, where the longitudinal field is depicted, the Rayleigh wave is clearly visible in the thin layer below the surface of the solid, located at r ≈ 35 , decaying exponentially in z. In panel b, where the transversal field is represented, in addition to the Rayleigh waves, the von Schmidt head wave is visible, whose oblique wavefront originates in the interaction of the longitudinal wave with the free surface ( r ≈ 80 ), and matches the shear wavefront in the bulk There, the longitudinal (a) and transverse (b) fields show a region of deformation limited to a thin layer below the interface. This corresponds to a comparatively weaker Rayleigh wave, which is traveling at a speed smaller than, but close to the speed of the transverse waves. The head wave instead is found in the zone ahead of the transverse disturbance generated by the load, in panel (b). The aperture angle of its Mach cone (straight boundary of the darker blue region on the rightmost part) is slightly, but detectably, smaller than that produced by the load (straight boundary of the reddish region fading to white) since c L ∕c T > c a ∕c T .

Discussion and conclusions
The collapse of a cavitation bubble close to a planar surface produces, as shown, extremely rich and complex effects. Although a basic understanding of the process can be gained by using the celebrated Rayleigh-Plesset equation, the presence of the solid wall makes the phenomenology much more complex, requiring to resort to solving the field equation for the fluid velocity. Here, we have discussed an extended version of the Navier-Stokes equations with distributed capillarity, accounting for the full thermodynamics of the liquid/vapor phase transition. The preeminent features observed in the simulations consist of the formation of a strong jetting motion toward the wall with the simultaneous emission of a system of intense compression waves. The intensity of the jet and the amplitude of the compression waves depend on a number of parameters, among which bubble size, bubble distance from the wall, and overpressure triggering the collapse are particularly important. Here, rather than exploring the parameter space of the system, we focused on one specific example to showcase the complexity of phenomenology in the case of a microbubble. From the point of view of the applications, the unsteady load induced on the underlying solid by the collapsing bubble is crucial. The extremely fast time scale over which the phenomenon evolves calls for solving the unsteady version of the equations governing the elastic solid. It is shown that the deformation is particularly significant in correspondence with the jet impingement region. The system of elastic waves generated by the jet and the propagation of the compression wave in the liquid, though obeying a linear system of equations, is far from being trivial. Different kinds of elastic waves undergo a coupled propagation forced by the time-dependent load. The details of the resulting pattern depend on the relative magnitude of the speed of the involved waves, compression (acoustic shock) wave in the liquid, longitudinal and transverse acoustic waves in the elastic solid, and Rayleigh-type surface waves. Space limitations did not allow us to discuss all the possibilities, see [31] for a more complete account. Here, to illustrate the richness of the phenomenology, we focused on the interesting case where the acoustic wave is slower than the longitudinal but faster than the transverse elastic waves.
We conclude by observing that the strong, very localized, and extremely fast deformation associated with the pressure wave impacting on the solid wall explains why the material damage appears in the form of localized pits formed just below the collapsing bubbles. Such effects can be computed considering an elastoplastic model for the solid as discussed in [31].