Numerical modelling and theoretical analysis of the acoustic attenuation in bubbly liquids

Abstract

The propagations of acoustic waves in bubbly liquids have been extensively investigated through experimental and theoretical methods, a computational fluid dynamics (CFD) method was introduced to investigate the acoustic attenuation through bubbles in this study. Numerical ‘measurements’ of attenuation coefficient (α) and phase velocity (V) were conducted using a homogeneous cavitation model, which were modeled from experimental schemes and compared with the theoretical results. At extremely small bubble volume fraction (β around 10⁻¹¹), new formulas of the theoretical α (α_theo) were proposed respectively for linear and transient oscillations, and a new formula of the numerical α (α_num) was proposed for transient oscillations. Results showed that α_num matched precisely with α_theo for linear, nonlinear and transient oscillations. At medium β (around 10⁻⁴), the relative difference of α_num between the VOF and present methods was less than 1.6%, while it reached 15.4% after replacing the bounded Keller-Miksis equation (KME) in the present method with the KME. However, the traditional theoretical α and V matched precisely with the predictions by the present method with the KME. Thus new theoretical α and V were proposed based on the bounded KME, and the relative differences between α_theo and α_num were less than 1%. It can be concluded that the bounded KME should be used in both numerical and theoretical predictions.

KEYWORDS:

• Ultrasonic cavitation
• numerical modelling
• attenuation
• phase velocity
• transient analysis

Nomenclature

α	=	attenuation coefficient
αtheo	=	theoretical α
αlin	=	αtheo of linear oscillations
αnonlin	=	αtheo of nonlinear oscillations
αtran	=	αtheo of transient oscillations
αm	=	αtheo at medium β using the KME
αmb	=	αtheo at medium β using the bounded KME
αnum	=	numerical α
αnumt	=	numerical α of transient oscillations
αexp	=	experimental α
αrad	=	α due to radiation
β	=	bubble volume fraction
σ	=	surface tension coefficient
σs	=	scattering cross section
δ	=	damping coefficient
d	=	acoustic path length
n	=	bubble number density
pA	=	amplitude of the acoustic excitation

1. Introduction

Ultrasonic cavitation has been widely applied in biomedical ultrasonic diagnosis and treatment (Helfield et al., 2021; Stride et al., 2020). The former application (contrast-enhanced ultrasound imaging) uses ultrasound contrast agents (UCAs) to increase the echo difference between the examination part and the surrounding tissues, so as to obtain clearer ultrasonic images (Yusefi & Helfield, 2022). The latter application often introduces UCAs as the cavitation nuclei and needs a high-intensity ultrasound to cause violent collapse, so as to achieve the purpose of drug/gene delivery (Alishiri et al., 2021; Yang et al., 2020), cancer treatment (Clark et al., 2021), sonothrombolysis (Nederhoed et al., 2021), and so on. UCAs are coated microbubbles with diameters ranging from 0.7–10 µm. They are made of high-molecular-weight inert gases and coated by lipid, protein, or polymer shells (Frinking et al., 2020) to improve their stability.

The Rayleigh-Plesset (RP) equation (Plesset & Prosperetti, 1977) provided the theoretical prediction of the dynamics of a spherical uncoated bubble in an unbounded liquid. Several extensions of the RP equation have been proposed to take into account the liquid compressibility (Brenner et al., 2002), such as the famous KME (Keller & Miksis, 1980). Bubble oscillations exhibit nonlinear behaviours under middle or high-intensity ultrasound, and the nonlinear characteristics are enhanced by the bubble shells of UCAs (Gümmer et al., 2021; Omoteso et al., 2021). It has been more than two decades of progress in the modelling of UCAs (Helfield, 2019; Versluis et al., 2020). By adding pressure terms to the RP type equations, many models (Allen et al., 2002; Cowley & McGinty, 2019; de Jong et al., 1994; Marmottant et al., 2005; Paul et al., 2010) have been proposed to consider the influence of bubble shells. The influence is mainly described by two terms: the shell elasticity and the dilatational viscosity (κ_s). The Marmottant model (Marmottant et al., 2005) is widely used and can well predict the nonlinear oscillation behaviours, such as the ‘compression-only’ and ‘thresholding’ (Versluis et al., 2020). In this model, the shell elasticity was embedded in the surface tension (σ), which was discontinuous when the buckling or rupture occurred. Paul et al. (2010) adopted an exponential elasticity model to make a smooth transition from σ_break−up to σ_L. κ_s increases with the ambient bubble radius (R₀), and it is usually given as a constant in the range of 10⁻¹⁰ kg/s to 10⁻⁸ kg/s (Helfield & Goertz, 2013; Tu et al., 2009; Tu et al., 2011).

The measurements of α can be used to obtain the shell parameters of UCAs. UCAs are around resonances under the most common frequencies of medical ultrasound (1-7 MHz) (Gorce et al., 2000), thus their dynamics are greatly dependent on R₀. For small-amplitude oscillations, RP type equations can be simplified to linear equations of motion, based on which α_theo can be obtained. Gorce et al. (2000) divided UCAs into groups and theoretically analysed the influence of R₀ on echogenicity and attenuation based on the simplification of linear oscillation. Segers et al. (2016) measured the scattering and attenuation of sorted UCAs at higher driving pressure, and calculated α_theo by solving the RP type equation to capture the nonlinear behaviours. At medium β, the dependences of α and V on β need to be considered. Silberman (1957) conducted experimental measurements of α and V in a standing wave tube at β above 0.03%. Wijngaarden (1972) proposed a continuum approach to model the wave propagation in bubbly liquids. Commander and Prosperetti (1989) developed a linear model to predict α and V at small driving pressure and medium β, and good agreements were obtained with the experiments by Silberman (1957) et al. Leroy et al. (2009) experimentally measured the transmission coefficients of ultrasound through a single layer of bubbles and found that the resonant frequency increased with the bubble concentration, which was contrary to the conclusion of traditional theoretical analyzes (Ida et al., 2007; Yasui et al., 2008; Yasui et al., 2009).

Numerical simulations have been widely used to investigate cavitation (Fuster, 2019; Modarreszadeh et al., 2021; Si et al., 2023; Subroto et al., 2021; Suo et al., 2021; Yasui, 2021). The interface capturing methods, such as the VOF method (Chen et al., 2021; Lee et al., 2021; Li et al., 2022; Yamamoto & Komarov, 2021; Ye et al., 2021b), can accurately predict bubble dynamics. However, numerous grids are needed to capture each bubble, and it is difficult to model the bubble shell of UCAs. Homogeneous cavitation models are almost the only way for the predictions of cavitation with numerous bubbles. Ye et al. (Ye et al., 2020a; Ye et al., 2020b) proposed a homogeneous cavitation model which can accurately predict the acoustic cavitation of the equal-sized spherical bubbles in a stationary liquid. Due to the dependency of bubble dynamics on R₀, Ye et al. (2021a) proposed a homogeneous grouping cavitation model which divided bubbles into groups, and the dynamics of each group were separately calculated.

