A unifying Rayleigh-Plesset-type equation for bubbles in viscoelastic media Open Access

The growth and collapse of bubbles in viscoelastic materials is a physical process that features in a variety of environmental, food processing, and industrial settings. From the bubbles formed during volcanic eruptions (Ichihara, 2008) to those used in the production of polymer foams (Andrieux , 2018; Everitt , 2003), understanding the influence of viscoelasticity on the bubble dynamics is a vital step for technological innovation. In particular, the crucial importance of viscoelasticity has been a topic of focus for biomedical applications (Dollet , 2019; Gaudron , 2015; Yang and Church, 2005). Indeed, micrometer-size bubbles are frequently used in ultrasound as contrast agents to visualize organ perfusion and blood flow down to the smallest capillary vessels (Frinking , 2020). They are also heavily investigated as therapeutic agents for cancer treatment, blood-brain-barrier opening (Alonso, 2015; Beccaria , 2013; Choi , 2011; Sheikov , 2008; Sulheim , 2019), sonoporation (Deprez , 2021; Lentacker , 2014; Sirsi and Borden, 2012), lithotripsy, histotripsy (Xu , 2021), or even sonothrombolysis (Bader , 2016; de Saint Victor , 2014; Miller , 2018). These therapeutic applications often require inducing controlled microbubble oscillations within tissue to impart treatment while preventing collateral damage to the surrounding healthy tissue. However, this embedding in a polymeric liquid or a soft solid has a dramatic effect on bubble dynamics (Dollet , 2019). More specifically these materials have viscoelastic properties and provide solid-like or liquid-like resistance to bubble motion, depending on the deformation time scale. Therefore, the ability to control and utilize bubbles in biomedical settings crucially relies on our understanding of the micro and macro-rheology of the surrounding medium and its impact on bubble oscillations.

Various biological materials exhibit viscoelastic properties that come with intricate behavior of stress and strain. A simple method to understand such responses consists of modelling the viscoelastic material as a combination of Hookean springs and Newtonian dashpots (Bird , 1987; Bland, 2016; Kelly, 2013). The spring introduces elastic characteristics to the material's response, while the dashpot adds a viscous contribution. Models built from a combination of springs and dashpots typically lead to exponential relaxation of stress or strain. Even though the spring-dashpot approach accurately describes various viscoelastic liquids and solids, it fails to capture the rheology of multi-scale materials such as gels or elastomeric networks (Karpitschka , 2015; Long , 1996; Winter and Chambon, 1986), which often come along with power-law stress relaxation rather than exponentials. An alternative and more general approach consists of characterizing the material via its complex modulus $μ ( ω ) = G ′ ( ω ) + i G ″ ( ω )$⁠, which combines an elastic storage modulus $G ′ ( ω )$ and viscous loss modulus $G ″ ( ω )$⁠. These moduli can be experimentally obtained by measuring the deformation of the material in response to a periodic excitation, and are functions of the imposed angular frequency ω. As an example, the complex modulus of gels typically exhibits a non-integer scaling behavior with respect to frequency, which cannot be represented by (a finite number of) springs and dashpots. An equivalent formulation of the rheology in the temporal domain can be achieved via the stress relaxation function $ψ ( t )$⁠. This function describes how the stress relaxes after a material is subjected to a step-strain, and can be related to the complex modulus as $μ ( ω ) = i ω ∫ 0 ∞ d t ψ ( t ) e − i ω t$ (Carcione , 2019); hence, $μ ( ω )$ and $ψ ( t )$ contain the exact same rheological information. Within this framework, $ψ → 0$ at large times for viscoelastic liquids, reflecting how the stress fully relaxes after stopping the flow. Viscoelastic solids, by contrast, exhibit a finite ψ at large times, and thus a residual elastic stress, expressing a long-term memory of the reference configuration at rest.

