Numerical investigation of acoustic vaporization threshold of microdroplets

Highlights • A numerical method is presented for the acoustic vaporization threshold of microdroplets by improving the Rayleigh-Plesset equation to properly treat the supercritical state occurring at rapid bubble collapse.• The van der Waals equation of state is employed to more accurately consider the supercritical state instead of the ideal gas equation.• The effective ADV threshold is observed to increase as the acoustic frequency increases and the droplet radius decreases, as experimentally observed in the previous studies.• The present numerical predictions of ADV threshold comparable to the experimental data.


Introduction
Acoustic vaporization of volatile submicron droplets has received attention due to its medical potentials [1]. Due to the submicron size and stability of droplets, they can reach cancerous tissues through the endothelial cell wall [2] and vaporize into microbubbles that can be used as ultrasound contrast agents. The acoustic droplet vaporization (ADV) enables the targeted drug delivery and precise ultrasound imaging [3]. The large volume change and strong flow due to vaporization can be used for various medical therapies such as a gas embolotherapy and thermal ablation [4]. Despite of these potentials for medical applications, the use of ADV in the human body is challenging. The ADV requires a strong acoustic pulse beyond the critical amplitude, called the ADV threshold, to vaporize droplets, which can induce negative biological effects [5,6].
Many experimental studies have been conducted for clarifying the ADV mechanism. Based on the previous experimental data for the acoustic threshold, Aliabouzer et al. [7] observed that the ADV threshold pressure decreases as the droplet size, ambient temperature, and number of acoustic cycles increase [8][9][10]. However, the effect of acoustic frequency on the ADV threshold was not consistent in the literature. Some studies showed that the ADV threshold increases with acoustic frequency [11,12] whereas others showed the opposite trend [13,14]. While analyzing the experimental observation that the acoustic threshold decreases as the acoustic frequency increases, Shpak et al. [15] noted the superharmonic effect that a spherical droplet concentrates an incident acoustic wave at its inner point. Their theoretical analysis showed that the superharmonic focusing effect, triggering the generation of bubbles, increases as the acoustic frequency and the droplet size increase. However, the superharmonic effect was observed to disappear in small droplets of 2 μm or less.
The ADV thresholds reported in the literature differ among researchers. It can be attributed not only to differences in threshold measurement techniques, such as hydrophone and imaging-based methods, but also to differences in experimental setups containing droplets, such as gel phantoms and tubes. The tube setups can cause uncertainties in threshold measurements through the scattering and reflection of acoustic waves. Aliabouzer et al. [7,16] conducted experiments for the ADV threshold of microdroplets in a tubeless setup to minimize the effects of tube wall. Their results showed that the ADV threshold for dodecafluoropentane (DDFP) droplets increases with the acoustic frequency.
Theoretical studies have also been carried out for predicting the ADV threshold. Vlaisavljevich et al. [12] applied classical nucleation theory to prediction of the ADV threshold for perfluorohexane (PFH) nanodroplets used in histrotripsy. Miles et al. [17] extended the classical nucleation theory to the ADV threshold for DDFP microdroplets by combining with the superharmonic effect. Their predicted threshold pressure is about 9.3 MPa for a droplet of 10 μm radius, which is much higher than the experimental data of Aliabouzer et al. [16]. As another predictive model for the ADV threshold, one-dimensional numerical models for bubble growth have been applied to determine the conditions for the survival of a small bubble assumed to nucleate. The Rayleigh-Plesset (RP) equation for spherical bubble growth was extended for ADV by Qamar et al. [18] considering a tube geometry, Shpak et al. [19] including the effect of thermal boundary layer, Pitt et al. [20] investigating the effects of acoustic parameters, Doinikov et al. [21] computing the vapor generation rate from the energy equation, and Guédra and Coulouvrat [22] and Lacour et al. [23] including the effect of droplet encapsulation. The RP based numerical model has the advantage over the multidimensional numerical models for ADV developed in Refs. [24][25][26] that it can use a very fine mesh needed for the calculation of the initial bubble nucleus. Recently, Lacour et al. [27] applied the RP based numerical model to investigate the ultimate fate of a DDFP bubble nucleus forming in a DDFP droplet and determine the ADV threshold in various acoustic frequencies. The bubble growth pattern is observed to have three regimes: irreversible collapse without growth (regime I), direct vaporization without rebound (III) and intermediate behavior with collapses, rebounds and oscillations (regime II). They provided a phase diagram for the three regimes of bubble growth in ADV, which is useful to understand the complex bubble/droplet behaviors in ADV. However, their model needs to be modified in that the bubble has a very high temperature, to the order of a few MKs, while it collapses with oscillations in the intermediate regime III, and the ADV threshold is much higher than the experimental data of Aliabouzer et al. [16].
In this work, the RP based numerical model for predicting the ADV threshold is improved by properly treating the supercritical state occurring at rapid bubble collapse. The van der Waals (VDW) equation is employed to more accurately consider the supercritical state instead of the ideal gas equation, which was used in most of the previous studies referenced above for ADV. The present model is tested on the previous experimental data obtained in a tubeless setup [16]. The effects of fluid properties and droplet radius on the bubble growth pattern and the ADV threshold are investigated. Fig. 1 shows the one-dimensional (spherical) schematic of ADV used in the present work. A spherical DDFP droplet with R d is immersed in ambient water, and a vapor bubble with R b is assumed to nucleate at the center of the droplet when the acoustic pressure pulse p a is excited from the boundary of the domain. The DDFP droplet and ambient water are considered incompressible fluids, and the DDFP bubble is considered a VDW gas rather than an ideal gas to more accurately consider the supercritical state occurring at bubble collapse. The VDW equation of state is expressed as

