Energy transfer through boiling-cavitation interaction in high-intensity focused ultrasound

2.1. Nonlinear sound field and bubble swarm homogenization scattering correction

The fluctuation of the bubble swarm integral considerably modulates the sound propagation process under the operation of a strong sound field. The updated Westervelt equation with the bubble source term is presented in order to explain the harmonic dispersion and the dynamic shielding impact of the bubble swarm on the sound field [13]:

The sound pressure ($p$), sound velocity ($c$) in the liquid phase, sound diffusivity ($\delta$), nonlinear coefficient ($\beta$), medium density ($\rho$), and local bubble volume fraction (${\alpha }_v(r,t)$) are all represented in Eq. (1). When combined with Wood’s homogenization theory, the abrupt decrease in $c_m$ results in a significant acoustic impedance mismatch that builds a backscattering barrier as the bubble concentration rises [14].

2.2. Vaporization phase change heat transfer model based on pure enthalpy method

This study builds a thermodynamic model using the pure enthalpy method instead of the conventional sensible heat equation in order to precisely represent the energy conservation and temperature clamping effect under superboiling point conditions [15]. With the medium's transient volumetric enthalpy $H(r,t)$ serving as the fundamental solution variable, the energy equation is rewritten as follows:

At the mathematical level, the energy conservation of the phase transition process is rigorously guaranteed by Eq. (2). The temperature $T$ and the volumetric enthalpy $H$ in this analytical mapping logic satisfy the following state relationship: the medium is a pure liquid phase when $H<\rho C_p{T}_{sat}$; it enters the boiling latent heat plateau when the local total enthalpy $H$ approaches the liquid phase limit ($\rho C_p{T}_{sat}$). At this point, the remaining acoustic energy will be entirely transformed into molecular phase transition potential energy, and the local temperature is strictly limited to $T\equiv T_{sat}$ and maintained there until the latent heat of vaporization is fully depleted. The enthalpy then enters the superheated vapor phase when it surpasses the volumetric latent heat of vaporization barrier $\mathrm{\Delta }{H}_{vap}=\rho L$. Here, $T_{sat}$ stands for the medium’s equilibrium saturation temperature, which in this paper's atmospheric pressure model is strictly assumed to be 100 °C. At the mathematical level, this treatment accurately describes the isothermal plateau characteristic that is specific to the boiling phase transition. It should be mentioned that there is a tight mathematical coupling between thermodynamics and acoustic dynamics because the local volumetric heat generation rate of the coke area and the effective sound intensity strictly satisfy the energy conservation mapping: $Q_{acoustic}=2{\alpha }_{eff}{I}_{eff}$, where ${\alpha }_{eff}$ is the nonlinear equivalent sound absorption coefficient.

The three-dimensional heat transport is reduced to a zero-dimensional lumped parameter model in order to get around the computational stiffness of very long-term simulations. The tissue thermal relaxation time constant can be mathematically determined as $\tau ={\omega }_0^2/\left(4{\alpha }_T\right)$ based on the characteristic waist radius of the HIFU Gaussian focusing spot, ${\omega }_0=$ 1.0 mm. The thermal diffusivity ${\alpha }_T=k/\left(\rho C_p\right)$ ≈ 1.26×10⁻⁷ m²/s and $\tau$ ≈ 2.0 s are determined by substituting the tissue thermal conductivity and specific heat capacity parameters from Table 1.

The macroscopic energy conservation of the focal region can be reduced to an ordinary differential form, $dH/dt=Q_{acoustic}-\rho C_p(T-T_{init})/\tau$, if this relaxation time approximation is introduced. In addition to overcoming computational stiffness, this dimensionality reduction perfectly closes the mathematical self-consistency of heat transport dissipation at the physical level. Additionally, this model ignores the endothermic reaction of tissue protein denaturation in order to maintain physical self-consistency because it is almost two orders of magnitude lower than the latent heat of vaporization (2.26×10⁶ J/kg).

2.3. Bubble dynamics equation under viscoelastic constraints in soft tissue

