Numerical investigation of acoustic cavitation and viscoelastic tissue deformation

Highlights • A numerical method for acoustic cavitation and viscoelastic tissue deformation is presented.• Tissue deformation is maximum at shear modulus in the range of 1–10 MPa.• Tissue deformation generally increases with tissue bulk modulus and density.• Liquid jet formation has a significant effect on large tissue deformation with perforation.


Introduction
Acoustic cavitation, which can be practically defined as the growth and subsequent collapses of pre-existing bubble nuclei by ultrasound radiation [1,2], is a promising tool in a variety of engineering applications with including water decontamination [3,4], surface cleaning [5,6], and targeted drug delivery [7,8].The liquid jet generated during asymmetric bubble collapse mechanically destroys surrounding organic matter, including viruses, enabling high-quality wastewater treatment [4].The microstreaming, caused by a liquid jet spreading over the surrounding biofilm, creates localized shear stresses and removes contaminants attached to the underlying substrate [6].Oscillating mic robubbles can exert mechanical force on nearby tissues, promoting local drug uptake and improving drug penetration efficacy into brain or tumor tissues [7].Ultrasonic cleaning and drug delivery should be done while reducing substrate erosion and biological damage, and wastewater treatment should be performed to maximize inactivation of bacterial cells or viruses [6,7].Understanding the interactions between cavitation bubbles and viscoelastic tissues or cells is important for various applications of acoustic cavitation.
Extensive theoretical models for acoustic cavitation have been developed based on the Keller-Miksis (KM) equation including the liquid compressibility effect [2].Yasui [9] analyzed the influence of phase change due to thermal conduction on acoustic cavitation dynamics by modifying the KM equation and using a temperature boundary layer approximation.The results considering thermal conduction were found to be in better agreement with experimental data than the results without conduction.Storey and Szeri [10] investigated cavitation bubble collapse through a comprehensive analysis including the effects of phase change, diffusion and chemical reactions.They reported that chemical reaction of vapor inside the bubble can significantly lower the bubble temperature.Hong and Son [11] developed a general numerical model by modifying the theoretical equation and determining the evaporation rate directly from the temperature gradients at the interface, and investigated the cavitation threshold by varying noncondensible bubble nucleus radius and acoustic frequency.Merouani's group [12][13][14] also extended the KM model to include the phase change and chemical reaction effects.They analyzed the influences of methanol concentration and acoustic conditions on the generation of hydrogen and reactive species in methanol pyrolysis [12,13], and analyzed the water dissolution mechanism depending on the ultrasonic frequency and number of acoustic cycle [14].Although the theoretical models have been widely applied to understand acoustic cavitation phenomena, they have limitations in considering non-spherical bubble behavior and bubble-solid interaction.
Experimental researches on the interactions between single cavitation bubbles and nearby elastic boundaries have been performed using laser-induced or spark-induced bubble generation techniques [15][16][17].Brujan et al. [15] found that a mushroom-shaped bubble appears during the bubble contraction stage due to the repulsion of the elastic boundary deformed by the laser-generated bubble.Ma et al. [17] compared the bubble dynamics using spark-induced cavitation bubbles under two different boundaries with the same elastic modulus but different thickness.However, few experimental studies have been performed under ultrasonic cavitation conditions due to the difficulties in controlling the bubble generation position and size [18,19].
Alternatively, numerical simulations have been performed for bubble-tissue interactions under ultrasonic conditions.Fong et al. [20] analyzed the bubble growth and collapse near various biological tissues by combining a boundary element method based on inviscid and incompressible flow assumptions with a model that treats the biological tissue as an infinite fluid with elastic properties at the fluid-fluid boundary.They studied the influences of acoustic frequency and elastic modulus on the interaction between bubble and various biomaterials.Kobayashi et al. [21] numerically simulated shock-bubble interactions for various biomaterials by combining an equation of state to incorporate the compressibility of tissues and an improved ghost fluid technique for reliable calculations of compressible two-phase flows.However, their method has limitations in that it treats the biological tissue as a compressible fluid without elasticity.Guo et al. [22] studied the deformation of spherical cells depending on initial bubble generation size and location under ultrasonic conditions by modeling the elasticity of cell membrane with surface tension proportional to the surface area.The above numerical approaches to mimic biological cell or tissue have limitation in reflecting viscoelastic feature [23,24] and accurately predicting later stage of tissue perforation [25].
More realistic simulations have been conducted by calculating tissue deformation by the finite element method (FEM) with a tissue model that regards the tissue as an elastic solid.Chahine et al. [26] investigated the effects of secondary collapse of toroidal-shape bubble and liquid jet on the cleaning of soft materials by considering the tissue as a solid with elastoplastic properties and using the FEM to calculate tissue deformation.Zevnik and Dular [27] analyzed bacterial cell damage during bubble collapse by modeling the cell as an elastic solid with multilayered shell structures and using the FEM for cell deformation.However, only a few numerical simulations have been performed for the interactions between cavitation bubbles and viscoelastic solid tissues under ultrasonic conditions.
Full Eulerian methods [28][29][30][31], which are advantageous for describing large solid deformation without mesh regeneration [32,33], have also been applied to various biological tissues.Ii et al. [29] simulated the complex migration and deformation patterns of blood cells in practical intravascular situations by modeling the multiple red blood cells as biconcave capsules and employing the full Eulerian method.Recently, Koukas et al. [30] investigated the interactions between shock wave and bubble near biological tissues with two different shear moduli by extending the fluid tissue model to include the elastic features of solid tissues.Their computations successfully represented a variety of lithotripsy situations, neglecting the viscosity and surface tension effects which can be important in the acoustic cavitation dynamics of preexisting nuclei or microbubbles [27,34].Park and Son [31] performed a preliminary study on ultrasound-driven dodecafluoropentane (DDFP) bubble motion and tissue deformation, which is applicable to medical imaging and therapy, by improving a level-set (LS) method for bubble dynamics under ultrasonic conditions to include the elastic properties of tissue using the full Eulerian method.They investigated the effects of ultrasonic parameters, bubble-tissue distances and shear moduli on bubble-tissue interactions.However, their computations were limited to only thin incompressible solid tissues, which can not take into account ultrasound transmission through the solid regions.
In this work, a more general numerical method is developed for acoustic cavitaion bubble motion near a viscoelastic tissue by extending the LS method to treat the elastic properties and acoustic impedance of tissues as well as the effects of liquid and gas compressibilities.We also perform a more comprehensive analysis of acoustic cavitation and tissue deformation by varying tissue properties, such as shear modulus, bulk modulus and density, as well as bubble-tissue distances.