Mathematical modeling
where R g is the gas constant and the constants a and b are evaluated as [28] a = 27 64 Here, the subscript c denotes the critical state.

Governing equations
Assuming that the vapor pressure, density and temperature are uniform, the one-dimensional conservation equations in the liquid (droplet and ambient water) regions as well as the mass equation in the bubble region are expressed as    (6) where the subscripts b and l indicate bubble and liquid regions, respectively. The conservation equations are solved using the matching conditions at the bubble surface (r = R b ) with phase change,

∂T l ∂t
∂T d ∂r (10) and the matching conditions at the droplet surface (r = R d ) without phase change, Here, the subscripts d and w denote droplet and ambient water regions, respectively, G is the phase-change mass flux, Ṙ b is the bubble surface velocity, and U b and U d are the vapor and liquid velocities at r = R b , respectively. The boundary condition at r = ∞ is specified as where p ∞ and p a are ambient and acoustic pressures, respectively. The vapor and liquid velocity profiles are solved from Eqs. (3) and (4) with the boundary conditions as Assuming ρ d ≫ρ b , integration of the momentum equations over the droplet and ambient water regions with the boundary conditions result in the Rayleigh-Plesset equation [21,26], where A fourth-order Runge-Kutta (RK) method with an adaptive computational time step is used to solve U d from Eq. (19), which includes G, whose calculation procedure will be described later, R b to be determined from Eq. (7), ρ b from Eq. (23) for the bubble mass conservation, R d from Eq. (24) for the droplet mass conservation, p b and T b from the VDW Eq.
where T b,∞ is the saturation temperature at p ∞ . It is noted that Eq. (23) has an advantage in preserving the bubble mass over Eq. (26) used in many previous studies [21][22][23]27], which can be derived from Eq. (23) and the thermodynamic relation To calculate the phase-change mass flux G from Eq. (10), the liquid-side