Since the bubble oscillations can be well predicted by the homogeneous cavitation model (Ye et al., 2021a), the acoustic attenuation through bubbles should also be well predicted. Numerical simulations can not only investigate the complex cavitating flows, but also help to improve the basic theory. This study conducted theoretical and numerical investigations on the attenuation and phase velocity through uncoated and coated bubbles. At extremely small β, the inter-bubble interaction can be neglected. α_theo caused by oscillating bubbles in the linear (α_lin2) and transient (α_tran) states were obtained, and several simplified forms of α_theo for linear and nonlinear oscillations were introduced in Section 2. The numerical ‘measurements’ of α were implemented by comparing the recorded pressure with and without bubbles in the acoustic beam, which were modeled from experimental schemes, and comparisons were made with α_theo in Section 3. At medium β, the present method was validated by comparing with the VOF method, and then comparisons were made with the experimental (Silberman, 1957) and theoretical (Commander & Prosperetti, 1989) results of millimetre bubbles, based on which the theoretical formulas were improved (Section 3.4). At last, α and V through microbubbles were numerically measured at medium β, and the influences of the nonlinearity and transience were numerically analysed.

2. Methods

2.1. Bubble dynamics equation

Neglecting the liquid compressibility, the phase change, and the diffusion of non-condensable gas into liquid, the dynamics of a spherical bubble under far-field acoustic excitation can be described by the RP type equation (Plesset & Prosperetti, 1977): (1) $\begin{array}{rl}R\ddot{R}+\frac{3}{2}{\dot{R}}^2& =\frac{1}{{\rho }_{\text{L}}}[p_{\text{G}}-\frac{2\sigma }{R}-\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}})\dot{R}}{R}\\ & -p_0-p_{\text{a}}(t)]\end{array}$(1) (2) $\begin{array}{rl}{p}_{\text{G}}& =p_{\text{G0}}{\left(\frac{R_0}{R}\right)}^{3\kappa }\end{array}$(2) where R is the bubble radius at time t, dots denote time derivatives, ρ_L is the liquid density, p_G is the gas pressure inside the bubble, the subscript 0 denotes the initial equilibrium value without excitation, σ is the surface tension coefficient, µ_L is the liquid viscosity, µ_th is the effective thermal viscosity, p_a(t) is the acoustic excitation, and κ is the polytropic exponent. The famous KME (Keller & Miksis, 1980) considering the liquid compressibility was given as: (3) $\begin{array}{rl}& \left(1-\frac{\dot{R}}{c_{\text{L}}}\right)R\ddot{R}+\frac{3}{2}\left(1-\frac{\dot{R}}{3c_{\text{L}}}\right){\dot{R}}^2\\ & =\frac{1}{{\rho }_{\text{L}}}\left(1+\frac{\dot{R}}{c_{\text{L}}}+\frac{R}{c_{\text{L}}}\frac{d}{dt}\right)\\ & \times \left(p_{\text{G}}-\frac{2\sigma }{R}-\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}})\dot{R}}{R}-p_0-p_{\text{a}}(t)\right)\end{array}$(3) where $c_{\text{L}}$ is the speed of sound in liquid. Substituting Eq. (1) into Eq. (3), another form of the KME can be obtained (Zhang, 2013): (4) $\begin{array}{rl}& \left(1-\frac{2\dot{R}}{c_{\text{L}}}+M\right)R\ddot{R}+\left(\frac{3}{2}-\frac{2\dot{R}}{c_{\text{L}}}-M\right){\dot{R}}^2\\ & =\frac{1}{{\rho }_{\text{L}}}[p_{\text{G}}-\frac{2\sigma }{R}-\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}})\dot{R}}{R}\\ & -\left(3\kappa p_{\text{G}}-\frac{2\sigma }{R}\right)\frac{\dot{R}}{c_{\text{L}}}-\frac{{\dot{p}}_{\text{a}}(t)R}{c_{\text{L}}}-p_0-p_{\text{a}}(t)]\end{array}$(4) (5) $\begin{array}{rl}& M=\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}})}{{\rho }_{\text{L}}{c}_{\text{L}}R}\end{array}$(5) M ≪ 1 for most oscillations, and $|\dot{R}|/c_{\text{L}}$≪ 1 for small-amplitude oscillations, thus a popular simplified form can be obtained by further neglecting the time differentials of the surface tension and $p_{\text{a}}(t)$ in Eq. (3) (Brenner et al., 2002): (6) $\begin{array}{rl}R\ddot{R}+\frac{3}{2}{\dot{R}}^2& =\frac{1}{{\rho }_{\text{L}}}[p_{\text{G}}\left(1-\frac{3\kappa \dot{R}}{c_{\text{L}}}\right)-\frac{2\sigma }{R}\\ & -\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}})\dot{R}}{R}-p_0-p_{\text{a}}(t)]\end{array}$(6) Marmottant et al. (2005) proposed an equation for coated bubbles based on Eq. (6) (Segers et al., 2016; Segers et al., 2018): (7) $\begin{array}{rl}R\ddot{R}+\frac{3}{2}{\dot{R}}^2& =\frac{1}{{\rho }_{\text{L}}}[p_{\text{G}}\left(1-\frac{3\kappa \dot{R}}{c_{\text{L}}}\right)-\frac{2\sigma (R)}{R}\\ & -\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}}+{\kappa }_{\text{s}}/R)\dot{R}}{R}-p_0-p_{\text{a}}(t)]\end{array}$(7) (8) $\begin{array}{rl}{p}_{\text{G0}}& =p_0+\frac{2\sigma (R_0)}{R_0}\end{array}$(8) (9) $\begin{array}{rl}\sigma (R)& =\{\begin{array}{l}\chi max\left({\displaystyle \frac{R^2}{R_{\text{b}}^2}}-1,0\right)R\le R_{\text{ruptured}}\sigma R>R_{\text{ruptured}}\end{array}\end{array}$(9) where κ_s is the dilatational viscosity of the shell, χ is the elastic modulus, R_b is the buckling radius, and R_ruptured is the rupture radius. Eqs. (6) and (7) are named as the simplified KME in this study. The complete KME for coated bubbles can be obtained by substituting the Marmottant model (Marmottant et al., 2005) into Eq. (3): (10) $\begin{array}{rl}& \left(1-\frac{2\dot{R}}{c_{\text{L}}}+M\right)R\ddot{R}+\left(\frac{3}{2}-\frac{2\dot{R}}{c_{\text{L}}}-M-\frac{4{\kappa }_{\text{s}}}{{\rho }_{\text{L}}{c}_{\text{L}}{R}^2}\right){\dot{R}}^2\\ & =\frac{1}{{\rho }_{\text{L}}}[p_{\text{G}}-\frac{2\sigma (R)}{R}-\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}}+{\kappa }_{\text{s}}/R)\dot{R}}{R}\\ & -\left(3\kappa p_{\text{G}}+\frac{2\sigma (R)}{R}+\frac{4\chi }{R}\right)\frac{\dot{R}}{c_{\text{L}}}-\frac{{\dot{p}}_{\text{a}}(t)R}{c_{\text{L}}}\\ & -p_0-p_{\text{a}}(t)\phantom{\frac{4{\kappa }_{\text{s}}}{{\rho }_{\text{L}}{c}_{\text{L}}{R}^2}}]\end{array}$(10) (11) $\begin{array}{rl}& M=\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}}+{\kappa }_{\text{s}}/R)}{{\rho }_{\text{L}}{c}_{\text{L}}R}\end{array}$(11)

2.2. Linear equation and damping coefficients

