A Model for the Dynamics of Stable Gas Bubbles in Viscoelastic Fluids Written in Bubble Volume Variation

1. Introduction

Acoustic cavitation, the formation, growth, and oscillation of gas bubbles under the influence of an acoustic field, is a fundamental phenomenon in a wide range of scientific, biomedical, and industrial applications. Technologies such as ultrasound imaging, targeted drug delivery [1], sonochemistry [2], and materials processing rely on the controlled behavior of bubbles to achieve both efficacy and safety.

Bubble dynamics has traditionally been studied in Newtonian fluids, such as water, for which there are well-established theoretical and numerical models [3]. However, real-world applications, including cavitation-enhanced therapies, food processing, cavitation erosion, and oceanic microbubble scattering, often involve bubbles oscillating in non-Newtonian media. This has driven a growing interest in understanding the dynamics of bubbles in these fluids, particularly with the advent of new materials, the inclusion of polymeric additives in many fluids [4], and the expansion of biomedical technologies into soft media. These fluids include polymeric liquids, gels, colloidal suspensions, and, notably, viscoelastic biological tissues such as cells and tissues, which exhibit mechanical behavior lying between that of an ideal fluid and a solid. These materials are characterized by viscosity and elasticity. For example, Jamburidze et al. [5] characterized the resonant behavior of ultrasound in isolated microbubbles embedded in agarose gels, commonly used as tissue-mimicking phantoms for biomedical applications. They found that the resonance frequency of the bubbles increased with the shear modulus of the medium. This implies that imaging and therapeutic ultrasound protocols must be optimized on the basis of the type of tissue in which the bubbles are embedded.

Moreover, substituting Newtonian liquids with viscoelastic media, such as polymeric gels, in experimental setups has proven to be an effective approach to overcome challenges such as uniformity of bubble sizes and bubble retention [6]. For example, gas-seeded polymeric gels have been used as analog models for oceanic bubble clouds in geophysical acoustics research [7], providing a more controlled platform for exploring complex acoustic interactions.

The viscoelasticity and its impact on the dynamics of spherical bubble is the main focus of this paper. Several studies have investigated bubble behavior in viscoelastic media, leading to models based on extensions of the Rayleigh–Plesset equation, accounting for viscosity, elasticity, and compressibility [8]. The early work by Fogler & Goddard [9] combined the Maxwell model with the Rayleigh-Plesset equation to study bubble collapse in viscoelastic fluids, finding that elasticity can retard collapse in certain parameter ranges. Tanasawa & Yang [10] applied the three-parameter Oldroyd model to bubble oscillations, showing that elasticity reduces viscous damping compared to pure fluids. Allen & Roy [11] used the Maxwell and Jeffreys models to study bubble oscillations, revealing that elasticity alters the oscillation phase and harmonic structure. Yang & Church [12] extended the Keller-Miksis equation by incorporating the Kelvin-Voigt model, showing that elasticity increases the threshold pressure for inertial bubble oscillation. Hua & Johnsen [13] coupled the standard linear solid model with the compressible Keller-Miksis equation, and found that nonlinear effects are significant when the stress tensor of the medium is prominent, affecting the bubble dynamics. Warnez & Johnsen [14] simulated the bubble dynamics in various viscoelastic media, finding that relaxation time increases bubble growth. Both Zilonova et al. [15] and Filonets et al. [16] employed the Gilmore-Akulichev-Zener model to study bubble dynamics in viscoelastic soft tissues. Zilonova et al. focused on high-intensity focal ultrasound thermal therapy, showing that elasticity and viscosity dampen oscillations, with relaxation time affecting bubble behavior depending on ultrasound frequency. Similarly, Filonets et al. studied the inertial cavitation threshold in viscoelastic soft tissues using a dual-frequency signal. Finally, Murakami [17] experimentally observed and theoretically analyzed the oscillations of a spherical bubble in a gelatin gel under ultrasound irradiation. They compared the finite-amplitude oscillations they found with nonlinear Rayleigh–Plesset calculations, suggesting the need to include the gel elasticity in the calculations to accurately reproduce the nonlinear bubble dynamics.

The study of bubble dynamics in viscoelastic media primarily involves integrating different constitutive viscoelastic models into the equations that govern the bubble dynamics. This approach aims to generate models applicable to various tissues and conditions, allowing for accurate predictions of bubble behavior. The correct choice of constitutive equations is crucial for describing the rheological properties of the surrounding medium and capturing the complex interactions between the bubble and its environment. Several linear viscoelastic models, such as the Maxwell, Kelvin-Voigt, and Zener models, are commonly used in the literature to simulate bubble oscillations in viscoelastic materials. They are effective when the deformations within the material are relatively small. However, in many practical scenarios, especially under large-amplitude acoustic forcing, viscoelastic media respond nonlinearly, significantly altering bubble dynamics. Nonlinear viscoelastic models capture these effects and predict complex phenomena such as aperiodic oscillations, modulation, and chaos under intense excitation. [18].

The most popular nonlinear equation to describe the behavior of a gas bubble in an acoustic field is the Rayleigh–Plesset equation [19,20]. Most differential equations developed in this context adopt the bubble radius as the dependent variable (radius framework) [21], and successfully describe the nonlinear oscillations and collapse of the bubble. Alternatively, models based on the bubble volume as the dependent variable also exist. A classic second-order approximation was first derived by Zabolotskaya and Soluyan in the 1960s–1970s [22,23], followed by a third-order approximation in the 1990s [24]. These models describe the bubble behavior under the assumption that volume variations remain small relative to the initial volume of the bubble, making them suitable for moderate-amplitude oscillations. A fourth-order approximation was later introduced within the Rayleigh–Plesset framework [25], enhancing the accuracy of the bubble response under acoustic forcing by including a fourth-order term in the adiabatic gas law in the development of the volume-frame equation. Reformulating the Rayleigh–Plesset equation in terms of bubble volume variation yields a physically equivalent model that retains the same dynamics as the classical radius-based form, offering a complementary perspective for analyzing bubble behavior [26]. Any variations in their predictions can be ascribed to factors such as numerical precision or differences in the truncation of asymptotic expansions, rather than differences in the fundamental physics.