∂T l ∂t
where the subscript l represents the droplet or ambient water region when the computational cell is a single-phase cell that does not include the droplet surface (r = R d ). For a droplet-water mixture cell, (ρc) l , μ l and λ l are interpolated from the droplet and water properties, as done in Ref. [25] using a level-set method. Eq. (27) is spatially discretized using second-order essentially non-oscillatory (ENO) and central difference schemes for the convection and diffusion terms, respectively. A fourthorder RK method with an adaptive time step is applied to solve the transient differential equation. While a bubble collapses rapidly in ADV, its pressure (or temperature) can exceed a critical value. To simulate the supercritical process (p b ⩾p c or T b ⩾T c ), where liquid and vapor phases are not distinguished, we assume that the phase-change mass flux is zero (G = 0), as done by Akhatov et al. [29] for computing a laser-induced bubble, and the bubble state varies in an isentropic process, which is expressed as [28] c Integration of Eq. (28) using the VDW Eq. (1) yields The bubble temperature is determined from Eq. (29) if ρ b ⩾ρ c and from the combined Eqs. (1) and (25) if ρ b < ρ c .

Computational conditions
We consider acoustic vaporization of a DDFP droplet in water at p ∞ =1 atm and T ∞ = 310 K, as depicted in For computation of ADV, we use the following acoustic pressure pulse.
Here, A a and f a vary in this work keeping the number of acoustic cycles as N a = 8, referring to the experimental condition of Aliabouzer et al. [16]. The minus sign is chosen considering that bubble nucleation occurs in the negative pressure period, as experimentally observed in Ref. [31].
As the initial conditions, we choose R do = 0.47 μm as the droplet radius and T l = T ∞ = 310 K as the liquid temperature profile, based on Ref. [16]. The initial bubble radius R bo is unknown and can be considered an external parameter. In this work, we choose R bo = 80 nm as a base case, which was also used by Lacour et al. [23,27] while predicting the ADV threshold of DDFP droplets. Introducing the time t bo at bubble nucleation after the acoustic pulsing starts and considering the Laplace pressure caused by the surface tensions at r = R bo and r = R do , the initial bubble pressure p bo can be estimated as The bubble nucleation time t bo is another parameter affecting the bubble growth in ADV. In the case of t bo = 0, which means that bubble nucleation occurs when the acoustic pulsing starts, we obtain p bo = 0.552 MPa from Eq. (31) and the corresponding bubble temperature of T bo = 363 K from Eq. (25). This case has a mismatch of T l ∕ = T b near the bubble surface (r = R b ) and thus the initial droplet temperature profile needs to be modified as done in the previous studies [19,22] introducing a thin liquid layer for smooth transition from T b to T ∞ . As another case of t bo , we assume that a bubble nucleates when its nucleus is in thermal equilibrium with the surrounding liquid at T ∞ . Assuming T bo = T ∞ = 310 K, the corresponding bubble pressure is estimated as p bo = 0.131 MPa from Eq. (25). The acoustic amplitude required for bubble nucleation is A a > 0.42 MPa and given A a , the bubble generation time t bo is determined from Eq. (31).
The liquid-side energy Eq. (27) formulated in a moving coordinate is solved in a computational domain of 0 < ξ < L (or R b < r < L − R b ). We use uniform fine meshes of size Δr in a region of 0 < ξ < 15R do and nonuniform meshes in the outer region of 15R do < ξ < 10 6 R do .
Time and space resolutions are tested to determine the proper time step Δt and mesh size Δr, as plotted in Fig. 2. The bubble growth computed with different time steps has no difference while using adaptive time steps with a maximum of Δt⩽0.2 ns. In Fig. 2b, the results obtained with different mesh sizes show that the bubble oscillations converge as the mesh size decreases below Δr = 10 nm. Therefore, we use Δr = 10 nm and the adaptive time step in the present computations.  Fig. 3a, when the bubble temperature is reduced to T ∞ = 310 K by the acoustic pressure. The bubble grows during a low pressure pulsing period of t bo < t < 0.165 μs and reaches its maximum radius of R b = 0.13 μm. As the bubble pressure p b decreases with bubble growth and the corresponding T b decreases below T ∞ , the bubble mass m b increases, as seen in Fig. 3c, because of the evaporative heat transfer from the droplet bulk. Thereafter, while the acoustic pressure increases, the bubble shrinks rapidly and p b increases. As T b becomes higher than T ∞ , m b decreases because of to the condensing heat transfer from the bubble surface. The computation is carried out until the bubble and its mass become very small to