Numerical analysis
Fig. 1 shows the schematic for analysis of acoustic cavitation of a preexisting air bubble near a viscoelastic solid tissue.Computations are conducted under the assumptions that the fluid flow and solid motion axisymmetric, and all physical properties for each phases are constant except density.The compressibility effects of water and solid are treated by evaluating the water and solid tissue pressures from the Tait equations, where p ∞ = 101.3kPa and γ w = γ s = 7.15, regarding that the volume fraction of water is nearly 80% of the cell tissue [35].To more accurately consider the high-pressure bubble upon collapse, the bubble pressure is determined from the van der Waals equation for isothermal conditions, where Here, the subscript c denotes the critical state of air bubble and R g = 287 J/kgK.

Governing equations
The LS method described in our previous studies for various twophase flow applications [35][36][37] is extended to investigate acoustic cavitation and resulting tissue deformation.The unified governing equations for bubble, water and tissue regions are expressed as ∂B ∂t where Here, B, τ e , G are left Cauchy-Green deformation tensor, an elastic stress tensor for neo-Hookean materials, shear elastic modulus, respectively, and α is a smoothed Heaviside function that varies linearly along the interface thickness of 1h where h is a grid size.
The LS functions ϕ b and ϕ s for tracking the bubble-water and watertissue interfaces are advanced and reinitialized with the following equations: where t * is an artificial time for iteratively solving Eq. ( 15) and ϕ b/s are not reinitialized in the vicinity of the interface (|ϕ b/s | < 0.5h) for better volume conservation.
The governing equations for air bubble, water and tissue phases are spatially discretized in a staggered grid system.Further details for spatial and temporal discretization of the governing Eqs. ( 4)-( 6) and validation tests of the LS formulation for solid deformation and ultrasound-driven bubble motion have been addressed in our previous studies [31,37].