For the small-amplitude oscillation about the equilibrium radius R = R₀(1 + x) (x ≪ 1), the KME (Eq. (4) or (10)) can be linearised to the classic equation by neglecting high-order terms (Leighton, 1994; Versluis et al., 2020): (12) $\ddot{x}+{\delta }_{\text{tot}}{\omega }_0\dot{x}+{\omega }_0^{2}x=-\frac{p_{\text{A}}}{{\rho }_{\text{L}}{R}_0^2}\sin (\omega t)$(12) where δ_tot is the total damping coefficient, p_A is the amplitude of the acoustic excitation, and ω₀ is the angular eigenfrequency. ω₀ of uncoated and coated bubbles can respectively be obtained as: (13) $\begin{array}{rl}{\omega }_0^2& =\frac{3\kappa p_0+2(3\kappa -1)\sigma /R_0}{{\rho }_{\text{L}}{R}_0^2(1+M(R_0))}+\frac{{\delta }_{\text{tot}}{\omega }_0{\omega }^2{R}_0}{c_{\text{L}}}\end{array}$(13) (14) $\begin{array}{rl}{\omega }_0^2& =\frac{3\kappa p_0+2(3\kappa +1)\sigma (R_0)/R_0+4\chi /R_0}{{\rho }_{\text{L}}{R}_0^2(1+M(R_0))}\\ & +\frac{{\delta }_{\text{tot}}{\omega }_0{\omega }^2{R}_0}{c_{\text{L}}}\end{array}$(14) where ${\delta }_{\text{tot}}{\omega }_0{\omega }^2{R}_0/c_{\text{L}}$ is from the linearisation of ${\dot{p}}_{\text{a}}(t)R/c_{\text{L}}$, and δ_totω₀ is independent on ω₀. δ_tot is a summation of the damping coefficients due to the heat diffusion (δ_th), the re-radiation of sound (δ_rad), the viscosities of liquid (δ_vis) and shell (δ_shell). The summation of δ_vis, δ_th, and δ_shell can be obtained as: (15) ${\delta }_{\text{vis}}+{\delta }_{\text{th}}+{\delta }_{\text{shell}}=\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}}+{{\kappa }_{\text{s}}/R}_0)}{{\rho }_{\text{L}}{\omega }_0{R}_0^2(1+M(R_0))}$(15) and δ_rad of both uncoated and coated bubbles are obtained as (Leighton, 1994; Li et al., 2018): (16) ${\delta }_{\text{rad1}}=\frac{{\omega }_0{R}_0}{c_{\text{L}}}+\frac{({\omega }^2-{\omega }_0^2)R_0}{{\omega }_0{c}_{\text{L}}}=\frac{{\omega }^2{R}_0}{{\omega }_0{c}_{\text{L}}}$(16) where $({\omega }^2-{\omega }_0^2)R_0/{\omega }_0{c}_{\text{L}}$ is from the linearisation of ${\dot{p}}_{\text{a}}(t)R/c_{\text{L}}$, and ${\omega }_0{R}_0/c_{\text{L}}$ is from the linearisation of the other radiation terms on the right side of the KME. If ${\dot{p}}_{\text{a}}(t)R/c_{\text{L}}$ in the KME is neglected, another δ_rad can be given as (Helfield, 2019): (17) ${\delta }_{\text{rad2}}=\frac{{\omega }_0{R}_0}{c_{\text{L}}}$(17)

If the differential of the surface tension is further neglected, i.e. the linearisation of the simplified KME (Eq. (6) or (7)), one more expression of δ_rad can be obtained (Paul et al., 2010): (18) ${\delta }_{\text{rad3}}=\frac{3\kappa p_{\text{G}0}}{{\rho }_{\text{L}}{\omega }_0{c}_{\text{L}}{R}_0}$(18) where p_G0 was replaced by p₀ in Refs. (Segers et al., 2016; Segers et al., 2018; Versluis et al., 2020).

2.3. Theoretical attenuation coefficient

2.3.1. Linear analysis

As mentioned above, the attenuation of a plane wave travelling through bubbles is due to the liquid viscosity, heat diffusion, and scattering. The scattered pressure at the bubble wall can be obtained as (Segers et al., 2016; Versluis et al., 2020): (19) $p_{\text{s}}(t)={\rho }_{\text{L}}(R\ddot{R}+2{\dot{R}}^2)$(19) The energy loss due to the scattering of a single bubble is usually expressed by the scattering cross-section σ_s, which is defined as the ratio of the scattering energy to the acoustic intensity. For the small-amplitude oscillation, the linearised σ_s is given as (Segers et al., 2016; Versluis et al., 2020): (20) ${\sigma }_{\text{s}}=4\pi R_0^2\frac{{|p_{\text{s}}|}^2}{p_{\text{A}}^2}$(20)

The amplitude of the scattered pressure $|p_{\text{s}}|$ can be approximately obtained as (Helfield, 2019): (21) $|p_{\text{s}}|\approx {\rho }_{\text{L}}{\omega }^2{R}_0^2|x|$(21) where the relative amplitude $|x|$ can be obtained by solving the analytic solution of the linearised equation of motion (Eq. (12)). Substituting Eq. (21) into Eq. (20), the linearised σ_s can be rewritten as (Versluis et al., 2020): (22) ${\mathrm{\Omega }}_{\text{s}}=\frac{4\pi R_0^2}{{\left(\frac{{\omega }_0^2}{{\omega }^2}-1\right)}^2+{\left(\frac{{\delta }_{\text{tot}}{\omega }_0}{\omega }\right)}^2}$(22)

The attenuation by a single bubble can be obtained by multiplying Ω_s by δ_tot/δ_rad. Neglecting the inter-bubble interaction, α_theo of the small-amplitude wave travelling through bubbles in unit distance can be obtained by the summation of the attenuation caused by the bubbles in the acoustic beam (Versluis et al., 2020): (23) ${\alpha }_{\text{lin1}}=\frac{10}{\ln 10}\sum _{R_0}\left(n(R_0){\mathrm{\Omega }}_{\text{s}}\frac{{\delta }_{\text{tot}}}{{\delta }_{\text{rad}}}\right)$(23) where n(R₀) is the bubble number density.

In addition, the energy loss rate of an uncoated bubble also can be obtained by integrating the work of all damping forces (Jamshidi & Brenner, 2013; Leighton, 1994): (24) $\begin{array}{rl}{\dot{W}}_{\text{disp1}}& =f{\int }_0^{1/f}[(\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}})}{R}\dot{R}+\frac{3\kappa p_{\text{G}}-2\sigma /R}{c_{\text{L}}}\dot{R}\\ & +\frac{{\dot{p}}_{\text{a}}(t)R}{c_{\text{L}}})4\pi R^2\dot{R}dt]\\ & =2\pi {\rho }_{\text{L}}{\delta }_{\text{tot}}(1+M(R_0)){\omega }^2{R}_0^5|x{|}^2\end{array}$(24) Then α_theo of the small-amplitude oscillations also can be obtained by the ratio of ${\dot{W}}_{\text{disp1}}$ and sound intensity: (25) $\begin{array}{rl}{\alpha }_{\text{lin2}}& =\frac{10}{\ln 10}\sum _{R_0}\left(\frac{n(R_0){\dot{W}}_{\text{disp1}}}{p_{\text{A}}^2/2{\rho }_{\text{L}}{c}_{\text{L}}}\right)\\ & =\frac{10}{\ln 10}\sum _{R_0}\left(\frac{n(R_0){\mathrm{\Omega }}_{\text{s}}{\delta }_{\text{tot}}{\omega }_0{c}_{\text{L}}}{R_0{\omega }^2}\right)\end{array}$(25)