Results and discussion
bo . This result means that the bubble nucleus collapses with no rebound at A a = 0.43 MPa. The irreversible bubble collapse occurs in the acoustic amplitude range of A a ⩽0.435 MPa, which is referred to bubble growth regime 1 in the present work. In this regime, the droplet size hardly changes with acoustic pulsing, which was also observed in the previous experimental studies [8,32] using high speed microscopy to investigate the ADV behavior.
When the acoustic amplitude increases to 0.44 MPa⩽A a ⩽0.48 MPa, a bubble nucleates earlier and grows larger during the first negative pulsing period (t < 0.5/f a = 0.22 μs). However, while the positive pulsing is applied (0.22 μs < t < 0.44 μs), the bubble collapses and rebounds, which is caused by the unbalance between liquid inertia and bubble pressure. The number of bubble collapse and rebound increases with A a . In Fig. 3d, while the bubble collapses due to the inertia of the surrounding droplet and ambient water, its temperature exceeds the critical temperature of T c = 430 K, which is evaluated from Eq. (25) with p c = 2.045 MPa for DDFP. In this work, we assume G = 0 at the supercritical state. As the bubble collapse builds up the bubble pressure and also the liquid pressure, the bubble rebounds rapidly. During the positive pulsing period, the bubble at the subcritical state (T b < T c ) has a higher temperature than T ∞ and thus m b decreases to zero. The acoustic amplitude range of 0.44 MPa⩽A a ⩽0.48 MPa, where the irreversible bubble collapse occurs after multiple collapses and rebounds, is referred to bubble growth regime 2. This regime corresponds to the experimental observation of Sheeran et al. [32] that some droplets expand during the negative cycle of acoustic pulsing and then disappear during the positive pulsing period.
When the acoustic amplitude further increases to 0.49 MPa⩽A a ⩽0.59 MPa, the negative pulsing causes the bubble to grow enough to survive during the following positive pulsing period, as seen in Fig. 4b. It is noted that whenever the bubble collapses and rebounds, its temperature exceeds the critical temperature (Fig. 4d). and thus proper modeling of the supercritical process is essentially important for simulating bubble behavior after collapse and rebound. When the second negative pulsing period begins at t = 0.44 μs, the bubble grows again and then reaches its maximum radius of R b,f = 2.26 μm at the time of complete vaporization. The acoustic amplitude range of 0.49 MPa⩽A a ⩽0.59 MPa, where the droplet vaporization is complete after bubble collapses and rebounds, is referred to bubble growth regime 3. This bubble growth behavior was observed in several experimental studies [8,[31][32][33][34]. The minimum acoustic pressure required for complete droplet vaporization is A a = 0.49 MPa. As A a increases above 0.6 MPa, the droplet vaporization is directly complete during the first negative pulsing period, which is called bubble growth regime 4. The minimum pressure for direct droplet vaporization without any collapse and rebound is called the direct threshold in Refs. [23,27]. This bubble growth regime was observed by Li et al. [34] while visualizing the formation of toroidal bubbles from high acoustic pressure. Fig. 5 presents the first local maximum bubble radius R b,lmax , which is observed during the first negative pulsing period, and the global maximum bubble radius R b,gmax in the vaporization period. It is seen that R b,lmax increases linearly with A a . The bubble growth regimes are distinguished as follows: regimes 1 and   [16]., respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) The effective ADV threshold A eth is defined in this work as the minimum acoustic amplitude that satisfies R b,gmax (A a ⩾A eth ) = R b,f . The effective ADV threshold is A eth = 0.49 MPa and the direct ADV threshold is A dth = 0.60 MPa in the base case. Fig. 6 shows bubble growth at an acoustic frequency increased to f a = 10 MHz. For A a = 0.49 MPa, which is the effective ADV threshold at f a = 2.25 MHz, a bubble immediately collapses without rebound because the first negative acoustic pulsing period is significantly reduced in this high frequency case. As A a increases, the maximum bubble radius increases and the bubble oscillations including collapses and rebounds become pronounced. For A a = 0.76 MPa, a bubble grows with oscillations during the acoustic pulsing periods of 8/f a and thereafter it shrinks with multiple collapses and rebounds and finally disappears.