Computational conditions
In this calculations, we consider a stationary air bubble and nearby viscoelastic tissue immersed in ambient water at p ∞ = 1atm, referring to The sinusoidal ultrasonic pulse p u +p ∞ is applied to the upper boundaries of y = Y d,n as follows.
where A u and f u are the ultrasonic pulse amplitude and frequency, respectively.In the present computations, we keep A u = 0.3MPa, f u = 1MHz and γ s = γ w while changing l o , G, Π s and ρ s .The symmetric condition is imposed at r = 0 and the convective pressure condition at y = Y d,s .It is noted that the domain sizes of R d,e and Y d,s are chosen to be nearly four times the wavelength (a w /f u ) to ensure that the bubble motion is not disturbed by any reflection from the boundaries.

Acoustic cavitation
Calculations are first performed for acoustic cavitation bubble without including surrounding tissue at A u = 0.3 MPa and f u = 1 MHz.The results are presented in Fig. 2. The microbubble begins to grow at t = 0.32 μs when an ultrasonic pulse traveling through the ambient water at the speed of sound reaches the bubble surface.During the negative pressure pulse stage (0.32 μs < t < 0.82 μs), the bubble grows and pushes out the surrounding liquid, forming a source flow outward from the bubble center.While the positive pressure pulse is transferred to the bubble (t > 0.82 μs), the bubble expands further due to the fluid inertia around its surface and then begins to contract at t = 0.910 μs.As the bubble shrinks rapidly, a high pressure field is created near the bubble interface, and when its volume becomes minimum, a shock wave forms, as seen at t = 1.130 μs.Thereafter, the bubble collapses and rebounds into a non-spherical shape at t = 1.133 μs, propagating a shock wave through the ambient liquid.This is coincident with previous numerical studies showing that the pressure gradient and direction of the acoustic wave determine non-spherical bubble dynamics even in the absence of surrounding walls or obstacles [41], but the asymmetric nature of the flow field is not large enough to cause liquid jet formation.After re-expansion, the bubble returns to a nearly spherical shape with subsequent oscillations and reaches an equilibrium state due to the surface tension and liquid viscous effects.
For validation of the obtained results, we consider the following Keller-Miksis equation [42] for the spherical bubble motion under a plane pressure wave: where Here, the bubble pressure p b is obtained from Eq. ( 2) and the bubble mass balance of ρ b = ρ bo R 3 bo /R 3 b , and C w is the sound speed in liquid water.Fig. 3 compares our numerical results with the theoretical prediction from the KM Eq. ( 17).The volume-averaged bubble radius R b from the present simulation agrees well with the theoretical prediction for the growth and contraction period of the bubble.This indicates that  the present numerical model is consistent with the theoretical cavitation model when the bubble-tissue distance is long.