The above result is also applicable for coated bubbles. It can be seen that α_lin1 = α_lin2 when δ_rad = δ_rad1.

2.3.2. Nonlinear analysis

With the increase of p_A, the bubble oscillations gradually exhibit nonlinear behaviours, which make α_lin1 and α_lin2 inapplicable. In order to predict α_theo considering the nonlinearity, Segers et al. (2016) solved the simplified KME numerically using the ode45 solver in MATLAB, and then calculated $|p_{\text{s}}|$ and σ_s, respectively, by Eqs. (19) and (20). Although the damping coefficients are dependent on p_A, Segers et al. supposed that δ_tot/δ_rad was independent on p_A, and α_nonlin was obtained as (Segers et al., 2016): (26) ${\alpha }_{\text{nonlin1}}=\frac{10}{\ln 10}\sum _{R_0}\left(n(R_0){\sigma }_{\text{s}}\frac{{\delta }_{\text{tot}}}{{\delta }_{\text{rad}}}\right)$(26) It should be noted that Eq. (20) is simplified for small-amplitude oscillations (R ≈ R₀). If the relative amplitude $|x|$ is considerable, the fluctuation of the bubble surface area needs to be considered in Eq. (20). Moreover, the supposition of δ_tot/δ_rad independent on p_A is in doubt.

After solving the bubble dynamics equation numerically, as mentioned in Section 2.3.1, the energy loss rate can be obtained by integrating the work of all damping forces (Leighton, 1994): (27) $\begin{array}{rl}{\dot{W}}_{\text{disp2}}& =f{\int }_0^{1/f}(\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}})}{R}\dot{R}+\frac{3\kappa p_{\text{G}}-2\sigma /R}{c_{\text{L}}}\dot{R}\\ & +\frac{{\dot{p}}_{\text{a}}(t)R}{c_{\text{L}}}-\left(\frac{2\dot{R}}{c_{\text{L}}}-M\right)R\ddot{R}\\ & -(\frac{2\dot{R}}{c_{\text{L}}}+M){\dot{R}}^2)4\pi R^2\dot{R}dt\end{array}$(27) The energy loss of a coated bubble can be similarly obtained. Then α_nonlin can be obtained by: (28) ${\alpha }_{\text{nonlin2}}=\frac{10}{\ln 10}\sum _{R_0}\left(\frac{n(R_0){\dot{W}}_{\text{disp2}}}{p_{\text{A}}^2/2{\rho }_{\text{L}}{c}_{\text{L}}}\right)$(28)

In addition, the energy loss should equal the energy input to the bubble from the acoustic excitation, thus ${\dot{W}}_{\text{disp}}$ also can be obtained as (Leighton, 1994): (29) ${\dot{W}}_{\text{disp3}}=f{\int }_0^{1/f}[-p_{\text{a}}(t)4\pi R^2\dot{R}dt]$(29) Then α_nonlin also can be obtained by: (30) $\begin{array}{rl}{\alpha }_{\text{nonlin3}}& =\frac{10}{\ln 10}\sum _{R_0}\left(\frac{n(R_0){\dot{W}}_{\text{disp3}}}{p_{\text{A}}^2/2{\rho }_{\text{L}}{c}_{\text{L}}}\right)\\ & =\frac{10}{\ln 10}\sum _{R_0}(-\frac{n(R_0)\overline{p_{\text{a}}(t)4\pi R^2\dot{R}}}{p_{\text{A}}^2/2{\rho }_{\text{L}}{c}_{\text{L}}})\end{array}$(30) where the overbar denotes the time average value over one cycle.

2.3.3. Transient analysis

The above analyzes are available for bubble oscillations in the steady (periodic) state. During experimental measurements, a dozen cycles are transmitted and then recorded, and α of these cycles is different from that in the steady state. Usually, p_A of these cycles gradually increases at the beginning and decreases in the final stage (as shown in Figure 4). At the beginning, the energy absorbed by the bubble is not only dissipated but also converted to mechanical energy. In the final stage, the mechanical energy stored by the bubble is released. If the energy loss is calculated by integrating the work of all damping forces (Segers et al., 2016), for example by Eq. (27), the mechanical energy is neglected. According to our simulations, the energy loss should be calculated by Eq. (29) to include the mechanical energy. The average sound intensity of these cycles before the attenuation is ${\int }_0^t{p}_{\text{a}}{(t)}^{2}dt/t{\rho }_{\text{L}}{c}_{\text{L}}$. Thus, α_theo of these cycles can be obtained as: (31) ${\alpha }_{\text{tran}}=\frac{10}{\ln 10}\sum _{R_0}\left(\frac{-n(R_0){\int }_0^t{p}_{\text{a}}(t)4\pi R^2\dot{R}dt}{{\int }_0^t{p}_{\text{a}}{(t)}^{2}dt/{\rho }_{\text{L}}{c}_{\text{L}}}\right)$(31)

2.3.4. Multibubble analysis

The obtained α_theo in Sections 2.3.1–2.3.3 are only applicable at extremely small β. At medium β, the wave number in the mixture of a monodisperse bubble population deduced by Commander and Prosperetti (1989) was given as: (32) $k_{\text{m}}^2=\frac{{\omega }^2}{c_{\text{L}}^2}+\frac{4\pi {\omega }^{2}n(R_0)R_0}{{\omega }_0^2-{\omega }^2+{\delta }_{\text{tot}}{\omega }_0\omega i}\frac{1}{1+M(R_0)}$(32) where M(R₀) was neglected in all previous researches. α_m and V_m in the mixture were respectively given by (Commander & Prosperetti, 1989; Trujillo, 2018): (33) $\begin{array}{rl}{\alpha }_{\text{m}}& =20\log e\cdot \mathrm{Im}(-k_{\text{m}})\end{array}$(33) (34) $\begin{array}{rl}{V}_{\text{m}}& =\frac{\omega }{Re(k_{\text{m}})}\end{array}$(34) where Im and Re are respectively the imaginary and real parts.

On the other hand, for a bubble cluster including several or dozens of bubbles, the bubble dynamics also can be predicted considering the inter-bubble interaction (Ida et al., 2007; Yasui et al., 2008; Yasui et al., 2009). For the bubble oscillation under far-field acoustic excitation, the scattered pressure at a distance of r from the bubble centre was usually obtained by (Ida et al., 2007; Yasui et al., 2008): (35) $p_{\text{s}}=p_r-p_{\mathrm{\infty }}=\frac{{\rho }_{\text{L}}(R^2\ddot{R}+2R{\dot{R}}^2)}{r}$(35) where the liquid compressibility was neglected. Assuming that the adjacent bubbles had the same dynamics, the inter-bubble interaction was usually considered by adding the scattered pressure caused by the surrounding bubbles to the right side of RP type equations (Ida et al., 2007; Yasui et al., 2008; Yasui et al., 2009): (36) $-\left(\sum _i\frac{R}{r_i}\right)(R\ddot{R}+2{\dot{R}}^2)$(36) where r_i is the inter-bubble distance. Adding Eq. (36) to the right side of Eq. (4), ω₀ can be obtained after the linearisation: (37) ${\omega }_0^2=\frac{3\kappa p_0+2(3\kappa -1)\sigma /R_0}{{\rho }_{\text{L}}{R}_0^2\left(1+M(R_0)+R_0\sum r_i^{-1}\right)}+\frac{{\delta }_{\text{tot}}{\omega }_0{\omega }^{2}R}{c_{\text{L}}}$(37) which decreases with the bubble concentration. However, according to the experimental investigation conducted by Leroy et al. (2009), ω₀ increased with the bubble concentration. In our opinion, the main reason might be that the acoustic excitation is supposed to be acting at far field in Eq. (4), and the liquid compressibility is neglected in Eq. (35).