Effect of frequency of acoustic pulse
The acoustic amplitude increased to A a = 0.80 MPa causes complete vaporization after multiple bubble oscillations, but a higher amplitude of A a = 1.02 MPa causes irreversible bubble collapse after multiple rebounds. This indicates that the ADV becomes more complicated under high frequency conditions and the bubble growth regimes 2 and 3 are not clearly distinguished by the acoustic amplitude, which was also observed by Lacour et al. [27]. The complete vaporization occurs in the acoustic amplitude range of A a ⩾1.03 MPa, as depicted in Fig. 7a. The effective ADV threshold of A eth = 1.03 MPa and the direct ADV threshold of A dth = 3.95 MPa are much higher than those at f a = 2.25 MHz. Fig. 7b presents the acoustic pressure-frequency diagram for bubble growth regimes. In this work, the bubble growth regimes are divided into four regimes indicating (1) irreversible collapse with no rebound, (2) irreversible collapse after multiple rebounds, (3) complete vaporization after collapses and rebounds and (4) complete vaporization with no collapse. The acoustic pressure-frequency diagram for bubble growth regimes seems to be similar to that proposed by Lacour et al. [27], but there is a significant difference in regime 3 for complete vaporization after collapses and rebounds. Considering the acoustic amplitude range of A eth ⩽A a ⩽A dth at a high frequency of f a = 10MHz, where bubble growth is large in the first negative pulsing period, as seen at A a = 3.50 MPa of Fig. 6, modeling of bubble collapse is expected to be very important for the subsequent bubble behavior. The present analysis of the rapid collapse of a large bubble, which occurs at the supercritical state (Fig. 4d), results in the rebound and complete vaporization of the bubble, whereas the analysis of Lacour et al. [27], extending the phasechange formulation to the supercritical state, resulted in the irreversible collapse of the bubble with a very high temperature of a few MKs (regime IIb of Ref. [27]). The experimental data for the ADV threshold of Aliabouzer et al. [16] are within complete vaporization regime 3 and are more comparable to the present prediction of A eth , indicated by the red line in Fig. 7b, than the direct ADV threshold A dth . This indicates that the present numerical model for ADV is improved by properly treating the supercritical state at bubble collapse. Fig. 8 shows ADV at a bubble nucleation time changed to t bo = 0, which was used in the previous studies [23,27]. If a bubble nucleates at t bo = 0 when the acoustic pulsing starts, its pressure is estimated from Eq. (31) including the Laplace pressure contribution as p bo = 0.552 MPa. The corresponding bubble temperature of T bo = 363 K evaluated from Eq. (25) is higher than T ∞ = 310 K, and thus the nucleated bubble initially condenses until its temperature drops below T ∞ by negative acoustic pulsing.

Effects of bubble nucleation time and initial bubble radius
For a low frequency of f a = 2.25MHz, the bubble growth regimes are sharply distinguished as depicted in Fig. 8a. The bubble collapses with no rebound in the acoustic amplitude range of A a ⩽1.18 MPa, whereas the bubble directly vaporizes with no collapse for A a ⩾1.19 MPa. The effective and direct ADV thresholds are the same as A eth = A dth = 1.19 MPa, which is much higher than p eth = 0.49 MPa for ADV at the bubble nucleation time determined in Section 3.1.
The acoustic pressure-frequency diagram for bubble growth regimes  [16]., respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 9. Effect of initial bubble radius on the ADV threshold with R do = 0.47 μm. The circle symbols represent the experimental data [16].
is plotted in Fig. 8b. Compared with Fig. 7b, the bubble growth regimes are quite different at low frequencies although the difference becomes small at high frequencies. The prediction of A eth with t bo = 0 is deviated from the experimental data of Aliabouzer et al. [16] for ADV at a low frequency. This indicates that the present assumption of bubble nucleation when the bubble nucleus is in thermal equilibrium with the surrounding liquid is reasonable for analysis of ADV. Fig. 9 presents the effect of R bo on the ADV threshold. The ADV threshold pressure is observed to slightly increase as R bo decreases. However, in the range of 40nm⩽R bo ⩽120nm, the ADV threshold does not change much, and the tendency of the threshold increase with frequency does not change. This means that the choice of R bo = 80nm is not a very limiting condition for investigating the ADV threshold.