The dynamical behavior of spherical bubbles is typically described using the Rayleigh-Plesset equation, initially developed for spherical bubbles oscillating in a free-field. Possible ways of adding the effects of the viscoelastic surrounding on the bubble dynamics in the context of the Rayleigh-Plesset equation have been the focus of many studies, starting with the pioneering work of Fogler and Goddard (Fogler and Goddard, 1970, 1971). Assuming a linear Maxwell model, the authors investigated the effects of elasticity in the cavitation of a spherical void. Subsequent studies extended this approach by considering different constitutive models of the surrounding material (Allen and Roy, 2000a; Cunha and Albernaz, 2013; Hamaguchi and Ando, 2015; Shima , 1986; Tanasawa and Yang, 1970; Yang and Church, 2005), including those applicable to large-amplitude bubble deformations (Allen and Roy, 2000b; Estrada , 2018; Gaudron , 2015; Jiménez-Fernández and Crespo, 2005; Kim, 1994; Naude and Mendez, 2008; Papanastasiou , 1984; Ting, 1975; Warnez and Johnsen, 2015). Yet, a unifying equation valid for large deformations and applicable to materials with arbitrary complex rheology [equivalently, arbitrary $ψ ( t )$], is still lacking.

Two remarks are in order here. First, in the limit of small deformations where $t ′ → t$⁠, the time derivative of the Finger tensor reduces to $∂ B ( t , t ′ ) / ∂ t ′ = − ϵ ˙ ( t )$⁠, where $ϵ ˙ ( t ) = ∇ v + ( ∇ v ) T$ is the rate of strain tensor ( Appendix B, see also Essink, 2022). Inserting this expression into Eq. (4), we obtain the conventional memory integral for the stress at small deformations (Bird , 1987), as is frequently used in the context of bubbles (Dollet , 2019; Fogler and Goddard, 1970). In general, however, $∂ B ( t , t ′ ) / ∂ t ′$ is not equal to $ϵ ˙$⁠, and this distinction is essential at large deformations. Second, Eq. (4) reflects the deformation history of a certain material point. By consequence, the integral must be carried out at a constant material point and calls for a Lagrangian description of the problem.

To define the Finger tensor $B ( t , t ′ )$ we introduce a Lagrangian description of the deformation, which features prominently in finite-strain theory. Specifically, this description involves relating the Eulerian position $x = ξ ( X , t )$ in the current configuration at time t to a Lagrangian material point X. In solids, one naturally defines X as the coordinates in the reference configuration, but, in general, one can define $X$ from the configuration at some arbitrary time t₀. The mapping between the two states ξ can be used to evaluate the deformation gradient tensor $F ( t , t 0 ) = ∂ x / ∂ X$⁠. This tensor describes the material stretching through the transformation of a line element $d X$ (i.e., the distance between two material particles) in the reference configuration at time t₀ to the same material line element $d x$ at time t. Using the deformation gradient, the Finger tensor is defined as $B ( t , t 0 ) = F ( t , t 0 ) · F T ( t , t 0 )$⁠.

We contextualize the derived Rayleigh-Plesset equation by considering a set of special cases for $ψ ( t )$⁠, enabling a benchmark with existing literature. The relaxation functions of various models encountered in the literature are summarized in Table I. Some of these can be represented schematically by a spring-dashpot, in which case it is possible to articulate the constitutive relation as a differential equation (see Table I). Note that, gel-like materials that exhibit a power-law relaxation function cannot be represented by a differential constitutive law, and thus must be treated with an integral formulation.

We first consider a Newtonian fluid of viscosity η, for which the relaxation function takes the form $ψ ( t ) = η δ ( t )$⁠, where $δ ( t )$ is the Dirac delta function. Using the convolution of the delta function, the integral in the Rayleigh-Plesset equation reduces to $− 4 η R ̇ / R$⁠, which is the standard viscous contribution. Second, we turn to the Neo-Hookean solid, which is a purely elastic medium with shear modulus G, obtained when the relaxation function is constant $ψ ( t ) = G$⁠. In this case, we can carry out explicitly the memory integral in Eq. (11). Applying the initial condition $R ( t ′ → − ∞ ) = R 0$⁠, with R₀ the radius in the reference configuration, yields an elastic contribution to the Rayleigh-Plesset equation of the form $( G / 2 ) [ 5 − 4 ( R 0 / R ) − ( R 0 / R ) 4 ]$⁠. This term exactly corresponds to the one previously obtained for a bubble inside a Neo-Hookean solid (Gaudron , 2015; Warnez and Johnsen, 2015). Third, we consider a bubble inside an Oldroyd-B fluid, for which the relaxation function takes the form $ψ ( t ) = η δ ( t ) + G exp ( − t / λ )$⁠, where η is the solvent viscosity, λ the relaxation time, and $G = η p / λ$ the polymer's elastic modulus (see also  Appendix C). Here, as well, the resulting Rayleigh-Plesset equation is identical to that reported by Ting, who considered the Oldroyd-B model in the context of bubble cavitation (Ting, 1975).