Influence of bubble-tissue distance
Fig. 4 shows the acoustic cavitation near a viscoelastic tissue and the associated liquid flow fields at l o = 5 μm, G = 1 MPa, Π s = Π w and ρ s = ρ w .While the bubble grows due to the negative acoustic pressure, it pushes away the surrounding liquid water and nearby deformable tissue.The compressed tissue rebounds due to its elastic restoring force and pushes the surrounding liquid upward at t = 0.920 μs, transferring its momentum to the bottom of the bubble surface.This perturbation propagates upward from the bubble bottom along the interface, forming a relatively fast velocity field around the bubble bottom, as observed at t = 0.990 μs.This results in a larger interfacial curvature in the bubble bottom part compared to the top part, as seen at t = 1.090 μs.As the bubble shrinks, the surrounding tissue expands upward by a sink flow developed toward the bubble center, which can be confirmed from previous numerical investigation [21] based on an improved ghost fluid method.The resistance of the underlying tissue causes the upper part of the bubble to contract more rapidly than the lower part, and an annular inflow parallel to the tissue boundary induces a large curvature at the top part, as depicted at t = 1.117 μs.The resulting asymmetric flow field between the bubble top and bottom parts induces a liquid jet formation and downward bubble migration.The bubble is pierced by a downward liquid jet at t = 1.201 μs as the bubble re-expands, but is not strong enough to cause significant deformation of the tissue.
Fig. 5 presents details of the bubble and tissue interfaces in Fig. 4 during the first period of bubble motion.While the bubble initially expands (t < 0.91 μs), it pushes out the surrounding tissue, causing the tissue to undergo compressive deformation until t = 0.78 μs, as presented in Fig. 5(a).As the bubble expansion rate decreases, the elastic restoring force in the compressed tissue exceeds the liquid momentum transferred to the tissue surface, raising the tissue position on the central axis by 30% of the maximum compressive deformation at t = 0.91 μs.During the bubble contraction phase (0.91μs < t < 1.12 μs), the tissue position returns to y = 0 at t = 0.94 μs and is gradually attracted toward the bubble.The bubble sides contract rapidly, which induces an elongated shape of the bubble along the vertical axis.Thereafter, as the bubble rapidly collapses with downward migration, the tissue expands to reach its highest position at t = 1.12 μs.During the bubble rebound phase (t > 1.12 μs), the downward movement of the re-expanding bubble once again lowers the tissue surface below y = 0 with rapid compressive deformation.
When l o decreases to 4 μm, the resulting flow fields are plotted in Fig. 6.As the liquid between the bubble and tissue interfaces is released during the bubble growth period and a thin liquid layer develops beneath the bubble, the bubble-tissue interaction becomes stronger, resulting in an inverted cone-shaped bubble at t = 1.090 μs, which was also observed in the previous experimental study [17].Due to the flow resistance of nearby tissue, a highly asymmetric velocity field develops along the y-axis and the upper part of the bubble surface rapidly contracts compared to the lower part.At t = 1.114 μs, the annular inflow towards the central axis impinges on the top of the bubble surface, forming an extreme interfacial curvature, which causes the liquid to be squeezed into the bubble in the form of a strong vertical jet.After the bubble collapse (t > 1.118 μs), the bubble migrates rapidly downward and almost touches the tissue surface, resulting in compressive tissue deformation accompanied by bubble re-expansion.As the bubble recontracts and the tissue undergoes expansive deformation toward the nearby bubble, narrowing the pore entrance, a small portion of the bubble becomes trapped inside the tissue, causing long-term tissue deformation, as depicted at t = 1.262 μs.Fig. 7 shows the influence of l o on R b and tissue surface location y sc at the central axis (r = 0) for G = 1 MPa, Π s = Π w and ρ s = ρ w .When an ultrasonic pulse traveling through the liquid water reaches the bubble at t = 0.32 μs, R b begins to increase from R bo = 1 μm.During the ultrasound pulse period, bubble contraction proceeds faster than bubble expansion, which can be observed in the previous cavitation models [10,14].This is because the ultrasonic pulsed positive pressure is applied to a larger bubble area during expansion.As the initial bubbletissue distance decreases from l o = 7 μm to l o = 4 μm, the maximum bubble radius R b,max decreases slightly by about 3.3%.During the first cycle of bubble growth and contraction, the tissue that resists the surrounding fluid flow further impedes bubble motion as l o decreases.
Thereafter, the bubble at l o = 4 μm re-expands faster than the other cases, and the second maximum bubble radius is about 34% larger than that at l o = 7 μm.This can be explained by the fact that the liquid kinetic energy loss by shock wave emission is relatively small at l o = 4 μm, which is consistent with the previous statements that bubble collapse near a wall causes less energy release and larger bubble re-expansion than farther from the wall [43].The effect of l o on tissue deformation compared to bubble motion is more pronounced at l 0 = 4 μm, as seen in Fig. 7(b).The tissue surface reaches the first minimum location of y sc,min1 = − 0.27 μm at t = 0.76 μs and the maximum location of y sc,max = 1.01 μm at t = 1.12 μs, whose magnitudes are about 52% and 92% larger than those for l o = 7 μm, respectively.After bubble collapse, the tissue is rapidly compressed to reach the second minimum location of y sc,min2 = − 1.10 μm and then undergoes long-lasting deformation.This is because the tissue surface is not restored due to the bubble portion trapped inside the tissue, as depicted in Fig. 5.