To our knowledge, no previous model has described the dynamics of bubbles in viscoelastic media using a volume-based formulation of the Rayleigh–Plesset equation. Motivated by this gap and the aforementioned studies, we propose a new model that couples the Rayleigh–Plesset equation with the Kelvin–Voigt viscoelastic constitutive law, formulated in terms of bubble volume variation. This approach extends previous works in Newtonian fluids by incorporating the effects of elasticity. The remainder of the paper is organized as follows: Section 2 presents the derivation of the model; in Section 3, we validate the performance of the proposed model, present some numerical simulations, and analyze the impact of the viscoelasticity of the medium on the bubble dynamics; the main conclusions of this work are given in Section 4.

2. Materials and Methods

In this section, we derive a generalized Rayleigh–Plesset equation for bubble dynamics in an arbitrary viscoelastic medium, formulated in terms of bubble volume variation. Building on this framework, we adopt the Kelvin–Voigt constitutive relation and focus on how bubble oscillations are influenced by the linear viscoelastic response of the medium.

2.1. Governing Equations

2.1.1. Bubble Dynamics in an Arbitrary Viscoelastic Medium Written in the Volume-Variation Framework

The Rayleigh–Plesset equation describes the radial dynamics of a spherically symmetric gas bubble in an incompressible fluid and can be expressed in a general form valid for any rheological behavior of the surrounding medium [11]:

$$\rho \left(R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2\right)=P_b-P_{\infty }-{\displaystyle \frac{2\sigma }{R}}+3{\int }_R^{\infty }{\displaystyle \frac{{\tau }_{rr}}{r}}dr,$$

(1)

where $R\left(t\right)$ is the instantaneous radius of the bubble, $\rho$ is the density of the fluid, $\sigma$ is the surface tension, $P_b$ is the gas pressure inside the bubble, $P_{\infty }$ is the ambient pressure, and ${\tau }_{rr}$ is the radial component of the deviatoric stress tensor (in the radial direction r). Note that the overdots denote time derivatives. Assuming an ideal gas inside the bubble, its pressure is modeled using a polytropic relation:

$$P_b=\left(P_{b0}+{\displaystyle \frac{2\sigma }{R}}\right){\left({\displaystyle \frac{R_0}{R}}\right)}^{3\gamma },$$

(2)

where $\gamma$ is the polytropic exponent, $P_{b0}={\rho }_b{c}_b^2/\gamma$ is the gas pressure at equilibrium, and $R_0$ is the radius of equilibrium bubble. Here, ${\rho }_b$ and $c_b$ denote the density and the speed of sound of the gas at equilibrium, respectively. The ambient pressure is given by:

$$P_{\infty }=P_{b0}+p_a\left(t\right),$$

(3)

with $p_a\left(t\right)$ being the acoustic pressure that excites the bubble.

To derive a nonlinear ordinary differential equation governing bubble oscillations in terms of bubble volume variation, we begin with Eq. (1). Neglecting surface tension, the equation becomes:

$$R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2={\displaystyle \frac{1}{\rho }}\left(P_b-P_{\infty }+3{\int }_R^{\infty }{\displaystyle \frac{{\tau }_{rr}}{r}}dr\right).$$

(4)

Let us define the instantaneous radius and volume of the bubble at time t as $R\left(t\right)=R_0+r\left(t\right)$ and $V\left(t\right)=V_0+v\left(t\right)$, where $r\left(t\right)$ and $v\left(t\right)$ represent the radius and volume variations of the bubble, respectively. The equilibrium volume is given by $V_0={\displaystyle \frac{4\pi }{3}}{R}_0^3$. Multiplying both sides of Eq. (4) by $4\pi R_0$, and using the relation for $P_{\infty }$ in Eq. (3), we obtain:

$$4\pi R_0\left(R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2\right)={\displaystyle \frac{4\pi R_0}{\rho }}\left(P_b-P_{b0}\right)-{\displaystyle \frac{4\pi R_0}{\rho }}{p}_a+{\displaystyle \frac{4\pi R_0}{\rho }}3{\int }_R^{\infty }{\displaystyle \frac{{\tau }_{rr}}{r}}dr.$$

(5)

Following the procedure outlined in [25], we expand the right-hand side of Eq. (5) in terms of the bubble volume perturbation $v\left(t\right)$, obtaining a second-order nonlinear approximation, where ${\omega }_b^2=3\gamma P_{b0}/\rho R_0^2$ defines the Minnaert resonance frequency of the bubble in a Newtonian fluid [27], and $a=(\gamma +1){\omega }_b^2/2V_0$ is the second-order nonlinear coefficient (due to the nonlinearity of the adiabatic state equation of gas in the bubble):

$$4\pi R_0\left(R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2\right)=-{\omega }_b^{2}v+a{v}^2-{\displaystyle \frac{4\pi R_0}{\rho }}{p}_a+{\displaystyle \frac{4\pi R_0}{\rho }}3{\int }_R^{\infty }{\displaystyle \frac{{\tau }_{rr}}{r}}dr.$$