2.4. Experimental attenuation coefficient

The experimental α can be measured by two aligned transducers, where one transmits acoustic waves (usually a dozen cycles) and the other records the pressure (Li et al., 2018; Segers et al., 2016), as shown in Figure 1. α can be obtained by comparing the recorded pressure with and without bubbles in the acoustic beam (Segers et al., 2016): (38) ${\alpha }_{exp}=\frac{10}{d}{\log }_{10}\frac{{|p_{\text{ref}}|}^2}{{|p_{\text{bub}}|}^2}$(38) where $|p_{\text{bub}}{|}^2$ and $|p_{\text{ref}}{|}^2$ are respectively the amplitudes of the power spectra of the recorded pressure with and without bubbles, and d is the acoustic path length over which the wave interacts with the bubbles.

Figure 1. Schematic diagram of the measurements of α.

2.5. Numerical modelling

2.5.1. Numerical measurements of α and V

The numerical measurements were modelled from the experiments described in Section 2.4. Figure 2 shows the 2-dimensional computational grids and boundary conditions for the numerical measurements of α and V. The computational domain was 30λ (λ is the wavelength) in both length and height. The line at the bottom was the axis, and the acoustic excitation was applied on the left boundary. The non-reflecting pressure outlet was specified at the top and right boundaries. Bubbles distributed in the region marked red, whose left border was 8λ away from the left boundary. A monitoring line was created to record the pressure with and without bubbles in the acoustic beam. It should be noted that although the non-reflecting pressure outlet was applied, there was still pressure reflecting back from the right boundary at an amplitude of nearly 1%. Thus, the monitoring line was 15λ away from the right boundary, and only the wave before the 30th cycle was recorded to exclude the interference of reflecting pressure. In addition, the scattered pressure from bubbles was considerable at medium β, it reflected from the left boundary and interfered with the recorded pressure at the monitoring line. Therefore, the distance between bubbles and the left boundary was increased to 15λ at medium β (Section 3.4). Based on the recorded pressure, α_num was obtained similar to α_exp using Eq. (38), i.e. α_num = α_exp. Only V of steady oscillations at medium β were investigated in this study. The recorded pressure (reaching the steady state) with and without bubbles were respectively fitted by sine functions to obtain the time lags t_bub and t_ref, then V_num was obtained by: (39) $\frac{1}{V_{\text{num}}}=\frac{1}{c_{\text{L}}}-\frac{t_{\text{ref}}-t_{\text{bub}}}{d}$(39) The acoustic excitation was a sine wave multiplied by two coefficients C_y and C_p. C_y was dependent on the ratio of y coordinate to λ, as shown in Figure 3(a). Figure 3(b) shows the pressure distribution at p_A = 10 Pa after the acoustic wave traveling through the whole computational domain. It can be seen that the acoustic excitation concentrated around the axis, so that only these bubbles around the axis were irradiated, and the influence of scattered pressure from these bubbles on the recorded pressure was negligible. C_p was used to generate a 20-cycle or Gaussian wave when analysing α in the transient state, as shown in Figure 4. After the grid and time-step independences verification, 242 thousand quadrilateral grids were employed, and the simulation time-step size was 0.004/f. The grid was concentrated in the direction of wave propagation, and the grid number per wavelength was 50.

Figure 2. Computational domain and boundary conditions for the numerical ‘measurements’ of α and V, which were modelled from the experimental schemes. The acoustic excitation was applied on the left boundary. The area-weighted average pressure at the monitoring line were recorded with and without bubbles in the acoustic beam, based on which α_num and V_num were obtained.

Figure 3. (a) Dependence of C_y on y/λ, and (b) the pressure distribution at p_A = 10 Pa. The plane wave concentrated around the axis to reduce the scattered pressure from bubbles far from the axis.

Figure 4. Waveforms of the 20-cycle and Gaussian pulses at p_A = 10 Pa, which were used to analyse α in the transient state.

In numerical simulations, the far-field p_a(t) in the RP equation is unknown, and it needs to be replaced by the local pressure. Thus the bounded RP equation neglecting the liquid compressibility was used in Refs. (Ye et al., 2020b; Ye et al., 2021a), and the local pressure was supposed to act at the distance of R_e away from the bubble. R_e was the radius of the region occupied by the bubble, which means R/R_e = β^1/3. Considering the liquid compressibility, the bounded KME for the numerical simulations of uncoated bubbles was obtained as: (40) $\begin{array}{rl}& \left(1-\frac{R}{R_e}-\frac{2\dot{R}}{c_{\text{L}}}+M\right)R\ddot{R}\\ & +\left(\frac{3}{2}-\frac{2R}{R_e}-\frac{2\dot{R}}{c_{\text{L}}}-M\right){\dot{R}}^2\\ & =\frac{1}{{\rho }_{\text{L}}}[p_{\text{G}}-\frac{2\sigma }{R}-\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}})\dot{R}}{R}\\ & -\left(3\kappa p_{\text{G}}-\frac{2\sigma }{R}\right)\frac{\dot{R}}{c_{\text{L}}}-\frac{\dot{p}R}{c_{\text{L}}}-p]\end{array}$(40) where p is the local pressure. The bounded KME reduces to the KME when β₀ tends to 0. The equation for coated bubbles can be similarly obtained by adding R/R_e to Eq. (10).

The present homogeneous method was implemented in the commercial software ANSYS FLUENT. The density and viscosity of the mixture were defined as the volume-weighted average values of the two phases. The slip velocity between phases was neglected. The bubble wall velocity $\dot{R}$ was recorded by an added transport equation, whose source term was equal to $\ddot{R}$ calculated by Eq. (40). After obtaining the bubble dynamics, the change rate of β was given as (Ye et al., 2020b): (41) $\frac{d\beta }{dt}=4\pi n{R}^2\dot{R}$(41) The liquid compressibility was modelled by the Tait's equation of state for water (Lauer et al., 2012): (42) ${\left(\frac{{\rho }_{\text{L}}}{{\rho }_{\text{L0}}}\right)}^{7.15}=1+\frac{7.15(p-p_0)}{K_0}$(42) where ρ_L0= 1000 kg/m³ and K₀= 2.2 × 10⁹ Pa. The other physical parameters for microbubbles were as follows: p₀ = 100 kPa, ρ_G= 1 kg/m³, μ_L= µ_th = 0.001 Pa·s, σ = 0.072 N/m, χ = 0.5 N/m, κ_s = 10⁻⁸ kg/s, R_b = 0.97R₀ (σ(R₀) = 0.0314 N/m), R_ruptured = 1.2R₀, κ = 1.4 for uncoated bubbles and 1.07 for coated bubbles. All flows were considered to be laminar and the gravity was neglected.