To further validate our model and its numerical implementation, we compare it to the results obtained by Allen and Roy (2000b). Specifically, Allen and Roy considered large deformations of an acoustically driven micron-sized bubble. The oscillatory pressure field is $p ∞ = p 0 + p a sin ( 2 π f t )$⁠, where p₀ is the ambient pressure, p_a is the pressure amplitude, and f the driving frequency. The gas pressure takes the form $p g = ( p 0 + 2 γ / R 0 ) ( R 0 / R ) 3 k$⁠, to account for thermal dissipation through the polytropic constant k. The bubble is surrounded by virtual Upper Convected Maxwell fluids (UCM) with relaxation times $λ =$ 0, 0.5, and 1 μs. The ratio between the relaxation time and the characteristic time scale of the flow is expressed by the Deborah number, which takes the values $De = 2 π f λ =$ 0, 1, and 2 for the three relaxation times considered.

Our model is in excellent agreement with the results reported by Allen and Roy, as shown in Fig. 2. Importantly, the numerical solution in Allen and Roy (2000b) is not based on the integral form of the constitutive equation, but on the differential form. Consequently, the stress field outside the bubble must be solved numerically, and the spatial integral of the stress needed in the Rayleigh-Plesset equation must also be evaluated numerically. Hence, the agreement in Fig. 2 offers a nontrivial validation of the proposed modelling framework. We also wish to highlight that Eq. (1) has the advantage of providing an autonomous equation for the bubble radius R(t), which no longer relies on a separate (numerical) evaluation of the stress outside the bubble.

Having established the validity of our model and its numerical implementation, we now further explore the oscillatory motion of bubbles in viscoelastic materials. We study the motion of a 2 μm sized bubble driven by a Gaussian-tapered pressure pulse with 16 acoustic cycles [Fig. 3(a)], pressure amplitude of 50 kPa, and a frequency in the range $0.5 ≤ f ≤ 8$ MHz. For material properties, we select a liquid density $ρ = 1 , 000$ $kg / m 3$⁠, and for most cases, we fix the surface tension to a constant value γ = 72 mN/m. For viscoelastic liquids we consider the bubble to be surrounded by an Oldroyd-B liquid or a Critical Gel, whose relaxation functions are $ψ ( t ) = G exp ( − t / λ )$ and $ψ ( t ) = G ( λ / t ) 1 / 2$⁠, respectively (see Table I). For the Oldroyd-B liquid we set the solvent viscosity to η = 2  $mPa · s$ and we test the viscoelastic effects by varying the values of the shear modulus G and relaxation time λ. Specifically, we consider three shear moduli $G = 10 , 100 , and 500$ kPa and two relaxation times of $λ = 0.01 and 1.00$ μs. The same parameters are used for the Critical Gel to enable a direct comparison. As an additional perspective, we will also consider an example where the bubble is coated by a layer of lipids, modelled via additional mechanical properties at the interface (see the following section). By comparing the resonance behavior of coated and uncoated bubbles, one further appreciates the relative importance of the viscoelasticity of the surrounding medium.

The results of our simulations for the evolution of the bubble radius R(t) in a Newtonian and Oldroyd-B fluid are shown in Fig. 3. Considering first a bubble in a Newtonian liquid, the oscillation amplitude of the bubble increases initially as it follows the applied pressure pulse. The bubble radius reaches stable oscillations with a peak-to-peak amplitude $Δ R = 0.4 R 0$⁠, before decreasing to its initial radius, as the pressure pulse decays [Fig. 3(b)]. Switching to an Oldroyd-B fluid with relaxation time λ = 1 μs and shear modulus G = 500 kPa, the amplitude decreases to approximately 20% of R₀. We thus observe that the resistance of the viscoelastic stresses can significantly affect the bubble oscillation. To further examine the effects of the surrounding medium, it is more instructive to evaluate the resonance behavior of the bubble. Indeed, bubbles are known to act as damped harmonic oscillators when subjected to oscillatory pressure driving (Minnaert, 1933; Prosperetti, 1977).