Influence of shear modulus
Fig. 8 presents the liquid flow fields associated with acoustic cavitation and tissue deformation for two different shear moduli, G = 10 MPa and G = 0.1 MPa.The tissue with G = 10 MPa behaves like a rigid wall and little liquid momentum due to tissue rebound is transmitted to the bubble, as depicted at t = 1.000 μs in Fig. 8(a).The stiff tissue hinders the fluid flow around the lower part of the bubble, causing non-spherical bubble contraction at t = 1.122 μs, which was observed in the pervious numerical computation of near-wall bubble dynamics [43].The fast liquid inflow on the bubble top surface induces rapid downward bubble migration and strong liquid jet formation.The liquid jet tip close to the tissue transfers dense liquid momentum to the adjacent tissue surface, resulting in rapid compressive deformation at t = 1.138 μs.
Compared to the case of G = 10 MPa, the more flexible tissue with G = 0.1 MPa has larger rebound and causes faster contraction of the bottom bubble surface than the top surface, as seen in Fig. 8(b).This results in a violent liquid annular inflow parallel to the tissue surface and pronounced inverted cone-like bubble shape at t = 1.122 μs.The annular flow impinges the bubble sides and produces a mushroom-shaped bubble characterized by depressions of the lateral sides, which can be observed in the previous experimental and numerical researches for bubble-tissue interactions [34,44].At t = 1.127 μs, the bubble collapses into a thin, elongated shape and is driven away from the tissue surface by the asymmetric velocity field along the vertical axis and the large curvature of the bubble bottom.This indicates that the location of bubble collapse and liquid jet formation is influenced by tissue elasticity as well as bubble-tissue distance.
Fig. 9 shows the integrated influences of l o and G on the bubble maximum radius R b,max and tissue deformation while keeping Π s = Π w and ρ s = ρ w .As l o decreases, R b,max decreases while the tissue deformations |y sc,min1 |, |y sc,min2 | and y sc,max increase.The effect of G is more pronounced in the small l o range, where the liquid momentum is transferred more directly.For l o = 4 μm, R b,max decreases monotonically as the shear modulus increases from G = 0.1 MPa.Similarly, during the bubble growth phase, the first compressive tissue deformation |y sc,min1 | gradually decreases with increasing G.However, R b,max and y sc,min1 do not vary for G > 10 MPa due to the large tissue resistance.On the other hand, y sc,max does not monotonically increase or decrease with G and is largest near G = 1 MPa.This can be explained by the competing effects of tissue repulsion due to the elastic restoring force and tissue attraction toward the center of the contracting bubble.The second compressive tissue deformation |y sc,min2 | that occurs after bubble collapse is found to be clearly maximized near G = 10 MPa, which is the condition for strong liquid jet formation and resulting tissue compressive deformation.The deformation |y sc,min2 | decreases rapidly beyond the point of G = 10 MPa.This indicates that tissue shear modulus is a key parameter for tissue deformation large enough to cause perforation.) than the ambient water, corresponding to a negative reflection coefficient of (I s − I w )/(I s + I w ) = − 0.17.Since the incident and reflected waves are out of phase, the superimposed liquid pressure is less negative around the bubble, which causes the bubble to reach a smaller maximum size at t = 0.870 μs compared to the case of Π s = Π w .On the other hand, when the tissue has a larger acoustic impedance of Π s = 2Π w , the superimposed liquid pressure is more negative around the bubble because the reflected wave is in phase with the incident wave.This results in faster bubble growth and larger bubble size at t = 0.930 μs.

Effect of tissue bulk modulus
The effect of Π s on the temporal behavior of R b and y sc is plotted in Fig. 11.For Π s = 2Π w , while the ultrasonic negative wave overlaps with the reflected wave around the tissue surface, the bubble grows rapidly and reaches R b,max = 5.04 μm at t = 0.930 μs, which is 16.7% larger than for Π s = Π w , as seen in Fig. 11(a).After the bubble expansion stage, the bubble begins to shrink and collapses rapidly, and the ratios of the bubble collapse time Δt col from its maximum radius compared to that for Π s = Π w are estimated as 0.92 and 1.06 for Π s = 0.5Π w and 2Π w , respectively.This can be qualitatively comparable to the ratios, 0.94 and 1.04, of the Rayleigh collapse time [43] proportional to R b,max / ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ |p s,max − p ∞ |/ρ w √ , with slight deviations due to the presence of surrounding tissue.Here, p s,max is the pressure at (r, y) = (R d,e , 0) when the maximum positive pulse reaches to the tissue surface.The changes in tissue bulk modulus also affect the temporal behavior of the tissue.During the initial bubble expansion and contraction period, the first compressive and expansive tissue deformations y sc,min1 and y sc,max become more pronounced as Π s increases.For Π s = 2Π w , the tissue deformed to y sc,min1 = − 0.39μm rapidly expands, reaching y sc,max = 1.18 μm, which is 47.5% larger than that for Π s = Π w .After bubble collapse, the tissue is rapidly compressed and perforated, arriving at y sc,min2 = − 1.42 μm.This indicates that acoustic cavitation near a tissue with larger tissue bulk modulus, leading very asymmetric velocity fields, is more likely to cause larger tissue deformation with perforation.9.This is because the intensity of liquid jet formed by the non-spherical and rapid contraction of the bubble becomes stronger as R b,max increases, which is similarly shown in the numerical observation that the jet velocity is proportional to the bubble expansion size ratio [45].

