Numerical study of the effect of liquid compressibility on acoustic droplet vaporization

In acoustic droplet vaporization (ADV), a cavitated bubble grows and collapses depending on the pressure amplitude of the acoustic pulse. During the bubble collapse, the surrounding liquid is compressed to high pressure, and liquid compressibility can have a significant impact on bubble behavior and ADV threshold. In this work, a one-dimensional numerical model considering liquid compressibility is presented for ADV of a volatile microdroplet, extending our previous Rayleigh-Plesset based model [Ultrason. Chem. 71 (2021) 105361]. The numerical results for bubble motion and liquid energy change in ADV show that the liquid compressibility highly inhibits bubble growth during bubble collapse and rebound, especially under high acoustic frequency conditions. The liquid compressibility effect on the ADV threshold is quantified with varying acoustic frequencies and amplitudes.


Introduction
Acoustic droplet vaporization (ADV) is emerging in the field of noninvasive targeted therapy and diagnosis, as reviewed in Refs. [1,2]. Nanodroplets are vaporized into microbubbles when the amplitude of ultrasonic pulses exceeds a certain intensity called ADV threshold. An ultrasound-triggered bubble, which oscillates with alternating the positive and negative pressure portion of the pulses, can rapidly collapse with or without rebounds depending on the acoustic amplitude [3][4][5][6][7]. During the rapid bubble collapse, the liquid region near the bubble is also compressed to a high pressure above the critical pressure, resulting in a shock wave emission [8]. The liquid compressibility can affect the bubble behavior after collapse and thus ADV threshold, as indicated in Refs. [6,7,9]. However, a general numerical model including the liquid compressibility effect on ADV is lacking in the literature.
Numerical studies of ADV were performed using Rayleigh-Plesset (RP) equations [6,7,[10][11][12] and multidimensional numerical models [13][14][15]. Although the multidimensional models can be applied to nonspherical bubble growth in ADV and the interactions of ADV bubbles and walls, they are not efficient for predicting the ADV threshold because very fine computational meshes are required to treat the initial bubble nuclei. Most of numerical models for ADV were developed from onedimensional conservation equations for spherical bubble growth. The RP equation was derived by integrating the radial momentum equations for droplet and surrounding water regions with the interface conditions at the bubble and droplet surfaces. Assuming that the bubble was an ideal gas, the bubble growth and survival due to acoustic pulses were investigated. Recently, Park and Son [7] improved the RP equation for ADV using the van de Waals (VDW) equation of state applicable to the supercritical state of a collapsing and rebounding bubble. However, all studies of ADV mentioned above did not consider liquid compressibility.
The effect of liquid compressibility was included in several studies investigating the bubble growth and collapse in infinite water. Gilmore [16] and Keller-Miksis [17] equations, which are a kind of RP equation applied to a weakly compressible liquid, were widely used to investigate bubble motions in acoustic fields. These equations can be derived from the integration of the radial momentum equation between the bubble surface and the domain boundary introducing a velocity potential, as described in Ref. [17]. The velocity potential cannot be determined directly from the continuity equation as in the incompressible liquid case, but is expressed as a general solution of the wave equation approximated from the continuity and momentum equations. In the compressible RP equations, which are derived by eliminating the time derivative of the velocity potential, the liquid compressibility effect is considered in terms of the sound speed and the time derivative of the boundary pressure, which are not shown in the incompressible RP equation. Although the compressible RP equations were applied to bubble growth and collapse in infinite water, they were not applied to the ADV case, in which the compressible RP equations cannot be easily derived for the droplet and water regions without including unknown variables such as the time derivatives of the droplet surface pressures or the velocity potentials in the two regions. As an alternative to the compressible RP-based model, numerical models directly solving the conservation equation without using the weak liquid compressibility assumption were proposed in several studies [18][19][20][21] for bubble motion in a compressible liquid. They investigated the bubble behaviors triggered by acoustic pressure waves and the pressure wave transition during bubble collapse. However, the methods were not extended to ADV which includes the droplet and water regions.
In this work, a numerical method for bubble growth in ADV considering liquid compressibility is developed by extending our previous RP based model [7]. One-dimensional spherical conservation equations of mass, momentum and energy are directly solved for bubble motion in ADV. The present model is applied to the ADV experimental cases of Aliabouzar et al. [22].     compressible liquid (droplet and ambient water) regions. A spherical dodecafluoropentane (DDFP) microdroplet with a radius of R d and a boiling point of 29 • C at 1 atm is placed in the ambient water. An acoustic pressure wave of p a is generated at the domain boundary of r = L, and transmitted through the droplet at r = 0. A DDFP bubble with a radius of R b is assumed to grow at the droplet center.

Governing equations
Assuming that the flows are spherically symmetric, the conservation equations for compressible liquid regions are expressed as where the subscript l denotes the droplet region (d) if R b < r < R d and the ambient water region (w) if R d < r < L. In Eq. (1), consisting of the unsteady and mass flux terms, the liquid velocity cannot be easily evaluated as a function of the bubble radius R b due to the liquid compressibility condition (Dρ l /Dt ∕ = 0). Eq. (2) contains not only the inertia and pressure terms, but also two terms of viscous stress, which do not appear in the incompressible fluid case. Each term in Eq. (3) represents internal energy change, pressure work, heat transfer and viscous dissipation, respectively. The equations for each liquid region are combined by the following conditions at r = R d : Assuming uniform pressure and temperature distributions inside the bubble [6,[10][11][12], the boundary conditions at r = R b are expressed as In Eq. (8), u d is the droplet-side velocity at the bubble surface, which can be obtained from Eq. (2), and is not equal to the bubble surface velocity dR b /dt due to the phase-change rate G/ρ d . The phase-change mass flux G at r = R b can be written as (11) The boundary conditions at r = L are specified as where p a,L is the acoustic pressure imposed at the domain boundary.
Assuming that the bubble density is uniform inside the bubble, ρ b is determined by the bubble mass conservation, The VDW Eq. (15) is used to consider the bubble compressibility and supercritical state of a collapsing and rebounding bubble. It is combined with Clausius-Clapeyron Eq. (16) and ρ b from Eq. (14) to determine the bubble temperature T b and pressure p b , as described in our previous study [7].
where the constants a and b are obtained from the critical properties [23]. The bubble pressure p b is used as the boundary condition to determine the droplet-side pressure at the bubble surface, as given by Eq. (9). When the bubble pressure exceeds the critical value p c , we assume that G becomes zero and the bubble is isentropic during the supercritical state as done by Refs. [7,24].
The liquid density ρ l is obtained from Eq. (1) and the liquid pressure p l is determined from the Tait equation, where the constants are evaluated as Γ w = 331 MPa and n w = 7.15 for water [25]. For the DDFP liquid, Γ d = 36.4 MPa and n d = 7.23 are determined by fitting the NIST data [26] at 303 K as depicted in Fig. 2 and referring to the sound speed of C d = 406 m/s reported in Ref. [27]. Introducing the moving coordinate ξ = r − R b and using the level-set method [5,14,21] so that the conservation equations implicitly satisfies the interface conditions at r = R d , given by Eqs. (4)-(7), the conservation Eqs. (1)-(3) can be rewritten for the liquid regions as (ρc) l (

∂T l ∂t
where the relative velocity v l and the step function α are defined as For a two-phase mixture cell that includes the droplet surface (r = R d ), the properties ρ l , (ρc) l , μ l and λ l are evaluated by interpolation of the droplet and water properties using the level-set method [5,14,21]. We discretize the conservation equations in a staggered grid system using the second-order ENO and central difference schemes for spatial discretization of the convection and diffusion terms, and a third-order TVD Runge-Kutta method [28] with an adaptive time step for temporal discretization of the conservation equations and other unsteady equations.

Computational conditions
Computations of ADV are performed for a DDFP microdroplet immersed in the ambient water at p ∞ = 1 atm and T ∞ = 310 K. We choose the initial droplet radius of R do = 0.47 μm from the ADV experimental condition of Aliabouzar et al. [22] and a large computational domain of L = 2.74 mm to prevent reflection of acoustic waves at the boundary (r = L). The properties of DDFP and water are obtained from our previous work [7]. The acoustic pressure pulse imposed at the domain boundary (r = L) is expressed as Here, P a,L and frequency f a vary in this work keeping N a = 8 based on Ref. [22]. For a computational domain of R b < r < L (or 0 < ξ <L − R b ), uniform fine meshes of size Δξ o are chosen near the droplet (0 < ξ < 20R do ), nonuniform meshes with sizes increasing by a factor of 1.0004 in the middle region (20R do < ξ < 1600R do ), and uniform meshes in the outer region (1600R do < ξ < L − R b ).

Acoustic pressure amplification in compressible liquids
Computations are first performed on acoustic pressure waves traveling in the ambient water without a droplet to validate the present model for compressible flows. The domain size chosen with L = 2.74 mm is 4 times the wavelength (C w,∞ /f a ) for f a = 2.25 MHz and almost 18 times the wavelength for f a = 10 MHz. Fig. 4 shows the evolution of pressure waves generated under the acoustic conditions of P a,L = 10 kPa and f a = 2.25 MHz or 10 MHz at r = L. The pressure waves travel through the ambient water with C w,∞ = 1540 m/s and reach the domain center (r = 0) at t = 1.775 μs. Assuming that the water viscosity is negligible and the pressure amplitude is small, the analytic solution of the linearized pressure wave equation for a onedimensional spherical case can be derived as [21] p(t, r) = p ∞ + P a,L L r where The pressure at r = 0 is evaluated as where The present numerical results have no difference from the analytical predictions given by Eq. (26). As the acoustic wave approaches r = 0, its magnitude is amplified to 50P a,L for f a = 2.25 MHz and 220P a,L for f a = 10 MHz. Since the acoustic amplitude varies with the radial location in a compressible liquid, we choose the reference amplitude from the amplitude P a,0 at r = 0 in the absence of a droplet. In a typical experimental setup for ADV, the acoustic amplitude was measured using a hydrophone at the underwater focal spot of an ultrasonic transducer, and then a droplet was placed at the focal spot [29,30]. This selection of acoustic amplitude as P a,0 seems to be consistent with the hydrophone measurement method. It is noted that a large number of mesh points of N ξ = 25, 000 are used in the present computations to avoid numerical oscillations occurring at high acoustic frequencies, as depicted in Fig. 4d for f a = 10 MHz.
Calculations are next made for an acoustic pressure wave traveling through a DDFP droplet immersed in the ambient water without bubble nucleation. The results are plotted in Figs. 5 and 6. Although the presence of a droplet affects the pressure wave passing the droplet region during the first cycle, as seen in Fig. 5, its influence on the pressure p 0 at r = 0, except for the first cycle, is insignificant at f a = 2.25 MHz regardless of the droplet radius and acoustic amplitude, as depicted in Fig. 6a and 6b. The droplet effect on the pressure amplification is pronounced for f a = 10 MHz and R do = 10 μm, as seen in Fig. 6c. This is expected from Eq. (29), which indicates the pressure amplification is proportional to f a and inversely proportional to the sound speed. Since the sound speed of DDFP droplet is lower than that of ambient water, the pressure amplification in the droplet increases as the droplet becomes larger, especially for the high acoustic frequency case of f a = 10 MHz. This is consistent with the observation in Ref. [27] that the amplification of pressure wave is proportional to f a with increasing the droplet size.

Acoustic vaporization in compressible liquids
We now consider acoustic vaporization of a DDFP droplet with R do = 0.47 μm in the ambient water at T ∞ = 310 K and p ∞ = 1 atm. An acoustic pressure wave generated at the domain boundary (r = L) with P a,L = 10 kPa and f a = 2.25 MHz is amplified to P a,0 = 0.5 MPa when it reaches the domain center (r = 0), as described in Section 3.1. A bubble with an initial radius of R bo is assumed to form at the droplet center. We choose R bo = 80 nm, which was also used in Refs. [7,12]. Since the effect of R bo on ADV was investigated in our previous study [7], we focus on the effect of liquid compressibility on ADV in this work.
Considering the Laplace pressure and acoustic pulsing condition, the initial bubble pressure is evaluated as Fig. 7. ADV at f a = 2.25 MHz and P a,0 = 0.5 MPa: temporal variations of (a) bubble radius, pressure, (b) temperature and mass.
where p a,0 is the acoustic pressure at the droplet center, as depicted in Fig. 6. Assuming that bubble nucleation occurs in thermal equilibrium with the ambient condition, p bo = p sat (T ∞ ), and using the properties of is reduced below T ∞ , the evaporation at r = R b increases the bubble mass m b . Thereafter, during the positive pulsing period, R b decreases to 0.13 μm. While T b increases above T ∞ , the condensation at r = R b reduces m b . As the bubble further compresses or collapses due to the positive pulsing, p b exceeds the critical pressure p c = 2.25 MPa and then rapidly rises to a very high pressure at the supercritical state without phase change. When the liquid pressure increases near the collapsing bubble and balances with the liquid inertia compressing the bubble, the bubble starts to bounce back. After subsequent collapses and rebounds, R b and m b finally become near zero, indicating an irreversible collapse of the bubble.
The liquid pressure and temperature profiles are plotted in Fig. 8. While the bubble compresses due to the positive pulsing, its pressure increases until the bubble bounces back at t = 2.21 μs (Fig. 8a). The bubble rebound generates a shock wave at t = 2.23 μs. The shock wave travels outward and attenuates to 37MPa, 17 MPa and 3.3 MPa as it reaches r = 1 μm, r = 2 μm and r = 10 μm, respectively. This attenuation is nearly proportional to 1/r as described in Ref. [31]. Another shock wave is emitted at t = 2.28 μs after the second bubble rebounds.
The liquid pressure changes in a large portion of the liquid region, but the liquid temperature variation is limited in a narrow region near the bubble, as depicted in Fig. 8b.

Effect of liquid compressibility on ADV
In Fig. 9, we compare the ADV results for compressible and incompressible liquids. Considering that the pressure wave in the incompressible liquid case immediately reaches the droplet due to its infinite sound speed, and its magnitude is not amplified as P a,0 = P a,L , the results are plotted against the relative time t * to have the same initial bubble growth pattern for the compressible and incompressible liquid cases. The discrepancy in both results is significant after the first bubble collapse and rebound. Compared to the results for the incompressible liquid case, the rebounding bubble has the first local maximum radius that is about 44.6% smaller for the compressible liquid case, the bubble surface speed becomes lower, and the time interval between two subsequent rebounds is reduced, which was also reported in Ref. [9]. Since the bubble radius is smaller during the collapse and rebound periods, T b remains higher in the compressible liquid case (Fig. 9c) and the bubble condenses and eventually disappears, whereas the bubble completely vaporizes at t * = 0.68 μs under the incompressible liquid condition. Fig. 10 shows the liquid pressure, velocity and ρ l T l profiles near r = R b during the first bubble collapse and rebound period for the compressible liquid case. The profiles rapidly vary in a narrow liquid during the bubble collapse period, whereas their variations spread differently in a wider region with a shock wave traveling outward during the bubble rebound period. For the incompressible liquid case with no shock, however, the liquid pressure, velocity and ρ l T l similarly change in reverse order during the bubble collapse and rebound periods, as plotted in Fig. 11. To further clarify the effect of liquid compressibility on ADV, we analyze the temporal changes in the liquid energy while the bubble temperature T b is exceeds the critical temperature T c during the first collapse and rebound period of t 1 ⩽t⩽t 2 , where the collapse time t 1 = 0.3226 μs and the rebound time t 2 = 0.3249 μs for the compressible liquid case, and t 1 = 0.3172 μs and t 2 = 0.3175 μs for the incompressible liquid case. Neglecting the contributions of heat transfer and viscous work to the liquid energy change during the collapse and rebound period at the supercritical state, the conservation equation of liquid energy can be written as ∂r 2 ρ l E l ∂t + ∂r 2 u l ρ l E l ∂r = − ∂r 2 p l u l ∂r (31) where E l = c l (T l − T ∞ ) + u 2 l /2. The integration of Eq. (31) from r = R b to r = L ′ , which is chosen as 40R do or 20 μm, yields d dt Considering that the second term on the right hand side of the above equation is relatively negligible because u l = dR b /dt at r = R b and E l ≃ 0 at r = L ′ , we have where the liquid kinetic energy (KE), the internal energy (IE) and the flow work (FW) are defined as The results are plotted in Fig. 12. The loss of kinetic energy for the compressible liquid case during the bubble collapse and rebound is almost 60%, which is 6 times greater than for the incompressible liquid case. The loss of kinetic energy results in the amount of internal or thermal energy increase for both cases. For the compressible liquid case, the liquid internal energy is increases by 70% due to not only the bubble motion and shock wave but also the flow work of positive pressure pulsing during the bubble collapse and rebound period. Fig. 13 shows the ADV results at a higher frequency of f a = 10 MHz and P a,0 = 1.57 MPa. Compared to the case of f a = 2.25 MHz, the bubble radius varies more dynamically because the sign of pressure pulse changes rapidly with increasing f a . When the loss of liquid kinetic energy during the first collapse and rebound period is evaluated as done in Section 3.3, it is 4.5% for the incompressible liquid case and increases to 60% for the compressible liquid case. The liquid internal energy for the compressible liquid case is increased by 43%, which is caused by the flow work of negative pressure pulsing in the bubble collapse and rebound period. The bubble vaporizes completely at t * = 0.16 μs for the incompressible liquid case, but the bubble under the compressible liquid condition undergoes multiple collapse and rebounds and then completely vaporizes 0.82 μs later than in the impressible liquid case. Fig. 14a and b present the first local and global maximum bubble radii in ADV at f a = 10 MHz. For the incompressible liquid case, the first local maximum bubble radius increases gradually with P a,0 , but the  global maximum bubble radius increases abruptly near the ADV threshold of P a,0 = 0.95 MPa. The rapid bubble growth after collapse and rebound helps the bubble to completely vaporize, as described in Refs. [6,7]. On the other hand, the global maximum bubble radius for the compressible liquid case gradually increases near the ADV threshold of P a,0 = 1.57 MPa, as seen in Fig. 14b, because the liquid compressibility inhibits the bubble growth during the bubble collapse and rebound, as depicted in Fig. 13.

Effects of liquid compressibility and acoustic frequency
The combined effects of liquid compressibility and f a on the ADV threshold are shown in Fig. 14c. The liquid compressibility is observed to increase the ADV threshold, especially as f a increases. The predicted ADV thresholds considering liquid compressibility are 0.52 MPa for f a = 2.25 MHz and 1.57 MPa for f a = 10 MHz, which are 6% and 65% higher than for the incompressible liquid cases, respectively. The numerical predictions taking account liquid compressibility are more comparable to the experimental data except near f a = 15 MHz. The experimental data at f a = 15 MHz is close to the computed ADV threshold at a slightly lower frequency of f a = 13.5 MHz. This discrepancy may be due to the    effect of phospholipid coating or encapsulation of the droplet, as investigated in the previous study [7], or the nonlinear effect of acoustic wave transmission occurring under high-frequency conditions, as described in Refs. [27,32].

Conclusions
The effect of liquid compressibility on ADV of volatile microdroplets was numerically investigated by directly solving the conservation equations for the liquid droplet and surrounding water regions. The computations of acoustic pressure waves traveling in the compressible liquid regions showed that the magnitude of acoustic waves near the domain center is amplified in proportion to the acoustic frequency and inversely to the sound speed. The numerical results for ADV demonstrated that the bubble motion is significantly influenced by liquid compressibility during the bubble collapse and rebound periods. As the bubble compresses, the bubble pressure exceeds the critical pressure and rapidly rises to a very high pressure at the supercritical state without phase change. During the subsequent rebound period, the bubble generates a shock wave and the liquid region near the bubble loses about 60% of its kinetic energy under the compressible liquid condition, which is 6 times greater than under the incompressible liquid condition. This indicates that the liquid compressibility highly inhibits bubble growth during bubble collapse and rebound. As a result, the ADV threshold is higher for the compressible liquid case than for the incompressible liquid case, especially under high-frequency conditions. The numerical predictions of ADV thresholds are in better agreement with the experimental data, including the liquid compressibility effect.

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.