To this end, we obtain the response amplitude of the bubble through a fast Fourier transform and plot it as a function of the driving frequency. The fast Fourier transform of the oscillation amplitude is also indicative of the presence of higher or lower harmonics [Fig. 3(c)]. The resulting normalized resonance curves of $Δ R / R 0$ are shown in Fig. 4. For a 2 μm radius bubble in a Newtonian liquid and a polytropic constant k = 1.4, the resonance frequency is very close to the Minnaert frequency $f 0 ≈ 1 / ( 2 π R 0 ) ( 3 k p 0 / ρ ) 1 / 2 ≈ 1.6$ MHz (Minnaert, 1933), where the bubble displays a maximum relative amplitude of oscillation of 0.42 [Fig. 4(a)]. The viscoelastic effects become immediately apparent for the Oldroyd B liquid; with a relaxation time $λ = 0.01$ μs, the viscoelastic effects almost exclusively translate into a reduction of the oscillation amplitude. Note that G = 10 kPa does not give rise to a significant difference as compared to the Newtonian case. By contrast, shear moduli of 100 and 500 kPa decrease the response amplitude by 20% and 50%, respectively. In all three cases, the resonance frequency remains unaffected. The effects of viscoelasticity change qualitatively when increasing the relaxation time to λ = 1 μs [Fig. 4(b)]. Even though G = 10 kPa does not significantly change the bubble behavior, setting G = 100 and 500 kPa results in a drastic increase in resonance frequency (20% and 100%, respectively). In addition, there is a further decrease in amplitude as compared to the case of $λ = 0.01$ μs.

As an additional perspective, we consider the resonance behavior of coated microbubbles. Microbubbles can increase the contrast in ultrasound imaging or, once embedded in tissue, be used as controlled cavitation agents for therapy. By contrast to the free gas bubble considered until now, these bubbles must be coated with a shell to prevent gas diffusion driven by the excess Laplace pressure. The coating has two additional contributions to the bubble dynamics (see  Appendix D). First, the dynamic change in surface area induces molecular friction, which leads to a significantly increased damping contribution (Van der Meer , 2007). Second, the surface tension becomes a function of the packing density of the phospholipid molecules and, therefore, for a certain initial packing fraction, of the bubble size (Marmottant , 2005). When the bubble size decreases, the lipids get compressed leading to shell buckling like a spherical elastic shell. Using the viscoelastic Rayleigh-Plesset equation for coated microbubbles (see  Appendix D), we test the resonance curve of a coated microbubble inside an Oldroyd-B fluid with λ = 1 μs and G = 500 kPa. We take a surface elasticity χ = 1 N/m and a shell dilatational viscosity of κ_s =  $10 − 8$ kg/s. The resulting resonance behavior of the coated microbubble is shown in Fig. 5(b). Compared to the result for the uncoated bubbles [Fig. 5(b)], two features stand out. First, the resonances shift to slightly higher frequencies, which can be attributed to the additional stiffness induced by the surface elasticity of the coating. Second, the peaks are wider, and of smaller amplitude, due to the additional damping induced by the shell's viscosity. Yet, the overall features of the resonance curves, specifically the trends with λ, remain unaffected by the coating, showing that these are dictated by the viscoelastic properties of the surrounding medium.