Effect of tissue density
Fig. 13 presents the effect of tissue density ρ s on acoustic cavitation bubble growth and tissue deformation at l o = 5 μm, G = 1 MPa and Π s = Π w .During the first period of bubble growth and collapse, the temporal variations of R b for ρ s = 2ρ w , ρ w and 0.5ρ w in Fig. 13(a) are similar to those for Π s = 2Π w , Π w and 0.5Π w in Fig. 10(a), respectively.This can be expected form the fact that the acoustic impedance I s = ̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ (ργΠ) s √ is the same for each case.However, the acoustic cavitation bubble-induced tissue deformation pattern is affected differently by ρ s compared to Π s , as seen in Fig. 13(b).As the tissue inertia increases with ρ s , the tissue responds more slowly to the surrounding bubble motion.The tissue with ρ s = 2ρ w is compressed to reach |y sc,min1 | = 0.29 μm, which is 26% smaller than for Π s = 2Π w .Thereafter, the tissue surface rises due to tissue elasticity and sink flow toward the bubble center and reaches y sc,max = 0.86 μm, which is 27% lower than for Π s = 2Π w .However, during bubble collapse and re-expansion, the tissue undergoes violent compressive deformation by the strong liquid jet to |y sc,min2 | = 3.53 μm, which is 2.48 times lager than for Π s = 2Π w .This implies that tissue density is a more important factor in determining strong liquid jet formation and large tissue deformation than tissue bulk modulus.
The combined effects of ρ s and G on R b,max and tissue deformation at l o = 5μm and Π s = Π w are plotted in Fig. 14.As the acoustic impedance of the tissue increases with ρ s , R b,max increases over the entire range of G. Its dependence on ρ s and G is very similar to that on Π s and G.For ρ s = 2ρ w , during the initial bubble growth period, the tissue deforms less due to its significant inertia, so |y sc,min1 | decreases less from that for ρ s = ρ w than for Π s = 2Π w .After the first tissue compression period, the expansive tissue deformation y sc,max increases with ρ s over a wide range of G.However, in the less elastic range near G = 0.1MPa, y sc,max is 0.42 μm for ρ s = 2ρ w , which is smaller than the other cases for ρ s = ρ w and 0.5ρ w .The shear modulus that maximizes y sc,max is close to G = 1MPa and the shear modulus point that maximizes the compressive tissue deformation |y sc,min2 | at ρ s = ρ w is close to G = 10MPa, as when varying Π s .However, when the tissue density increases to ρ s = 2ρ w , | y sc,min2 | has large values in a wide range of G, except near G = 100 MPa, where the liquid jet cannot penetrate the stiff tissue.This indicates that the tissue density is also an important parameter determining strong liquid jet formation and large tissue deformation with perforation.

Conclusions
Acoustic cavitation of a pre-existing microbubble near a viscoelastic tissue was investigated by extending our level-set method for bubble dynamics under ultrasonic pulse conditions to include the viscoelastic and compressible properties of the nearby tissue.Our numerical model was validated through comparison with the predictions from the Keller-Miksis equation including the liquid compressibility and ultrasonic conditions of plane pressure wave.We demonstrated various bubbletissue interactions, including inverted cone-shape bubbles, bubble migration, liquid jet formation, compressive and expansive tissue deformation, and tissue perforation.The bubble was observed to grow larger with increasing tissue bulk modulus and density, and it decreases as the initial bubble-tissue distance decreases below a certain point.The maximum expansive and compressive tissue deformation generally increases with decreasing initial bubble-tissue distance and with increasing tissue bulk modulus and density, except in the case of highdensity, weakly elastic tissues (ρ s = 2ρ w and G = 0.1MPa).The tissue shear modulus conditions that maximize expansive and compressive tissue deformation are close to G = 1 MPa and G = 10 MPa, respectively, unless the tissue density is very large (ρ s < 2ρ w ).This occurs due to the competing effects of tissue elasticity and bubble-induced flow.The present numerical simulations demonstrated that initial bubble-

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.J. Park and G. Son

Fig. 3 .
Fig. 3. Acoustic bubble growth without including nearby tissue compared with the Keller-Miksis prediction at A u = 0.3 MPa, f u = 1 MHz.