Effects of water properties and ambient temperature
Computations are performed to investigate the effects of water surface tension σ dw and water viscosity μ w on bubble growth in ADV. Such effects are implicitly included in the vaporization of droplets, which are typically phospholipid-coated or encapsulated in ADV applications. The encapsulation tends to reduce the surface tension and increase the viscosity effect due to the additional shell viscosity [35]. The analysis of Doinikov et al. [21] for ADV showed that the numerical predictions of bubble growth can be matched with the experimental data by varying the water surface tension and water viscosity to account for the encapsulation effect. Fig. 10a  The ADV results at four times increased water viscosity to μ w = 4 × 10 − 3 Pas are plotted in Fig. 10b. Compared to Fig. 7b for ADV at a lower water viscosity, the increased viscosity has a significant effect on A eth , but its effect is weak for regime 1 of irreversible bubble collapse without rebound. This is expected from the fact that the viscosity effect is remarkable in regimes with bubble oscillations of collapse and rebound.
The effective ADV threshold A eth for μ w = 4 × 10 − 3 Pas increases by about 140 %. Fig. 11a shows the effect of water density ρ w on the effective ADV threshold A eth . When f a ⩽7 MHz, A eth is not sensitive to ρ w . However, in a higher frequency range, A eth increases with ρ w . This means that the energy required for the bubble to push the heavier surrounding liquid increases under higher frequency conditions with more bubble oscillations. The influence of ambient temperature T ∞ on A eth is depicted in Fig. 11b. As T ∞ increases from 303 K to 310 K and 318 K, the effective the present numerical predictions of the ADV threshold and the experiment data [16]., respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 11. Influences of water density ρ w and ambient temperature T ∞ on the ADV threshold with R do = 0.47 μm. The circle symbols represent the experimental data [16].
ADV threshold decreases by about 10 %, respectively. This indicates that A eth decreases as the evaporation heat transfer increases with T ∞ . Fig. 12a shows bubble growth for ADV at the initial droplet radius increased to R do = 6 μm. As R do increases, the period of bubble growth is not complete during the acoustic pulsing periods of 8/f a . After the pulsing period, the bubble is observed to grow slowly because of the evaporative heat transfer from the surrounding liquid.

Effect of initial droplet radius
The effects of initial droplet radius R do and water viscosity μ w on the effective ADV threshold A eth are plotted in Fig. 12b. As the acoustic frequency f a increases, A eth is observed to increase regardless of the magnitude of R do and μ w . A eth increases with increasing μ w in the case of a small droplet, but it has little effect on μ w for a large droplet of R do = 6 μm. When R do increases from 0.47 μm to 6 μm, A eth decreases, as observed in the experiment of Aliabouzer et al. [16]. The present predictions of A eth are more comparable to the experimental data [16] at R do = 6 μm than R do = 0.47 μm. The comparison is better when μ w is increased by a factor of four.

Conclusions
A numerical model for the ADV threshold of microdroplets was developed by improving the Rayleigh-Plesset equation to properly treat the supercritical state occurring at rapid bubble collapse and by employing the van der Waals equation to consider the supercritical state.
The assumption of bubble nucleation when the bubble nucleus is in thermal equilibrium with the surrounding liquid is observed to be reasonable for analysis of ADV. The numerical results showed four distinguished bubble growth regimes of irreversible collapse or complete vaporization with no rebound or multiple rebounds. The effective ADV threshold is observed to increase as the acoustic frequency increases and the initial droplet radius decreases, as experimentally observed in the previous studies. The ADV threshold increases with increasing the water viscosity in the case of a small droplet, but it has little effect on the water viscosity for a large droplet. The present numerical predictions of ADV threshold comparable to the experimental data.

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.