Polymer networks typically exhibit a broad spectrum of relaxation times, and therefore cannot be described by a simple exponential decay. For example, silicone gels are well-described by the so-called Chasset-Thirion model (cf. Table I), which exhibit an algebraic decay in the range $0 < n < 1$ (Karpitschka , 2015; Long , 1996). Here, we focus on the special case of a material at the gelation point, for which the equilibrium shear modulus vanishes. Following Winter and Chambon (1986), the Critical Gel corresponds to the case where $G ′ = G ″$ at all frequencies [Fig. 4(d) inset]. This criterion is achieved only for the special case of $ψ = G ( λ / t ) 1 / 2$⁠, which gives $G ′ = G ″ ∼ G λ 1 / 2 ω 1 / 2$⁠. Such a material at the gelation point has the peculiar property that, at low frequency, it has a vanishing shear modulus (⁠ $lim ω → 0 G ′ = 0$⁠) and an infinite steady viscosity (⁠ $lim ω → 0 G ″ / ω = ∞$⁠). Conversely, at large frequency, the storage modulus diverges (⁠ $lim ω → ∞ G ′ = ∞$⁠) while the effective viscosity vanishes (⁠ $lim ω → ∞ G ″ / ω = 0$⁠).

The resonance behavior of free gas microbubbles inside a Critical Gel is indeed very different from that in an Oldroyd-B fluid. We find a stronger reduction in the amplitude of oscillation. In addition, the Critical Gel also induces a much stronger resonance frequency shift. This shift is not even visible for the lowest relaxation time of $λ = 0.01$ μs [Fig. 4(c)], which did not induce a significant shift for the Oldroyd-B fluid. This is a direct consequence of the absence of a characteristic time scale in the Critical Gel. In fact, even though we separately defined G and λ, the gel is entirely characterized by the single combination $G λ 1 / 2$⁠, which is also known as the strength of the gel (Winter and Chambon, 1986). Therefore, the results for a larger relaxation time of λ = 1 μs [Fig. 4(d)], shows similar trends as those for $λ = 0.01$ μs. More quantitatively, we can see in Figs. 4(c) and 4(d) that the combination of [G = 100 kPa, $λ = 0.01$ μs] [green curve in Fig. 4(c)] and [G = 10 kPa, λ = 1 μs] [red curve in Fig. 4(c)] result in identical curves, since the product of $G λ 1 / 2$ is the equal. Thus, for the gel, in contrast to the Oldroyd-B fluid, increasing λ has a similar effect as increasing G.

The previous results illustrate the importance of the rheology for the resonance behavior of microbubbles; in particular when comparing exponentially decaying relaxation functions to gel-like materials that exhibit a power-law spectrum. The fact that the strength of a critical gel involves the combination $G λ 1 / 2$ prevents the existence of an “elastic limit” like for the Oldroyd-B fluid since the equivalent solid would naturally have an infinite stiffness. The absence of an elastic limit at $λ → ∞$ is not limited to the Critical Gel, but also applies to the Chasset-Thirion model that is used to describe silicone gels.

Microbubbles driven at sufficiently high pressure amplitudes can undergo highly non-linear oscillations, generating higher order harmonic responses of the bubble (Lauterborn, 1976; Prosperetti, 1974), as well as subharmonic behavior (Neppiras, 1969). In the context of viscoelastic fluids, Allen and Roy showed that the presence of elasticity can enhance the second harmonic response of an Upper Convected Maxwell fluid (Allen and Roy, 2000b). In this section, we briefly comment on the second harmonic by examining the fast Fourier transform of the bubble response [Fig. 3(c)]. We record the oscillation amplitude $Δ R 2 H$ at twice the driving frequency and directly compare the results with respect to the fundamental frequency in the context of the Oldroyd-B and Critical Gel of Fig. 4. The resonance curve of the second harmonic is shown in Fig. 6. For the Oldroyd-B fluid with a relaxation time of $λ = 0.01$ μs we observe the appearance of two peaks for each curve [Fig. 6(a)]. The main peak in the response occurs at approximately half the fundamental frequency, and a smaller one at the fundamental frequency. The higher peak at half the fundamental frequency reflects the preference of the bubble to oscillate at its eigenfrequency. An increase in the shear modulus dampens the curves and leads to smaller amplitudes. The relaxation time is too short as compared to the oscillation period, and thus excites a viscous response from the Oldroyd-B fluid. Upon increasing the relaxation time to λ = 1 μs, the resonance frequency is shifted for the higher shear moduli [Fig. 6(b)]. The two peaks still appear at the fundamental and half the fundamental frequency but are shifted depending on the elasticity. We thus find that the effects of the Oldroyd-B fluid on the second harmonic are similar to those observed for the fundamental frequency. The enhancement of the second harmonic response reported by Allen and Roy could therefore be attributed to the non-monotonic behavior of the Oldroyd-B model in the elastic limit. Finally, the second harmonic of the Critical Gel also shares the same characteristics, as the trends follow the behavior of the fundamental frequency [Fig. 6(c)]. Two peaks also appear for certain values of $G λ 1 / 2$⁠, with the amplitude decreasing as the shear modulus is increased.