2.5.2. Validation at medium β

In addition to the numerical measurements described in Section 2.5.1, comparisons were made with the predictions by the interface capturing method VOF (also using ANSYS FLUENT). Figure 5 shows the computational grids and boundary conditions of the validation cases. Acoustic excitation was applied on the left boundary and traveled through a pipe. 49 bubbles were distributed around the outlet of the pipe, and they were arranged neatly (7 × 7) on the cross-section parallel to the left boundary at a distance of 8λ. The inter-bubble distance was 10⁻⁴ m, which means n = 10¹² m⁻³. A monitoring face was created at a distance of 2λ from bubbles to obtain α_num, and only the wave before the 16th cycle was recorded. Two symmetry planes were applied as shown in Figure 6, and only one-eighth of the bubbles were inside the computational domain. For the VOF method, special contractions of the mesh were applied around bubble surfaces, about 1 million tetrahedral grids were employed around the bubbles, and 0.5 million hexahedral grids were employed in the other region; the gas inside the bubbles obeyed the ideal-gas law. For the present method, the grid number around bubbles was decreased to 0.1 million, while the grids in the other region were kept the same. It is complex for the VOF method to handle the surface tension of microbubbles, especially for coated microbubbles. It is also complex for the homogeneous models to precisely predict the heat transfer. For a better comparison, the surface tension was neglected, and the flow was assumed to be isothermal, which is the same as in Refs. (Ye et al., 2020b; Ye et al., 2021a). These two simplifications were only used in this case.

Figure 5. Computational grids and boundary conditions for the measurement of α through 49 bubbles when using the VOF method.

Figure 6. (a) Distribution of bubbles on the cross-section when using the VOF method. (b) Special contractions of the mesh were applied around bubble surfaces.

In order to further validate the present method at medium β, the experimental measurements of α and V conducted by Silberman (1957) were simulated at β₀ = 3.77 × 10⁻⁴. Standing waves were established inside a vertical tube filled with air bubbles in the experiments. Since the pressure can be precisely obtained during numerical simulations, there is no need to generate standing waves, and the computational domain shown in Figure 2 was still used. Both R₀ and β₀ were medium in these cases, it took a long time for bubbles to reach the steady state, especially around resonance. For simplicity, bubbles were uniformly distributed in the whole computational domain, just as the sound generator and hydrophone were placed inside the bubbly liquid in the experiments (Silberman, 1957). And two monitoring points (points 1 and 2) were created close to the left boundary to record the pressure, as sketched in Figure 2. α_num was obtained by comparing the power spectra of the pressure at these two points. The pressure at these two points were also fitted by sine functions to obtain the time lags t₁ and t₂, V_num of these cases were obtained by: (43) $V_{\text{num}}=\frac{d_{12}}{t_2-t_1}$(43) The thermal damping was significantly dependent on f in these cases, µ_th and κ were obtained by (Prosperetti et al., 1988): (44) $\begin{array}{rl}{\mu }_{\text{th}}& =\frac{p_{\text{G0}}}{4\omega }Re\mathrm{\Phi }\end{array}$(44) (45) $\begin{array}{rl}\kappa & =\mathrm{Im}\mathrm{\Phi }/3\end{array}$(45) (46) $\begin{array}{rl}\mathrm{\Phi }& =\frac{3\gamma }{1-3(\gamma -1)iX[{(i/X)}^{1/2}\coth {(i/X)}^{1/2}-1]}\end{array}$(46) (47) $\begin{array}{rl}X& =\frac{D_{\text{G}}}{\omega R_0^2}\end{array}$(47) where γ = 1.4 is the ratio of specific heats, D_G = 2.1 × 10⁻⁵ m²/s is the thermal diffusivity of the gas.

3. Results and discussion

The acoustic attenuation caused by uncoated and coated bubbles were numerically ‘measured’ at p_A = 10 Pa in Section 3.1, and comparisons were made with the theoretical α_lin1 and α_lin2. Each α_lin was respectively predicted using δ_rad1 (α_lin1 = α_lin2), δ_rad2 (neglecting ${\dot{p}}_{\text{a}}(t)$), and δ_rad3 (neglecting ${\dot{p}}_{\text{a}}(t)$ and the differential of surface tension). Then p_A was increased to 10 and 20 kPa in Section 3.2, α_num were compared with α_nonlin1 (Segers et al., 2016), α_nonlin2 (work of damping forces), and α_nonlin3 (work of p_a(t)). In Section 3.3, the 20-cycle and Gaussian pulses were applied, and α_num of these pulses were compared with α_tran. α_nonlin and α_tran were obtained by solving the KME using the ode solver in MATLAB. For better comparisons with α_theo, n^2/3Ω_s was set to 10⁻⁵ in Sections 3.1-3.3, n values were mostly between 10⁴ and 10⁵ m⁻³ at f = 1 MHz (β₀ were around 10⁻¹¹). In Section 3.4, α and V at medium β₀ (the maximum n^2/3Ω_s exceeded 0.1) were numerically measured, and comparisons were made with the experimental (Silberman, 1957) and theoretical results.

3.1. Small-amplitude oscillations in the steady state

Figure 7 compares α_lin and α_num caused by uncoated bubbles at p_A= 10 Pa, f = 1 and 3 MHz. Firstly, the difference of α_lin1 predicted using δ_rad1, δ_rad2, and δ_rad3 is obvious. Especially when the differential of the surface tension is neglected, α_lin1 is obviously under-estimated around resonance, since δ_rad3 is overpredicted. Secondly, the deviations of α_lin2 predicted using δ_rad2 and δ_rad3 are much smaller. When δ_rad2 is used, deviations mainly appear for bubbles far from resonance. When δ_rad3 is used, the deviations increase for bubbles around resonance, taking f = 1 MHz for example, the relative deviation increases to 1% at resonance. Thirdly, α_num match precisely with α_lin2, both by using the KME (comparing with α_lin1-δ_rad1) and the simplified KME (comparing with α_lin2-δ_rad3). It can be concluded that α_lin1 is only applicable when δ_rad1 is used, and if the other δ_rad is used, the total energy loss should be directly obtained by integrating the work of all damping forces, not by multiplying the scattering energy by δ_tot/δ_rad. In other words, α_lin2 should be used if the differential of p_a(t) or surface tension is neglected.

Figure 7. Comparison of α_lin and α_num caused by uncoated bubbles at p_A= 10 Pa, f = (a) 1 and (b) 3 MHz. α_num predicted using both the KME and simplified KME match precisely with the corresponding α_lin2.

Figure 8. Comparison of α_lin and α_num caused by coated bubbles at p_A= 10 Pa, f = (a) 1 MHz and (b) 3 MHz. α_lin1 linearised from the simplified KME (δ_rad3) is obviously overpredicted, since the differential of the surface tension in the KME cannot be neglected for coated bubbles.

Figure 8 compares α_lin and α_num caused by coated bubbles at p_A= 10 Pa, f = 1 and 3 MHz. α_lin1 obtained using δ_rad3 is overpredicted by more than 100%, the reason is that the differential of the surface tension is dominant in δ_rad for coated bubbles. α_lin2 obtained using δ_rad3 is also obviously overpredicted around resonance, for example, it is overpredicted by 7.6% at resonance at f = 1 MHz (much larger than uncoated bubbles). α_num predicted using the two equations still match precisely with α_lin2. It can be concluded that the differential of the surface tension cannot be neglected for coated bubbles, and the proposed numerical measurements of α match precisely with α_theo for small-amplitude oscillations. In the following sections, α_lin2 obtained using δ_rad1 is used as a reference.