Fig. 4 .
Fig. 4. Liquid velocity field with direction vectors (left) and pressure field (right) associated with acoustic cavitation and tissue deformation at l o = 5 μm and G = 1 MPa.The gray and white regions represent the tissue and bubble, respectively.

Fig. 1 ,
and choose the physical properties of air, water and tissue as follows: σ = 7.28 × 10 − 2 N/m, μ w = μ s = 10 − 3 Pa⋅s, ρ w,∞ = 998 kg/m 3 , Π w = 3.31 × 10 8 Pa, μ b = 1.8 × 10 − 5 Pa⋅s, and ρ b,∞ = 1.2 kg/m 3 .A spherical microbubble with an initial radius R bo lies on the central axis at y = l o in a cylindrical domain of r⩽R d,e and Y d,s ⩽y⩽Y d,n , where R d,e = 5897R bo , Y d,s = − 5897R bo and Y d,n = 508R bo .In the cylindrical domain, uniform fine grids with h min = 0.1 μm are chosen around the bubble with r < 19.2R bo and |y| < 19.2R bo , whereas nonuniform coarse grids for the other regions.

Fig. 6 .Fig. 7 .
Fig. 6.Liquid velocity field with direction vectors (left) and pressure field (right) associated with acoustic cavitation and tissue deformation at l o = 4 μm and G = 1 MPa.

Fig. 10 .
Fig. 10.Instantaneous liquid pressure field without including bubble motion (two figures on the left) and instantaneous flow fields associated with acoustic cavitation and tissue deformation at l o = 5 μm, G = 1 MPa, ρ s = ρ w and different bulk moduli: (a) Π s = 0.5Π w and (b) Π s = 2Π w .

Fig. 11 .
Fig. 11.Effect of tissue bulk modulus on acoustic cavitation bubble and tissue behavior at l o = 5 μm, G = 1 MPa, and ρ s = ρ w : (a) temporal change of R b and (b) tissue deformation y sc at the central axis (r = 0).

Fig. 10
Fig. 10 describes the effect of tissue bulk modulus (γΠ) s on the liquid pressure field and bubble growth at l o = 5 μm, G = 1 MPa and ρ s = ρ w .It is noted that the acoustic impedance I s = ρ s C s , where C s is the speed of sound through the tissue, is approximated as I s = ̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ (ργΠ) s √ from Eq. (1).The ultrasonic wave traveling through the liquid water reaches the underlying tissue with Π s = 0.5Π w at t = Y d,n /C w = 0.33 μs, as seen in Fig.10(a).It is partially reflected by the tissue, which has a smaller acoustic impedance (I s = ̅̅̅̅̅̅̅ 0.5 √ I w ) than the ambient water, corresponding to a negative reflection coefficient of (I s − I w )/(I s + I w ) = − 0.17.Since the incident and reflected waves are out of phase, the superimposed liquid pressure is less negative around the bubble, which

Fig. 12
shows the combined effects of Π s and G on R b,max and tissue deformation at l o = 5 μm and ρ s = ρ w .As the acoustic impedance of the tissue increases with Π s , R b,max , |y sc,min1 |, |y sc,min2 | and y sc,max all increase over almost the entire range of G. Their dependence on G for various Π s ′s is almost similar to that on G for various l o ′s, except for R b,max in the range of G < 1 MPa, as seen in Fig. 9. Regardless of the value of Π s , the maximum expansive deformation y sc,max of the tissue occurs at G = 1MPa and the maximum compressive deformation |y sc,min2 | occurs at G = 10MPa.The deformation |y sc,min2 | is 4.57 μm at Π s = 2Π w and l o = 5 μm, which is 41% larger than |y sc,min2 | at Π s = Π w and l o = 4 μm in Fig.

Fig. 12 .
Fig. 12. Combined effects of tissue bulk modulus and shear modulus on bubble growth and tissue deformation at l o = 5 μm and ρ s = ρ w .

Fig. 13 .
Fig. 13.Effect of tissue density on acoustic cavitation bubble and tissue deformation at l o = 5 μm ,G = 1 MPa, and Π s = Π w : (a) temporal change of R b and (b) tissue deformation y sc at the central axis (r = 0).

Fig. 14 .
Fig. 14.Combined effects of tissue density and shear modulus on bubble growth and tissue deformation at l o = 5 μm and Π s = Π w .