In this article, we have derived a unifying Rayleigh-Plesset-type equation for bubbles in viscoelastic materials. Using concepts from finite linear viscoelasticity, we express the viscoelastic stresses as a function of the deformation history. This approach has the benefit of extending the Rayleigh-Plesset equation to a broad class of viscoelastic materials with arbitrary complex moduli through their relaxation function $ψ ( t )$⁠, while consistently accounting for large deformations. Not only were we able to capture the bubble dynamics in materials with well-defined constitutive relations, such as the Upper Convected Maxwell fluid, but also in gel-like materials that can only be characterized by their storage and loss moduli, and not via a constitutive differential equation. With these ideas in mind, we tested the resonance behavior of a microbubble in different types of viscoelastic materials when driven by a tapered sinusoidal pressure pulse. The broad spectrum of relaxation times of gel-like materials completely alters the resonance behavior of a microbubble for sufficiently stiff shear moduli. In contrast, the Olroyd-B fluid was found to behave both as a viscous liquid and an elastic solid, depending on the relaxation time. On the one hand, the Oldroyd-B fluid acts as a viscous liquid when the relaxation time is much faster than the typical oscillation period. Conversely, for relaxation times much slower than the typical oscillation period, the resonance curve of an Oldroyd-B fluid converges to that of a Kelvin-Voigt solid.

Our model is applicable to various viscoelastic materials; yet, it still has its limitations. Even though the bubble deformation is nonlinear in terms of the bubble radius R, our stress formulation of Eq. (3) can be classified as finite linear viscoelastic, as the memory integral is still based on a superposition principle, involving a linear operator in B. A possible extension to nonlinear response, while preserving the superposition principle, is by taking the relaxation function $ψ ( t )$ to be dependent on the deformation (Wineman, 2009). For the derivation of Eq. (1), however, the assumption that $ψ ( t )$ is independent of deformation is crucial in our approach as it enabled us to separate the spatial and temporal integrals of the stress that lead to Eq. (11). More generally, there are many constitutive equations that are not described by Eq. (3). Typical examples include models that exhibit strain stiffening (Yang , 2022; Yang , 2020) or finite extensibility, which limits the extent to which the material can be stretched, such as the Gent model for solids (Gent, 1996) or the FENE-P model for fluids (Warner, 1972). Furthermore, viscoplastic materials also do not fall into the class of Eq. (3) (De Corato , 2019). The key advantage of Eq. (1) therefore lies in the applicability to materials of arbitrary stress relaxation function, while preserving the possibility of finite deformations.

Another limitation of our model is that we consider strictly radial bubble motions, i.e., purely spherical volumetric oscillations. This assumption can be consequential for sufficiently large oscillations, as the shape of a microbubble is susceptible to parametric instabilities when driven at large acoustic pressures (Versluis , 2010). The assumption of purely radial motions is also crucial when extending our model to coated microbubbles, which are known to exhibit non-spherical oscillations either as a result of a symmetry breaking in the medium or as surface modes develop (Dollet , 2019; Dollet , 2008). As a future perspective, it would thus be worthwhile to focus on the non-spherical motion of coated and uncoated microbubbles in viscoelastic materials. Indeed, our theoretical approach lays a generalized framework for bubble dynamics in viscoelastic media, which, as we have exemplified, could also be extended to model viscoelastic dissipation at the interface of coated microbubbles.

There are no conflicts of interest to declare.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Here, we provide an overview of the kinematic relations that can be used to better interpret the stress Eq. (4) in certain limits. Recall that the deformation gradient tensor is defined as $F ( t , t ′ ) = ∂ x / ∂ x ′$⁠. Taking the time derivative $d / d t$ at a constant material point and using the chain rule, one can show that (Essink, 2022; Morozov and Spagnolie, 2015; Snoeijer , 2020; Stone , 2023)