Biological tissues’ viscoelastic matrix greatly prevents bubbles from expanding excessively during the negative pressure phase and gives the bubble wall a powerful elastic rebound driving force during the positive pressure contraction phase. The elastic deviatoric stress integral is strictly correlated with the shear modulus $G=E/\left[2\right(1+\nu \left)\right]$ of the medium because bubble expansion and contraction cause shear deformation on the surrounding tissue. The microcompressible Keller-Miksis equation with the Kelvin-Voigt solid constitutive term was developed in order to examine the dynamic properties of the collapse stage [16]:

Eq. (3) shows that the bubble’s total internal pressure is $p_B$=$(p_0+(2\sigma /R_0)-p_v)(R_0/R{)}^{3\gamma }+p_v$, where $p_v$ is the saturated vapor pressure, which exerts positive work outward during the negative pressure expansion phase. $\gamma$ is the polyhedral index of the gas inside the bubble; $\mu$ and $\sigma$ are the liquid phase’s dynamic viscosity and surface tension coefficient, respectively; $R$ is the bubble’s instantaneous radius; $\dot{R}$ and $\ddot{R}$ are the bubble wall velocity and acceleration, respectively; $E$ is the tissue’s Young’s modulus; $\upsilon$ is Poisson’s ratio (soft tissue is regarded as a completely incompressible solid matrix, strictly taken as 0.5); and $R_0$ is the bubble’s initial radius. It should be noted that the viscoelastic deviatoric stress, viscous dissipation, and surface tension terms in Eq. (3) adopt a low Mach number approximation, i.e., ignore their extremely weak acoustic radiation delay derivative terms, and concentrate the core nonlinear effect of acoustic damping on the bubble's internal pressure $p_B$ and the far-field driving pressure $p_{\infty }\left(t\right)$, in order to balance the numerical stability and physical fidelity of the strong nonlinear collapse period. This method carefully upholds the kinetic equation’s mathematical well-posedness while guaranteeing the accuracy of thermodynamic evolution. During the bubble growth phase, the deviatoric shear constraint term on the right side of Eq. (3) functions as an inward constraint resistance. However, the massive shear elastic potential energy that is instantly released by the highly stretched viscoelastic matrix during the collapse phase ($R<R_0$) reverses the dominance of viscous damping and is converted into strong kinetic energy that propels the bubble wall's centripetal acceleration. This is the primary dynamic source that causes the bubble wall to collapse supersonically ($Ma>1$).

2.4. Composite thermo-mechanical bimodal damage criterion

The bimodal damage field was treated to dimensionless normalization based on the membrane rupture mechanical critical value (${\widehat{\mathrm{\Omega }}}_{thermal}$ and ${\widehat{\mathrm{\Omega }}}_{mech}$, respectively) and the cell coagulation heat threshold (CEM_43℃) in order to reduce the disparity in the dimensions of thermo-mechanical integrals [12]. The disaster of computational dimensionality would result from directly doing discrete integration on the tiny collapse Mach number of a collection of billions of bubbles on a gigantic spatiotemporal scale. Thus, this work presents an energy-based continuous integral surrogate model based on cavitation dosimetry theory.

A highly sensitive continuous phase transition probabilistic surrogate model was developed in order to mathematically describe finite-rate condensation kinetics and prevent numerical divergence across boiling points:

The delayed condensation characteristics caused by the extreme thinning and collapse of the thermal boundary layer in a local subcooled liquid environment are characterized by the condensation threshold temperature $T_{th}=$ 98 ℃ (slightly lower than the standard pressure boiling point) in Eq. (4); the phase change sensitivity coefficient $k_s=$ 5.0$K^{-1}$ strictly maps the intrinsic time scale of finite-rate condensation of bubbles dissolving rapidly in a very narrow temperature range (≈±1 ℃). This proxy model firmly avoids the non-physical numerical divergence brought on by the temperature beyond the boiling point at the mathematical level, in addition to physically fitting the actual vapor bubble dissolving process. The continuous cavitation survival probability density $P_{cav}\left(T\right)$ modulates the mechanical damage equivalent, which is adjusted to be proportional to the square of the local nonlinear acoustic energy dissipation:

where, $\mathrm{\Pi }\left(t^{\text{'}}\right)$ is a rectangular wave gating function strictly modulated by a low-frequency pulse sequence (PRF); ${\widehat{\mathrm{\Omega }}}_{mech}$ is the phenomenological scaling coefficient, whose dimension is strictly limited to m⁴/(W².s), which is used to forcefully cancel the dimension of the integral term and convert it into a dimensionless relative damage equivalent. The fundamental mechanism of discrete step-like stagnation of mechanical damage is constructed at the mathematical level when $\mathrm{\Pi }\left(t^{\text{'}}\right)$ = 1 during the pulse-on period (ON-time) and $\mathrm{\Pi }\left(t^{\text{'}}\right)$ = 0 during the pulse-off period (OFF-time). The composite bimodal damage field is defined as the spatial extremum envelope following normalization of the two in order to precisely describe the lesion’s ultimate three-dimensional geometric topology:

The biological logic that tissue necrosis is declared when it crosses the deadly threshold under any single thermal or mechanical mode is physically represented by this definition, which carefully avoids the direct algebraic addition of many physical dimensions. The continuous wave (CW) control group in this study is defined as the continuous irradiation mode without duty cycle intervention. The maximum damage extreme value in the CW mode is used to normalize the mechanical damage equivalent.

3.1. The clamping effect of boiling phase change latent heat on focal temperature

The focal center temperature’s transient evolution under a high-intensity continuous sonic radiation load is shown in Fig. 2. After about a second of irradiation, the classical heat conduction model (red dashed line), ignoring phase transition, surpasses 100 °C. This is followed by a linear surge without physical boundaries, reaching a maximum temperature above 160 °C. Introducing the pure enthalpy phase transition coupling model of Eq. (2) (blue solid line), a significant amount of ultrasonic energy is required to overcome the hydrogen bond binding of water molecules and cross the volumetric latent heat of vaporization barrier when the temperature approaches the saturation temperature. The focal temperature is firmly restricted to about 100 °C, and the temperature curve enters a clear boiling plateau period. From a thermodynamic point of view, this mechanism verifies that the massive volumetric latent heat absorbed by the liquid-gas phase transition functions as a natural temperature buffer, thereby limiting the thermal runaway border. Note that the maximal acoustic energy deposition rate ($Q_{max}$ = 4.1×10⁸ W/m³) employed in this study is not based on real data. The equivalent sound absorption coefficient (${\alpha }_{eff}$) of the coke zone is significantly increased to roughly 54 Np/m under extreme sound pressure of $P\left(A\right)$ = 3.5 MPa due to the nonlinear distortion of the sound wave, which creates high-order harmonic cascades [19]. This extreme value can be rationally obtained by substituting it into the coke zone energy conservation equation $Q={\alpha }_{eff}{P}_A^2/\left(\rho c\right)$, which appropriately characterizes the explosive thermal deposition in the HPSP mode.

Fig. 2The clamping effect of latent heat of boiling phase change on the temperature of the HIFU coke zone

3.2. Accelerating effect of tissue viscoelasticity on bubble inertial collapse

In an isolated porcine liver tissue environment, Fig. 3 illustrates the dynamic properties of a single cavitation nucleus traveling through a full sonic cycle. The bubble grows to about 2.5$R_0$ during the negative pressure phase; during the positive pressure phase, dramatic inertial collapse takes place, with the bubble wall velocity surpassing the liquid phase sound velocity ($Ma$ > 1.2).

The nonlinear expansion and contraction properties of bubbles under HIFU negative pressure were validated by the single-bubble dynamics simulation (Fig. 3). This work examined the collapse Mach number under isolated porcine liver tissue constraint and degassed water free field settings in order to further assess the impact of the viscoelastic matrix of living tissue on the collapse process (Fig. 4). The numerical rigidity issue during the collapse peak was addressed using an implicit variable step-size technique based on the backward differential formula. The findings demonstrated that the peak collapse Mach number in the tissue environment was substantially larger than that in the degassed water free field, and that the tissue shear elastic potential energy released during the severe contraction phase reversed the dominance of viscous damping. This unusual “elastic rebound mechanism” challenges the conventional wisdom that “high viscosity necessarily inhibits cavitation”, implying that solid tumors with greater stiffness are more vulnerable to cavitation-related mechanical damage. The sensitivity of the collapse dynamics to the tissue shear modulus ($G$) needs to be further discussed in order to increase the universality of this statement. Although isolated porcine liver was represented in this study using a fixed Young’s modulus of $E$ = 15 kPa (equivalent to $G$ = 5 kPa for incompressible tissue with $\upsilon$ = 0.5), actual clinical solid tumors show significant biomechanical heterogeneity, with stiffness typically ranging from 5 kPa to 50 kPa. The shear elastic deviatoric stress is exactly proportional to the modulus according to the theoretical framework of Eq. (3). Therefore, a parameter sensitivity analysis shows that the viscoelastic matrix will store significantly more shear elastic potential energy during the negative-pressure expansion phase as the target tissue stiffness climbs near 50 kPa (e.g., highly fibrotic pancreatic or breast cancers). This results in an even higher peak collapse Mach number during the collapse phase due to a significantly increased elastic rebound driving force. Conversely, for softer tissues approaching 5 kPa, this acceleration effect lessens, and the collapse dynamics smoothly converge toward the fluid-like behavior of the water free field. In addition to confirming the robustness of the suggested mechanism, this strong parameter sensitivity offers an essential biomechanical theoretical foundation for customizing HIFU dosimetry depending on tumor stiffness.

Fig. 3Nonlinear inertial collapse of air bubbles in isolated porcine liver tissue under HIFU pulse irradiation

Fig. 4Comparison of the accelerating effect of tissue viscoelasticity on the inertial collapse phase of bubbles

3.3. Spatial dislocation of composite damage induced by cavitation bubble group acoustic shielding

Using a one-dimensional spatial coupling model, Fig. 5 depicts the axial separation characteristics of the dual-modal damage field: the thermal damage considerably decreases behind the focal point ($z$ > 10 mm) due to acoustic shielding, resulting in trailing edge thermal underkill, while the mechanical damage equivalent forms a peak region at the leading edge of the focal point (8~10 mm). It should be made clear that Fig. 5 performs morphological normalization of the thermal and mechanical damage equivalents based on their respective spatial peak values, which is only used to characterize the relative shift of the spatial envelope, in order to clearly reveal the topological dislocation law of the two damage modes along the spatial axis.

However, a global absolute normalization based on the same uniform death physical threshold is used in the following temporal development (Fig. 7) and two-dimensional distortion reconstruction (Fig. 6) to guarantee rigorous dimensional consistency of the dual-modal superposition. ${\alpha }_b=$ 2.0 Np/mm was chosen as the bubble's equivalent acoustic attenuation coefficient. The significant multiple scattering cross-section properties of dense cavitation clouds under the HPSP model were used to estimate this value, which is extremely comparable with the strong echo spot attenuation magnitude observed in current live tests [20].

Furthermore, the spatial damage proxy model is strictly unified as a nonlinear dependence on the square of the effective acoustic intensity ($I_{eff}^2$), which guarantees the absolute isomorphism of cross-dimensional simulation from the mathematical bottom layer, and the cavitation probability function in the one-dimensional and two-dimensional models shares the same set of intrinsic parameters (threshold $I_{th}$ = 0.35, steepness 20).

Fig. 5Axial dual-modal damage spatial misalignment induced by bubble cloud acoustic shielding

The damage field’s axial misalignment trend is evident in the one-dimensional axial distribution (Fig. 5), but this is not enough to describe the lesion’s two-dimensional geometry. This study developed an axisymmetric sound-thermal coupling model (Fig. 6), taking into account the asymmetric shielding effect of the bubble cloud in the direction of sound propagation. The contour lines of the composite damage field in Fig. 6(b) clearly show the swollen “tadpole head” (overkill zone at the front edge) and the narrow “tadpole tail” (underkill zone at the rear field); the area circled by the red dot is the mechanical overkill zone caused by cavitation collapse. Fig. 6(a) illustrates how the strong nonlinear sound attenuation of the bubble cloud causes the effective sound energy field to shift forward to the shallow part. The usual lesion morphology seen under clinical B-ultrasound imaging is in good accord with the numerical results [21].

Fig. 6Two-dimensional axisymmetric acoustic field shielding effect and spatial distortion reconstruction of “tadpole-shaped” composite damage

a) Acoustic shielding and energy shift

b) Tadpole composite lesion distortion

3.4. Thermo-mechanical dual-mode damage optimization based on dynamic duty cycle feedback control

This paper further suggests an optimization strategy for dynamic duty cycle control based on the previously mentioned discovery of the spatial misalignment problem of composite damage, namely the coexistence of mechanical overkill at the leading edge of the target area and thermal underkill at the trailing edge. The choice of a 20 % duty cycle is not based on empirical assumptions, but rather on a rigorous analytical solution of the system's thermal relaxation characteristics: $DC$ = 0.2 is the only optimal solution that precisely anchors the asymptotic temperature $T_{ss}$ within the tissue lethal range (60-80 ℃) while reserving enough pulse intervals to ensure smooth condensation of cavitation bubbles under the constraints of $\tau$ = 2.0 s and $Q_{max}$ = 4.1×10⁸ W/m³. In Fig. 7, a comparative scenario is simulated in which the system modulates the actual pulse repetition frequency (PRF = 10 Hz) and adds a heat conduction dissipation term based on the thermal relaxation time constant ($\tau$ = 2.0 s), automatically reducing the duty cycle from 100 % to 20 % once the focal temperature reaches the 100 °C threshold. The focal temperature displays a true microscopic sawtooth fluctuation following feedback control (Fig. 7(a)). The dynamic equilibrium between heat conduction and pulse heating causes the temperature to stabilize at about 78.4 °C (purple dashed line), which is still above the tissue coagulation and necrosis threshold (60 °C), demonstrating the maintenance of the effective thermal ablation effect. It should be mentioned that the focal temperature’s ultimate stabilization at 78.4 °C is not a numerical coincidence. The analytical steady-state asymptote of the system can be strictly written as $T_{ss}=T_{init}+DC\cdot {Qp}_{max}$ once the duty cycle ($DC$ = 0.2) is introduced, according to the lumped parameter energy conservation equation. The theoretical calculation value precisely points to 78.43 °C when the intrinsic parameters derived in this study ($\tau$ = 2.0 s, $Q_{max}$ = 4.1×10⁸ W/m³) are substituted, achieving perfect cross-scale self-consistency with the numerical answer. The drop in duty cycle successfully breaks the boiling-cavitation positive feedback and actively departs the boiling state, as seen by the mechanical damage accumulation curve in Fig. 7(b), which is globally normalized based on the highest damage extremum in CW mode. The mechanical overkill accumulation at the leading edge of the target region exhibits a smooth asymptotic decrease and tends to halt when the coke temperature falls outside of the boiling danger zone, resulting in accurate physical melting of mechanical damage. The normalized steady-state extremum of mechanical damage in CW mode is 1.0, as indicated by the quantitative integration results in Fig. 7(b); however, the steady-state final value in PWM mode is rigidly locked at 0.0057 after initiating 20 % dynamic duty cycle regulation. This indicates that our discrete integration technique based on transient sound pressure gating physically and precisely lowers (1 – 0.0057 = 99.43 %) the ineffective mechanical overkill when compared to the CW benchmark.

Fig. 7Comparison of optimization effects of dynamic duty cycle feedback control on thermo-mechanical dual-mode damage

a) Dynamic thermal equilibrium

b) Staircase arrest of mechanical damage

The mechanical damage curve following duty cycle intervention in Fig. 7(b) exhibits a microscopic step accumulation characteristic, which is important for assessing the intervention effect of dynamic duty cycle. Heat transfer lag causes a quick condensation period in the tissue temperature when the system initiates a 20 % duty cycle modification. The remaining vapor bubbles in this incredibly brief buffer window only undergo violent inertial collapse (shown as a vertical rise in the curve) within 0.02 seconds of pulse activation ($\mathrm{\Pi }\left(t^{\text{'}}\right)$ = 1), and the violent collapse of microbubbles momentarily stops (shown as an absolute horizontal line) within 0.08 seconds of pulse deactivation ($\mathrm{\Pi }\left(t^{\text{'}}\right)$ = 0) because the transient sound pressure returns to zero. Calculated by combining strict transient sound pressure gating with a highly sensitive phase change condensation model, this discrete step envelope profoundly reveals the true fluid dynamics signature of low-frequency pulsed ultrasound’s capacity to achieve “precise start-stop control” of tissue damage rather than being a numerical artifact. The need for microsecond-level wave pauses following the emergence of “strong echo spots” in clinical ultrasonography is quantitatively demonstrated by this study.