3.2. Nonlinear oscillations in the steady state

Figure 9 compares α_nonlin and α_num caused by uncoated bubbles at f = 1 MHz, p_A= 10 and 20 kPa. ω₀ decreases with increasing p_A (Helfield, 2019). α_num match precisely with α_nonlin2 and α_nonlin3, while α_nonlin1 (using δ_rad1) is obviously overpredicted. There are two reasons leading to the deviation of α_nonlin1. Firstly, σ_s (Eq. (20)) is inapplicable for nonlinear oscillations, since the bubble surface area cannot be treated as a constant. Secondly, the ratio of the total energy loss to scattering energy is dependent on p_A. To prove this, the attenuation due to radiation (α_rad) is obtained by only keeping the radiation damping force in α_nonlin2, and the total attenuation is obtained by multiplying α_rad by δ_tot/δ_rad. It can be seen from Figure 9 that although α_rad is precisely predicted, the attenuation is still overpredicted, which means that the ratio of the total energy loss to scattering energy is overpredicted by δ_tot/δ_rad.

Figure 9. Comparison of α_nonlin and α_num caused by uncoated bubbles at f = 1 MHz, p_A= (a) 10 kPa and (b) 20 kPa. α_num match precisely with α_nonlin2 (work of damping forces) and α_nonlin3 (work of p_a(t)).

Figure 10 compares α_nonlin and α_num caused by coated bubbles at f = 1 MHz, p_A= 10 and 20 kPa. α_num still match precisely with α_nonlin2 and α_nonlin3. α_nonlin1 and ${\alpha }_{\text{rad}}\cdot {\delta }_{\text{tot}}/{\delta }_{\text{rad}}$ are respectively smaller and larger than α_num, which means that the scattering cross-section is under-estimated by Eq. (20), and the ratio of the total energy loss to scattering energy is overpredicted by δ_tot/δ_rad.

Figure 10. Comparison of α_nonlin and α_num caused by coated bubbles at f = 1 MHz, p_A= (a) 10 and (b) 20 kPa.

3.3. Attenuation in the transient state

Figure 11 compares α_tran and α_num of the 20-cycle and Gaussian pulses through uncoated and coated bubbles at f = 1 MHz, p_A= 10 Pa. Firstly, both α_tran and α_num are smaller than α_lin2 around resonance while larger than α_lin2 far from resonance. The reason for the former is that the amplitude of R gradually increases for the oscillation around resonance. As shown in Figure 12, although the pressure amplitudes of cycles 4–10 have reached p_A, their relative amplitudes $|x|$ are smaller than that in the steady state, which reduces α. On the other hand, it takes much shorter time to reach the steady state for bubbles far from resonance, and some absorbed energy is converted to mechanical energy, which increases α. α_num predicted using cycles 14–17 agree much better with α_lin2 (see Figure 11), since these four cycles are close to the steady state. The Gaussian pulse only consists of several cycles as shown in Figure 11, larger deviations exist between α_num and α_lin2, the relative difference was 45% at R₀ = 3.72 µm. Secondly, α_num of the 20-cycle pulse does not match well with α_tran, especially around resonance, and the relative difference was 4.2% at R₀ = 3.72 µm. The reason is that the formula of α_num (Eq. (38)) is applicable for oscillations in the steady state. For transient oscillations, the sound intensities of these 20 cycles with and without bubbles can be directly obtained by integrations, then a new formula can be obtained: (48) ${\alpha }_{\text{numt}}={\alpha }_{exp\text{t}}=\frac{10}{d}{\log }_{10}\frac{{\int }_0^t{p}_{\text{ref}}{(t)}^{2}dt}{{\int }_0^t{p}_{\text{bub}}{(t)}^{2}dt}$(48) It can be seen from Figure 11 that α_numt match well with α_tran, and the relative difference decreased to 0.59% at R₀ = 3.72 µm. Thus the experimental α of transient oscillations should also be predicted using Eq. (48) instead of Eq. (38). Figure 13 compares α_tran and α_num of the 20-cycle pulse when p_A is increased from 10 Pa to 20 kPa, α_tran around resonance is still smaller than that in the steady state (α_nonlin3), and α_numt still match well with α_tran.

Figure 11. α_tran and α_num of the 20-cycle and Gaussian pulses through (a) uncoated and (b) coated bubbles at f = 1 MHz, p_A= 10 Pa. α_num is smaller than α_lin2 around resonance while larger than α_lin2 far from resonance, especially for the Gaussian pulse.

Figure 12. Evolution of x of the uncoated bubble with R₀ = 3.72 µm at f = 1 MHz, p_A= 10 Pa predicted using the ode solver in MATLAB.

Figure 13. α_tran and α_num of the 20-cycle pulse through (a) uncoated and (b) coated bubbles at f = 1 MHz, p_A= 20 kPa.

3.4. Attenuation and phase velocity at medium β

It has been shown that α_num match precisely with α_theo for linear, nonlinear, and transient oscillations at extremely small β. As described in Section 2.5.2, the present method was validated by compared with the results of the VOF method at medium β. Figure 14 compares α_num caused by 49 bubbles predicted by the VOF and present methods at f = 1 MHz, p_A= 10 Pa, n = 10¹² m⁻³ (β₀ around 10⁻⁴). It can be seen that good agreements are obtained between the VOF method and the present method with the bounded KME, and the relative differences are all less than 1.6%. Obvious deviations are obtained when replacing the bounded KME with the KME (the maximum relative difference reaches 15.4%), and R₀ corresponding to the peak of α decreases. Since the VOF method can precisely predict the dynamics of multibubbles, α was well predicted by the present method at β₀ around 10⁻⁴.

Figure 14. Comparison of α_num caused by 49 bubbles predicted by the VOF and present methods at f = 1 MHz, p_A= 10 Pa, n = 10¹² m⁻³. The relative differences between the VOF method and the present method (with the bounded KME) are all less than 1.6%, while obvious deviations exist using the present method with the KME.

After comparing with the VOF method, comparisons with the experimental and theoretical results were made. Figure 15 compares the experimental (Silberman, 1957), theoretical and numerical α and V through air bubbles of R₀ = 0.994 and 1.02 mm at β₀ = 3.77 × 10⁻⁴. It was found that α_num and V_num predicted using the KME matched precisely with α_m and V_m. As concluded in the previous paragraph, the bounded KME should be used in numerical simulations, thus we deduced the theoretical α and V based on the bounded KME (α_mb and V_mb). The linearisation of the bounded KME lead to: (49) $\begin{array}{rl}& {\omega }_{0\text{b}}^2=\frac{3\kappa p_0+2(3\kappa -1)\sigma /R_0}{{\rho }_{\text{L}}{R}_0^2(1+M(R_0)-{\beta }_0^{1/3})}+\frac{{\delta }_{\text{tot}}{\omega }_0{\omega }^{2}R}{c_{\text{L}}}\end{array}$(49) (50) $\begin{array}{rl}& {\delta }_{\text{vis,b}}+{\delta }_{\text{th,b}}=\frac{4({\mu }_{\text{L}}+{\mu }_{\text{th}})}{{\rho }_{\text{L}}{\omega }_0{R}_0^2(1+M(R_0)-{\beta }_0^{1/3})}\end{array}$(50)