(6)

To express the left-hand side of Eq. (4) in terms of the volume variation, we make use of the following relations:

$$V={\displaystyle \frac{4}{3}}\pi R^3,\dot{V}=4\pi R^2\dot{R},\ddot{V}=4\pi R(2{\dot{R}}^2+R\ddot{R}).$$

(7)

Therefore, in terms of V:

$$R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2={\displaystyle \frac{1}{4\pi }}{\left({\displaystyle \frac{4\pi }{3V}}\right)}^{{\displaystyle \frac{1}{3}}}\left(\ddot{V}-{\displaystyle \frac{1}{6}}{\displaystyle \frac{{\dot{V}}^2}{V}}\right).$$

(8)

Since $V=V_0+v$, we have $\dot{V}=\dot{v}$ and $\ddot{V}=\ddot{v}$. On the other hand, providing that $v\ll V_0$ (bubble oscillations assumed to be of moderate amplitude), which implies $V=V_0+v\approx V_0$, we can approximate some terms through the linear exponent Taylor series expansion:

$$V^{-1/3}=V_0^{-1/3}{\left(1+{\displaystyle \frac{v}{V_0}}\right)}^{-1/3}\approx {\displaystyle \frac{1}{V_0^{1/3}}}-{\displaystyle \frac{1}{3V_0^{4/3}}}v.$$

(9)

Assuming the error of approximation of order $O{(v/V_0)}^2$, Eq. (8) can be written:

$$\begin{array}{cc}\hfill R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2& ={\displaystyle \frac{1}{4\pi }}{\left({\displaystyle \frac{4\pi }{3}}\right)}^{{\displaystyle \frac{1}{3}}}\left({\displaystyle \frac{1}{V_0^{1/3}}}-{\displaystyle \frac{1}{3V_0^{4/3}}}v\right)\left(\ddot{v}-{\displaystyle \frac{{\dot{v}}^2}{6(V_0+v)}}\right).\hfill \end{array}$$

(10)

As seen in the previous paragraph, we have $v<<V_0$ and $V=V_0+v\approx V_0$, which allows us to write:

$$\begin{array}{cc}\hfill R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2& \approx {\displaystyle \frac{1}{4\pi }}{\left({\displaystyle \frac{4\pi }{3}}\right)}^{{\displaystyle \frac{1}{3}}}\left({\displaystyle \frac{\ddot{v}}{V_0^{1/3}}}-{\displaystyle \frac{{\dot{v}}^2}{6V_0^{4/3}}}-{\displaystyle \frac{v\ddot{v}}{3V_0^{4/3}}}+{\displaystyle \frac{v{\dot{v}}^2}{18V_0^{7/3}}}\right).\hfill \end{array}$$

(11)

Assuming the error of approximation $V=V_0+v\approx V_0$ and since $V_0=4\pi R_0^3/3$, the above expression becomes,

$$\begin{array}{cc}\hfill R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2& ={\displaystyle \frac{1}{4\pi R_0}}\left(\ddot{v}-{\displaystyle \frac{1}{6V_0}}({\dot{v}}^2+2v\ddot{v})+{\displaystyle \frac{v{\dot{v}}^2}{18V_0^2}}\right),\hfill \end{array}$$

(12)

and we neglect the third-order term to obtain,

$$4\pi R_0\left(R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2\right)=\ddot{v}-{\displaystyle \frac{1}{6V_0}}({\dot{v}}^2+2v\ddot{v}).$$

(13)

Substituting Eq. (13) into Eq. (6) allows us to express Eq. (4) entirely in terms of the volume variation v:

$$\ddot{v}-{\displaystyle \frac{1}{6V_0}}({\dot{v}}^2+2v\ddot{v})=-{\omega }_b^{2}v+a{v}^2-{\displaystyle \frac{4\pi R_0}{\rho }}{p}_a+{\displaystyle \frac{4\pi R_0}{\rho }}3{\int }_R^{\infty }{\displaystyle \frac{{\tau }_{rr}}{r}}dr.$$

(14)

After setting the nonlinear parameter, due to the dynamic nonlinearity of the bubble, $b=1/6V_0$ and the constant $\eta =4\pi R_0/\rho$, we finally obtain the following new Rayleigh-Plesset differential equation of second-order approximation, which governs the nonlinear behavior of a bubble in an acoustic field in terms of bubble volume variation in an arbitrary viscoelastic medium:

$$\ddot{v}+{\omega }_b^{2}v=b({\dot{v}}^2+2v\ddot{v})+a{v}^2-\eta p_a+3\eta {\int }_R^{\infty }{\displaystyle \frac{{\tau }_{rr}}{r}}dr.$$

(15)

2.1.2. Bubble Dynamics in a Kelvin-Voigt Viscoelastic Medium Written in the Volume-Variation Framework

Now consider a homogeneous, isotropic viscoelastic medium. In order to determine the stress term ${\tau }_{rr}$, it is necessary to choose a suitable viscoelastic model, i.e., a constitutive equation that relates the radial stress tensor to the spatial deformations in the fluid. The linear Kelvin–Voigt model [28] is chosen as a starting framework, as it is a simple yet effective model and has been shown to outperform other models, such as the Maxwell model, in capturing soft tissue dynamics [29]:

$${\tau }_{rr}=2G{\gamma }_{rr}+2\mu {\dot{\gamma }}_{rr},$$

(16)

where ${\gamma }_{rr}$ is the strain, ${\dot{\gamma }}_{rr}$ the strain rate in the radial direction, G is the shear modulus, and $\mu$ the shear viscosity of the medium. As $G\to 0$, the elastic contribution disappears, leaving only a viscous stress response. In this case, the classical Rayleigh–Plesset equation for a Newtonian fluid is recovered. Evaluating the integral in the right-hand side of Eq. (1), by taking Eq. (16) into account, yields [12]:

$$\rho \left(R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2\right)=P_b-P_{\infty }-{\displaystyle \frac{2\sigma }{R}}-{\displaystyle \frac{4\mu \dot{R}}{R}}-{\displaystyle \frac{4G}{3}}{\displaystyle \frac{R^3-R_0^3}{R^3}}.$$

(17)

We follow a procedure analogous to that described in Section 2.1.1 to incorporate the additional viscous and elastic contributions into the bubble dynamics equation, formulated in terms of volume variation. Thus, Eq. (17), once multiplied by $4\pi R_0$ and using the expression for the ambient pressure in Eq. (3), can be rewritten as:

$$4\pi R_0\left(R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2+{\displaystyle \frac{4\mu \dot{R}}{\rho R}}+{\displaystyle \frac{4G}{3\rho }}{\displaystyle \frac{R^3-R_0^3}{R^3}}\right)={\displaystyle \frac{4\pi R_0}{\rho }}\left(P_b-P_{b0}\right)-{\displaystyle \frac{4\pi R_0}{\rho }}{p}_a,$$

(18)

The first term in the right-hand side of Eq. (18) can be expressed, following the same procedure used from Eq. (5) to Eq. (6), leading to:

$$4\pi R_0\left(R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2+{\displaystyle \frac{4\mu \dot{R}}{\rho R}}+{\displaystyle \frac{4G}{3\rho }}{\displaystyle \frac{R^3-R_0^3}{R^3}}\right)=-{\omega }_b^{2}v+a{v}^2-{\displaystyle \frac{4\pi R_0}{\rho }}{p}_a.$$

(19)

Using the relations in Eq. (7) and assuming $v\ll V_0$, i.e., $V=V_0+v\approx V_0$, we can rewrite the new left-hand side in terms of V, retaining errors of order $O{(v/V_0)}^2$:

$$\begin{array}{cc}\hfill {\displaystyle \frac{4\mu \dot{R}}{\rho R}}+{\displaystyle \frac{4G}{3\rho }}{\displaystyle \frac{R^3-R_0^3}{R^3}}& ={\displaystyle \frac{4\mu }{3\rho }}{\displaystyle \frac{\dot{V}}{V}}+{\displaystyle \frac{4G}{3\rho }}{\displaystyle \frac{V-V_0}{V}}\hfill \\ \hfill & ={\displaystyle \frac{4\mu }{3\rho }}{\displaystyle \frac{\dot{v}}{V_0+v}}+{\displaystyle \frac{4G}{3\rho }}{\displaystyle \frac{V_0+v-V_0}{V_0+v}}\hfill \\ \hfill & \approx {\displaystyle \frac{4\mu }{3\rho }}{\displaystyle \frac{\dot{v}}{V_0}}+{\displaystyle \frac{4G}{3\rho }}{\displaystyle \frac{v}{V_0}}.\hfill \end{array}$$

(20)

Recalling that $V_0={\displaystyle \frac{4\pi }{3}}{R}_0^3$, we simplify further:

$$\begin{array}{cc}\hfill {\displaystyle \frac{4\mu \dot{R}}{\rho R}}+{\displaystyle \frac{4G}{3\rho }}{\displaystyle \frac{R^3-R_0^3}{R^3}}& ={\displaystyle \frac{\mu }{\pi \rho R_0^3}}\dot{v}+{\displaystyle \frac{G}{\rho }}{\displaystyle \frac{v}{\pi R_0^3}}.\hfill \end{array}$$

(21)

By substituting this along with the previous result in Eq. (13) into the left-hand side of Eq. (18), we obtain:

$$\begin{array}{cc}\hfill & 4\pi R_0\left(R\ddot{R}+{\displaystyle \frac{3}{2}}{\dot{R}}^2+{\displaystyle \frac{4\mu \dot{R}}{\rho R}}+{\displaystyle \frac{4G}{3\rho }}{\displaystyle \frac{R^3-R_0^3}{R^3}}\right)=\hfill \\ & =\ddot{v}-{\displaystyle \frac{1}{6V_0}}({\dot{v}}^2+2v\ddot{v})+{\displaystyle \frac{4\mu }{\rho R_0^2}}\dot{v}+{\displaystyle \frac{4G}{\rho R_0^2}}v.\hfill \end{array}$$

(22)

As a result, Eq. (19) can now be entirely reformulated in terms of the bubble volume variation v as:

$$\ddot{v}-{\displaystyle \frac{1}{6V_0}}({\dot{v}}^2+2v\ddot{v})+{\displaystyle \frac{4\mu }{\rho R_0^2}}\dot{v}+{\displaystyle \frac{4G}{\rho R_0^2}}v=-{\omega }_b^{2}v+a{v}^2-{\displaystyle \frac{4\pi R_0}{\rho }}{p}_a.$$

(23)

We now define the new elastic coefficient:

$$\overline{G}={\displaystyle \frac{4G}{\rho R_0^2}}.$$

(24)

Introducing the viscous damping coefficient $\delta =4\mu /\rho {\omega }_b{R}_0^2$ [25], and, as before, the nonlinear coefficient due to the dynamic nonlinearity of the bubble b and constant $\eta$, we finally obtain a novel second-order approximation Rayleigh–Plesset-type equation governing the nonlinear dynamics of a bubble in a viscoelastic medium modeled using the Kelvin–Voigt formulation:

$$\ddot{v}+\delta {\omega }_b\dot{v}+\left({\omega }_b^2+\overline{G}\right)v=b({\dot{v}}^2+2v\ddot{v})+a{v}^2-\eta p_a.$$

(25)

This model will be solved in Section 2.2 and used in Section 3.

2.2. Numerical Solution of the Bubble Equation

We solve the differential system formed by Eq. (25) and the Cauchy conditions, which state the natural rest of the bubble at the onset, from $t=0$ to the lifetime of the experiment, $t=T$:

$$\left\{\begin{array}{cc}\hfill & \ddot{v}+\delta {\omega }_b\dot{v}+\left({\omega }_b^2+\overline{G}\right)v=b({\dot{v}}^2+2v\ddot{v})+a{v}^2-\eta p_a,0<t<T,\hfill \\ \hfill & v\left(0\right)=0,\hfill \\ \hfill & \dot{v}\left(0\right)=0,\hfill \end{array}\right.$$

(26)

in which the acoustic excitation of the bubble is defined as continuous signal at the driven frequency f, and set at $p_a\left(t\right)=p_{0}sin\left(\omega t\right)$, where $\omega =2\pi f$. The second-order ordinary differential equation is reduced to a first-order system of two coupled ordinary differential equations by setting the new variable $u=\dot{v}$, which yields:

$$\left\{\begin{array}{cc}\hfill & \dot{v}=u,0<t<T,\hfill \\ \hfill & \dot{u}+\delta {\omega }_{b}u+\left({\omega }_b^2+\overline{G}\right)v=b(u^2+2v\dot{u})+a{v}^2-\eta p_a,0<t<T,\hfill \\ \hfill & v\left(0\right)=0,\hfill \\ \hfill & u\left(0\right)=0.\hfill \end{array}\right.$$

(27)

The resulting system is then solved using a fourth-order Runge-Kutta method [30] and numerical solutions are obtained for the instantaneous volume variation of the bubble under ultrasonic excitation. In the following, we use a time step of $\tau =1\times {10}^{-9}$s to ensure the convergence of the numerical solution.

3. Results

This section presents the main results obtained with the proposed model Eq. (25). We begin by validating the formulation through the prediction of the resonance frequency of bubbles in viscoelastic media. We then assess the utility of our model by analyzing the variation in bubble volume for several representative soft materials. Next, we examine the effect of medium elasticity on bubble dynamics, focusing on the influence of the newly introduced elastic term. Finally, we explore the broader impact of viscoelasticity on bubble behavior.

3.1. Validation of the Model

To validate our model Eq. (25), a validation study is conducted by comparing the predicted resonance frequency of the bubble with established theoretical results [28]. Specifically, we verify whether the model accurately reproduces the expected resonance frequency of bubble oscillations under linear conditions. The resonance frequency of gas bubbles in different media is computed as a function of their initial radius. Next, using the proposed model, a frequency sweep is performed to identify the excitation frequency at which the bubble exhibits its maximal response, defined as the frequency that induces the largest volume variation. To this end, we first use two distinct sets of media characterized by different ranges of elasticity and viscosity values to verify that the model accurately captures these variations. Subsequently, we employ a set of representative soft media to validate the model against realistic material properties and behavior.

3.1.1. Effect of Elasticity and Viscosity on Bubble Resonance

As a first step, we validate the model Eq. Eq. (25) against the theoretical expression of the natural frequency [28], given by:

$${\omega }_0^2={\displaystyle \frac{1}{\rho R_0^2}}\left[3\gamma P_{b0}+{\displaystyle \frac{2\sigma }{R_0}}(3\gamma -1)+4G\right].$$

(28)

When surface tension $\sigma$ and elasticity G are neglected, this expression reduces to the classical Minnaert frequency that was used to obtain Eq. (6).

For the given media listed in Table 1, characterized by their shear elasticity G, we solve the differential system (26) considering an ultrasonic perturbation of infinitesimal amplitude $p_0=1\mathrm{mPa}$. For simplicity, biological soft tissues are modeled as viscoelastic liquids [18,31,32,33,34]. In this study, standardized values of $\rho =1000{\mathrm{kg}/\mathrm{m}}^3$ for density and $\mu =0.001\mathrm{Pa}\text{·}\mathrm{s}$ for viscosity are used for all biological tissues. Although this simplified scenario does not fully reflect the actual properties of biological tissues [35], it is used to validate the inclusion of shear elasticity in the proposed model.

Figure 1 shows the variation in the linear resonance frequency $f_0={\omega }_0/2\pi$, as a function of the initial bubble radius $R_0$ for gas bubbles embedded in soft media with increasing shear modulus G (Table 1). Theoretical predictions (solid lines) based on Eq. (28) are compared with the results obtained using the proposed model Eq. Eq. (25) (star symbols).

They are all consistent, confirming that for bubbles with an initial radius greater than $1\mu \mathrm{m}$, an increase in the shear modulus G leads to a rise in the resonance frequency of the bubble, as previously reported in the literature [12,28,31]. For example, for a bubble with $R_0=1\mu \mathrm{m}$, Figure 1 displays the natural frequency $f_0$ for each medium considered, illustrating how the resonance frequency shifts with increasing shear modulus G. This upward trend reflects the role of medium elasticity in constraining the bubble’s radial oscillations, thus increasing the resonance frequency of the system. Previous studies investigated the effects of larger variations in shear modulus, specifically $G=0$, 0.5, 1, and 1.5 MPa, and clearly demonstrated that the natural frequency increases noticeably with increasing rigidity G, due to the increased stiffness of the system imparted by the solid-like elasticity of the surrounding medium. In contrast, in our case, we select real soft materials with smaller differences in shear modulus between media, resulting in a more moderate change in the resonance frequency.

Nonetheless, the results validate our model and confirm that assuming the medium behaves purely as a Newtonian fluid, while neglecting its elastic properties, can lead to substantial inaccuracies in the predicted resonance frequency.

Secondly, we consider the resonance frequency, which differs from the natural frequency due to damping effects, particularly viscous dissipation, given by [28,31]:

$${\omega }_{0^{*}}^2={\displaystyle \frac{1}{\rho R_0^2}}\left[3\gamma P_{b0}+{\displaystyle \frac{2\sigma }{R_0}}(3\gamma -1)+4G\right]-2{\left({\displaystyle \frac{2\mu }{\rho R_0^2}}\right)}^2.$$

(29)

To validate our model Eq. (25) using this expression, we consider the media listed in Table 2, which includes various aqueous glycerol solutions [36] with increasing viscosity. For this analysis, we assume $G=0$ and a constant density $\rho =1000{\mathrm{kg}/\mathrm{m}}^3$.

Figure 2 shows the variation in the linear resonance frequency $f_0^{*}={\omega }_0^{*}/2\pi$, as a function of the initial bubble radius $R_0$ for different viscoelastic media (Table 2). Theoretical predictions (solid lines) based on Eq. (29) are compared with the results obtained using the proposed model Eq. (25) (star symbols). As shown in this figure, the model results closely follow the solid lines, offering further validation of the proposed model equation.

For a bubble with a fixed initial radius, the resonance frequency decreases as the fluid viscosity increases. This trend is observed for $R_0=1\mu \mathrm{m}$, where the resonance frequency drops from $f_0^{*}=3.42$ MHz in GLY00 to about $f_0^{*}=2.75$ MHz in the most viscous solution considered. These results demonstrate that viscosity plays a non-negligible role in determining the natural frequency, even in moderately viscous fluids, due to its dissipative effect on bubble oscillations. This increased damping shifts the resonance frequency to lower values. Consequently, assuming that the fluid behaves as inviscid can lead to an overestimation of the natural frequency of the bubble, particularly in viscoelastic or biological media where viscous effects cannot be neglected.

3.1.2. Bubble Resonance in Representative Soft Viscoelastic Media

We now verify our results against a set of representative soft media. Figure 3 presents how the linear resonance frequency $f_0^{*}={\omega }_0^{*}/2\pi$ varies with the initial bubble radius $R_0$ for each medium listed in Table 3, characterized by their density $\rho$, shear elasticity G, and viscosity $\mu$. The numerical results of our model (depicted as star symbols) show a clear agreement with the theoretical predictions (solid lines) obtained from Eq. (29), which confirms the precision of our model in capturing the resonance behavior.

The results highlight the critical role of the viscoelastic properties of the surrounding medium in shaping the resonant behavior of the bubble.

3.2. Bubble Volume Dynamics over Time

3.2.1. Bubble Volume Variation in Representative Soft Media

To illustrate the capabilities of the model, we examine the behavior of an air bubble with an initial radius $R_0=4.5\mu \mathrm{m}$ excited by an ultrasonic wave at its natural frequency $f=f_0$. The different media analyzed in this study are those listed in Table 3.

Figure 4 presents the variation of the bubble volume obtained from the proposed model Eq. (25) as a function of time under two excitation conditions: (a) an infinitesimal source amplitude $p_0=1\mathrm{mPa}$ and (b) a high finite source amplitude $p_0=5\mathrm{kPa}$.

For infinitesimal amplitudes, the system exhibits a linear response characterized by a vertically symmetric sinusoidal waveform. A substantial difference is observed in the amplitude of the volume variation across the different media, with water exhibiting the largest oscillations because of the absence of viscoelastic resistance.

However, under high finite-amplitude excitation, the response becomes nonlinear in the purely liquid medium (water). The signal is visibly distorted and asymmetric with respect to the $v=0$ axis, indicating a nonlinear behavior. In the case of viscoelastic media, the response shows a more restrained increase of amplitude, and the overall behavior remains within the linear regime despite the higher excitation. This confirms that the viscoelastic properties of the medium significantly influence the dynamics of bubble oscillations. In particular, it mitigates certain nonlinear components of the bubble oscillation. Moreover, the results show that not only an increase in the shear modulus (rigidity), but also an increase in viscosity leads to a notable reduction in the amplitude of bubble oscillations, in agreement with previous studies such as [38].

The model captures qualitative differences in bubble dynamics by identifying distinct oscillation patterns across media with markedly different rheological properties. Despite the simultaneous variation in elasticity, viscosity, and density, the individual contribution of each property, particularly elasticity, remains unclear, which motivates the following section.

3.2.2. Effect of Medium Elasticity on Bubble Dynamics

A more detailed parametric study is conducted with our proposed model Eq. (25) to isolate the specific effects of elasticity on the dynamics of bubble. We analyze the oscillatory behavior of a bubble with an equilibrium radius $R_0=4.5\mu \mathrm{m}$, subjected to an infinitesimal acoustic amplitude at the source ($p_0=1\mathrm{mPa}$), on the one hand, and to a finite acoustic amplitude at the source ($p_0=5\mathrm{kPa}$), on the other. The shear modulus is varied within the range $G=0\text{–}100\mathrm{kPa}$, covering values representative of soft biological tissues [39]: $G=0$, $G=50$, and $G=100\mathrm{kPa}$. The viscosity and density are set to $\mu =0.0014\mathrm{Pa}·\mathrm{s}$ and $\rho =1000{\mathrm{kg}/\mathrm{m}}^3$, corresponding to water as the host medium.

The temporal evolution of the bubble volume variation obtained from the proposed model Eq. (25) is illustrated in Figure 5 and Figure 7, for infinitesimal and finite amplitudes, respectively, to evaluate how elasticity influences expansion behavior. To better capture both transient and steady-state dynamics, the initial and final segments of the simulation are shown in Figure 6a, Figure 6b (linear regime) and Figure 8a, Figure 8b (nonlinear regime). For comparison, results for water ($G=0$), representing the Newtonian medium case, are included as a reference.

In the linear regime (Figure 5, Figure 6a, Figure 6b), the oscillations exhibit a sinusoidal-like behavior and are progressively damped as elasticity increases, indicating a strong stabilizing influence. Transients decay rapidly, and a quasi-steady periodic state is reached, with amplitudes notably reduced in stiffer media. Additionally, the results reveal subtle yet important differences in the oscillation phase between bubbles oscillating in a viscoelastic medium and those in a Newtonian fluid. A slight phase shift is observed, with the peak oscillations in the viscoelastic case occurring slightly earlier than those in the Newtonian case. This shift becomes more noticeable as the elasticity of the medium increases. This phase advancement arises from the elastic component of the Kelvin–Voigt model, which introduces a restorative force that accelerates the return of the bubble to equilibrium, thus altering the timing of the oscillatory response. Whereas the amplitude behavior reflects energy dissipation and nonlinear behavior, the phase shift reveals the memory effects inherent to viscoelasticity.

In contrast, in the nonlinear regime (Figure 7, Figure 8a, Figure 8b), the bubble exhibits larger amplitudes, and elasticity modulates the nonlinear amplification effect. The initial stage (Figure 8a) shows significant transient growth, with bubbles reaching higher expansion rates in the water media. By the final stage of the simulation, periodic oscillations persist, but their amplitude and waveform vary with G, revealing a complex interplay between elasticity and nonlinear resonance. The differences between Newtonian and viscoelastic responses are more pronounced here, underscoring the importance of rheological effects in finite-amplitude acoustic forcing.

3.2.3. Bubble Behavior Across Shear Elasticity-Viscosity Parameter Space

Having examined the behavior of bubbles in well-characterized media and explored the individual effect of elasticity, we now turn to a more application-oriented analysis. This section investigates how the combined viscoelastic properties of soft media influence bubble oscillations, with the aim of optimizing the response depending on the type of tissue or gel in which the bubbles are embedded to achieve the desired effects. To this end, the shear modulus and viscosity are varied in the ranges $G=0$–$100,\mathrm{kPa}$ (as in Section 3.2.2) and $\mu =1$–$10\mathrm{mPa}·\mathrm{s}$, respectively.

Figure 9 illustrates the dependence of bubble volume oscillations on the viscoelastic properties of the surrounding medium, for a given excitation frequency $f=f_0$, by displaying the maximum variation of bubble volume obtained with the proposed model Eq. (25) for different combinations of shear modulus G and viscosity $\mu$. It must be noted that $f=f_0$ is set for each pair of ($G,\mu$), because it changes from one medium to the other. The horizontal axis on the right represents the shear modulus G, the horizontal axis on the left shows the viscosity $\mu$, and the vertical axis, with color scale, corresponds to the maximum variation in the volume of the bubble. The diagram (Figure 9a) corresponds to an infinitesimal excitation amplitude $p_0=1\mathrm{mPa}$ (linear regime), whereas the diagram (Figure 9b) displays the response under a finite excitation amplitude $p_0=5\mathrm{kPa}$ (nonlinear regime).

In both cases, increasing either the shear modulus or the viscosity results in a reduction in the amplitude of bubble volume oscillations. This trend reflects the enhanced mechanical resistance and energy dissipation provided by stiffer and more viscous media, which suppress bubble expansion. The effect is particularly pronounced along the viscosity axis, where even moderate increases in $\mu$ significantly attenuate the bubble response.

In the linear regime, for low-viscosity media ($\mu \le 2\mathrm{mPa}·\mathrm{s}$), a noticeable drop in the amplitude of the oscillation occurs beyond a threshold shear modulus of approximately $G=60\mathrm{kPa}$, indicating that elasticity strongly suppresses bubble vibration when the stiffness is sufficiently high. Moreover, in the same linear regime, when the viscosity exceeds approximately $\mu =3\mathrm{mPa}·\mathrm{s}$, viscous damping dominates the dynamics, making the influence of the shear modulus negligible. In contrast, in the nonlinear regime, the bubble exhibits significantly larger oscillation amplitudes across a broad parameter space, as expected, reflecting stronger volumetric deformation when the driving pressure is raised. Unlike the linear case, elasticity in low-viscosity media does not suppress bubble oscillations, and the influence of viscosity becomes noticeable only for $\mu >8\mathrm{mPa}·\mathrm{s}$, where stabilization occurs.

Given the significant variation in mechanical properties across biological tissues and synthetic gels, the dynamic response of the bubble is highly sensitive to the viscoelastic characteristics of the surrounding medium. The parametric map in the G–$\mu$ space (shear modulus and viscosity) reveals regions (three in Figure 9a and three in Figure 9b, delimited by dashed lines) where oscillation is maximized or suppressed. In particular, areas can be identified where viscosity dominates over elasticity, or vice versa, providing insight into the relative influence of these two parameters.

This information enables the identification of regions of practical interest. For example, areas with strong bubble response may be targeted for drug delivery applications, where inertial cavitation improves transport and release mechanisms [40]. Note that our bubble model is valid for stable cavitation bubbles that are generated just before the onset of inertial cavitation. However, regions of weak response may be preferable in contexts where tissue damage must be minimized [41].

If we extract specific profiles from this figure to study in detail the differences between the linear and nonlinear regimes, Figure 10 shows how the normalized maximum amplitude of bubble volume variation varies with increasing shear modulus G of the medium, for two significant values of viscosity $\mu$, under both excitation regimes. The volume response is normalized with respect to the maximum value at the lowest G, allowing comparison between regimes.

For a medium with low viscosity, such as $\mu =2\mathrm{mPa}·\mathrm{s}$, the normalized bubble response decreases with increasing shear modulus G in both linear and nonlinear regimes, due to the stiffening effect of elasticity. However, the nonlinear regime exhibits a noticeably stronger reduction in the oscillation amplitude, especially at low G. This indicates that elasticity exerts a stronger influence in mitigating large-amplitude volumetric oscillations, resulting in a clearer divergence between linear and nonlinear responses as G increases. On the other hand, for higher viscosity $\mu =8\mathrm{mPa}·\mathrm{s}$, both regimes exhibit a more gradual and similar decay in the amplitude of the oscillation with increasing G, suggesting that elasticity does not play a leading role in reducing the response of the bubble when the viscosity is raised.

Similarly, for two fixed values of the shear modulus G, we examine the evolution of the normalized maximum amplitude of bubble volume variation as the viscosity $\mu$ increases. Figure 11 shows that in both cases, $G=10\mathrm{kPa}$ and $G=60\mathrm{kPa}$, the bubble response decreases with increasing viscosity, reflecting the dissipative effect of viscosity on bubble oscillations. However, this decrease differs from the behavior observed when varying the shear modulus. Three distinct viscosity regions can be identified. In the low-viscosity range ($\mu \in [1,2]\mathrm{mPa}·\mathrm{s}$), the nonlinear regime exhibits a significantly higher normalized amplitude than the linear one, particularly at low elasticity. This suggests that nonlinear excitation is capable of enhancing the bubble response before viscous effects become dominant. In the intermediate-viscosity range ($\mu \in [2,5]\mathrm{mPa}·\mathrm{s}$), viscous damping becomes more influential, and the difference between linear and nonlinear responses begins to narrow, although nonlinear amplification remains visible. Finally, for high viscosities ($\mu >5\mathrm{mPa}·\mathrm{s}$), both regimes display strong attenuation and converge to similar values, indicating that nonlinear effects are largely suppressed and the response is governed by viscous dissipation.

At low elasticity ($G=10\mathrm{kPa}$), the nonlinear regime maintains a higher normalized amplitude throughout the viscosity range, highlighting its ability to partially counteract damping. However, the increasing viscosity progressively reduces this difference. In contrast, at higher elasticity ($G=60\mathrm{kPa}$), both regimes exhibit stronger attenuation and a more similar trend, with minimal nonlinear enhancement. This suggests that in media with high stiffness and viscosity, the oscillatory response of the bubble is significantly minimized, and nonlinear effects lose their ability to modify its dynamics.

Therefore, selecting the appropriate values of shear modulus G and viscosity $\mu$, for fixed excitation amplitude and frequency, turns out to be a critical design criterion. These parametric diagrams thus serve as a practical tool for tayloring the properties of the host medium to achieve the desired cavitation responses in soft tissues, hydrogels, or engineered viscoelastic materials.

4. Conclusions

We have proposed a nonlinear ordinary differential equation to model the nonlinear oscillations of a gas bubble subjected to an ultrasonic field, formulated in terms of bubble volume variation within a viscoelastic medium and valid in the stable cavitation framework. The model has been developed considering a second-order approximation of the adiabatic gas law in the Rayleigh–Plesset equation that is valid for an arbitrary viscoelastic medium. The Kelvin–Voigt model has been incorporated to account for the viscoelastic response of the surrounding medium, thus coupling the classical nonlinear Rayleigh–Plesset dynamics with a linear viscoelastic constitutive relation. The proposed model has been solved numerically and the results have been compared with theoretical resonance curves, showing close agreement and validating our formulation. As expected, the resonance frequency increases with the shear modulus due to enhanced stiffness, while increasing viscosity introduces damping, which slightly reduces the resonance frequency. This approach offers a consistent framework for analyzing how the viscoelastic properties of the medium, specifically shear elasticity, influence the volumetric dynamics of bubbles. This has revealed significant deviations from Newtonian predictions: in particular, the oscillation phase shifts slightly due to elastic energy storage, i.e., the maximum expansion of the bubble occurs earlier in viscoelastic media. These effects become more pronounced with increasing shear modulus. By analyzing bubble dynamics over the shear elasticity–viscosity parameter space, it is evident that both viscosity and elasticity play crucial roles in stable cavitation behavior. Elasticity has a pronounced effect in low-viscosity media, whereas viscosity emerges as the dominant factor modulating the amplitude of oscillations in both the linear and nonlinear regimes. These results emphasize that accurately modeling bubble volume variation requires accounting for the full viscoelastic behavior of the medium, which is essential for predicting and controlling stable bubble dynamics. This work thus provides a robust tool for analyzing how viscoelastic properties affect bubble dynamics, contributing to improved prediction and control of stable cavitation phenomena in complex media. Future work should explore the use of the fourth-order approximation in the Rayleigh–Plesset model [25] to approach the bubble dynamics in this context, which is expected to be useful for driven frequencies close to the bubble resonance. Also, the bubble equation developed here should be taken into account to model the mutual nonlinear influence of ultrasound and populations of multiple bubbles in viscoelastic media.