Figure 15. Comparisons of the experimental (Silberman, 1957), theoretical and numerical (a) α and (b) V caused by uncoated bubbles with R₀ = 0.994 and 1.02 mm at β₀ = 3.77 × 10⁻⁴. α_num predicted using the KME match precisely with α_m, while those using the bounded KME match precisely with the proposed α_mb.

The corresponding wave number for the monodisperse bubble population was obtained as: (51) $k_{\text{mb}}^2=\frac{{\omega }^2}{c_{\text{L}}^2}+\frac{4\pi {\omega }^{2}n(R_0)R_0}{{\omega }_{0\text{b}}^2-{\omega }^2+{\delta }_{\text{tot,b}}{\omega }_{0\text{b}}\omega i}\frac{1}{1+M(R_0)-{\beta }_0^{1/3}}$(51) α_mb and V_mb based on the bounded KME can be obtained by replacing k_m in Eqs. (33) and (34) with k_mb. It can be seen from Figure 15 that α_num and V_num predicted using the bounded KME match precisely with α_mb and V_mb, the relative differences of α were all less than 1%. It can be concluded that the bounded KME should be used in both numerical and theoretical predictions of α and V. This conclusion is consistent with the actual situation that the acoustic excitation is not at infinity, but near bubbles. f corresponding to the peak of α_lin2, α_m and α_mb are respectively named as: f_lin2, f_m and f_mb. f_m is 1.19% larger than f_lin2, while 3.74% smaller than f_mb since ω₀ is 3.85% smaller than ω_0b, which is consistent with the results in the previous paragraph: R₀ corresponding to the α peak is smaller using the KME. If the inter-bubble interaction is considered using Eq. (36), ω₀ will decrease but not increase. The difference between α_m and α_mb is negligible far below resonance, so is the difference between V_m and V_mb.

Figure 16 and Figure 17 show the influences of β₀ on α/n respectively caused by uncoated and coated microbubbles at p_A= 10 Pa, f = 1 and 3 MHz. R₀ corresponding to the peak of α_m increases with β₀, and it further increases when using the bounded KME (α_mb). α_num predicted using the bounded KME match precisely with α_mb. Figure 18 compares V_m, V_mb and V_num caused by uncoated and coated microbubbles at β₀ = 10⁻⁵, f = 1 MHz, p_A= 10 Pa. R₀ corresponding to the peak of V_mb is larger than that of V_m due to the increase of resonance frequency. V_num predicted using the bounded KME match precisely with V_mb.

Figure 16. Influences of β₀ on α/n caused by uncoated microbubbles at p_A= 10 Pa, f = (a) 1 MHz and (b) 3 MHz. α_num predicted using the bounded KME match precisely with α_mb.

Figure 17. Influences of β₀ on α/n caused by coated bubbles at p_A= 10 Pa, f = (a) 1 MHz and (b) 3 MHz.

Figure 18. Comparisons of V_m, V_mb and V_num caused by (a) uncoated and (b) coated microbubbles at β₀ = 10⁻⁵, f = 1 MHz, p_A= 10 Pa. V_num predicted using the bounded KME match precisely with V_mb.

Figure 19. Influence of d on α_num caused by uncoated bubbles at β₀ = 2 × 10⁻⁶, f = 1 MHz, p_A= 10 Pa.

Figure 20. Influences of d on α_num caused by (a) uncoated and (b) coated bubbles at β₀ = 2 × 10⁻⁶, f = 1 MHz, p_A= 20 kPa.

The advantage of the present CFD method is to deal with the complex conditions combining the nonlinearity, transience, medium β, focused wave and so on, which are difficult to handle by theoretical analyzes. Figure 19 shows the influence of d on α_num at β₀ = 2 × 10⁻⁶, f = 1 MHz, p_A= 10 Pa. It can be seen that the influence of d on the linear α is negligible. For nonlinear oscillations, α is dependent on p_A, therefore, it is also dependent on d at medium β. Figure 20 shows the influences of d on α_num at β₀ = 2 × 10⁻⁶, f = 1 MHz, p_A= 20 kPa. R₀ corresponding to the α peak increases obviously with the increasing d, since p_A decreases during the propagation. Especially for uncoated bubbles (Figure 20a), α_num becomes closer to α_lin2 with the increasing d. Figure 21 shows α_num and α_numt of the 20-cycle pulse at β₀ = 10⁻⁵, and those at small β₀ are repeated from Figure 11a as references. The transient α_num is 12.5% smaller than α_lin2 around resonance (R₀ = 3.72 µm) at small β₀, while it is only about half of α_mb around resonance (R₀ = 3.8 µm) at β₀ = 10⁻⁵. And the differences between α_num and α_numt become greater at β₀ = 10⁻⁵.

Figure 21. α_num and α_numt of the 20-cycle pulse at β₀ = 10⁻⁵, p_A = 10 Pa. The differences between the steady (α_mb) and transient α at medium β are much greater than those at small β. The differences between α_num and α_numt also increase.

4. Conclusions and prospects

Numerical measurements of α and V were implemented using a homogeneous cavitation model, which were modelled from experimental schemes, and comparisons were made with α_theo. At extremely small β, α_num match precisely with α_theo for linear, nonlinear, and transient oscillations. For linear oscillations, α_lin2 was obtained, which calculated the total energy loss by integrating the work of all damping forces, instead of by multiplying the scattering energy by δ_tot/δ_rad. For nonlinear oscillations, both the linearised σ_s (Eq. (20)) and δ_tot/δ_rad are dependent on p_A. For transient oscillations, the total energy loss should be obtained by integrating the work of p_a(t) when calculating the theoretical α (α_tran) to include the mechanical energy, and a new formula was proposed for the calculations of the numerical (α_numt) and experimental (α_expt) α by integrating the sound intensity of recorded pressure.

At medium β, the relative difference of α_num between the VOF and present methods was less than 1.6%, while it reached 15.4% after replacing the bounded KME in the present method with the KME. And simulations towards millimetre bubbles showed that α_num and V_num predicted using the KME matched precisely with the traditional theoretical results (α_m and V_m). Then new theoretical α_mb and V_mb were proposed based on the bounded KME, and the relative differences between α_mb and α_num (based on the bounded KME) were less than 1% at β₀ = 3.77 × 10⁻⁴. Therefore, the bounded KME should be used in both numerical and theoretical predictions, and the resonance frequency increases with β (Eq. (49)). In addition, the influences of nonlinearity and transience on α at medium β are quite different from that at small β, and the present method can be applied to investigate complex acoustic cavitation where theoretical analyzes are difficult to apply.

The bubble radius is usually distributed within a range. This study focused on α caused by equal-sized bubbles. The inter-bubble interaction between different-sized bubbles is more complex and needs to be considered by using the grouping cavitation model in the future.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available within the article.

Additional information

Funding

Financial support from the Zhejiang Provincial Natural Science Foundation of China (No. LQ21E060003, LY23E060001, LY22E080006) and Fundamental Research Funds of Zhejiang University of Science and Technology (No. 2021QN029).