1. Introduction
Motivated by negative outcomes such as vibrations, noise and erosion (Hammitt Reference Hammitt1980; Brennen Reference Brennen1995), cavitation research originally focused on hydrodynamically produced vapour bubbles in water in the context of naval engineering and turbomachinery, starting with the work of Rayleigh (Reference Rayleigh1917). Plesset (Reference Plesset1949) generalized this approach by incorporating a time-varying pressure, viscosity and surface tension into the celebrated Rayleigh–Plesset equation, which describes the dynamics of a single spherical bubble. Subsequent refinements to the theory include heat transfer (Hickling & Plesset Reference Hickling and Plesset1964; Stricker, Prosperetti & Lohse Reference Stricker, Prosperetti and Lohse2011), compressibility (Herring 1941; Gilmore Reference Gilmore1952; Keller & Miksis Reference Keller and Miksis1980; Prosperetti & Lezzi Reference Prosperetti and Lezzi1986) and other effects summarized in Plesset & Prosperetti (Reference Plesset and Prosperetti1977). More recently, advances in high-speed imaging and pressure transducers (reviewed in Lauterborn & Kurz Reference Lauterborn and Kurz2010) and numerical simulations, incompressible (e.g. Blake & Gibson Reference Blake and Gibson1987) and compressible (e.g. Johnsen & Colonius Reference Johnsen and Colonius2009), have expanded our understanding of cavitation in Newtonian liquids.
Interest in hydrodynamic cavitation in non-Newtonian liquids started when it was demonstrated that even a small amount of macromolecular additive inhibits cavitation inception (Ellis & Hoyt Reference Ellis and Hoyt1958). Theoretical studies of cavitation in viscoelastic fluids based on Maxwell models have further shown that viscoelasticity retards and damps the bubble collapse (Fogler & Goddard Reference Fogler and Goddard1970, Reference Fogler and Goddard1971; Shima & Tsujino Reference Shima and Tsujino1976), introduces memory effects that are dominant at early times (Papanastasiou, Scriven & Macosko Reference Papanastasiou, Scriven and Macosko1984) and reduces the growth of non-spherical perturbations by increasing viscosity or amplifies it by increasing relaxation (Hara & Schowalter Reference Hara and Schowalter1984). Experimental investigations have shown that polymer additives decrease the speed of the re-entrant jet in non-spherical collapse (Brujan et al. Reference Brujan, Ohl, Lauterborn and Philipp1996). Bubbles in shear-thinning non-Newtonian media like blood have been investigated in the context of intravascular gas embolism as well (Mukundakrishnan, Eckmann & Ayyaswamy Reference Mukundakrishnan, Eckmann and Ayyaswamy2009). While some direct simulations of the equations of motion have been performed (Foteinopoulou & Laso Reference Foteinopoulou and Laso2010; Lind & Phillips Reference Lind and Phillips2010), none have accounted for compressibility of the surrounding medium, which is expected to play an important role in inertial collapse.
More recently, acoustic cavitation in soft tissue has received interest, following progress in diagnostic and therapeutic ultrasound (Goldberg, Liu & Forsberg Reference Goldberg, Liu and Forsberg1994; Lingeman Reference Lingeman1997; Roberts et al. Reference Roberts, Hall, Ives, Wolf, Fowlkes and Cain2006). In these non-invasive procedures, cavitation generated by incoming pressure waves may occur directly in tissue (Coussios & Roy Reference Coussios and Roy2008), an inherently viscoelastic medium with complex constitutive relations (Fung Reference Fung1993), and cause damage to the surroundings. In particular, histotripsy (Parsons et al. Reference Parsons, Cain, Abrams and Fowlkes2006; Roberts et al. Reference Roberts, Hall, Ives, Wolf, Fowlkes and Cain2006) is a recently developed tissue ablation procedure that has been shown to be capable of treating, in animal models, benign prostatic hyperplasia, prostate cancer, renal masses and renal stones (Roberts Reference Roberts2014). Histotripsy is thought to be a non-thermal cavitation-based procedure, in contrast to high-intensity focused ultrasound in which boiling plays an important role (ter Haar, Sinnett & Rivens Reference ter Haar, Sinnett and Rivens1989); however, isolation of mechanical effects from thermal ones is challenging.
In the context of mechanical tissue fractionation, experiments on cavitation in hydrogels indicate that elasticity inhibits cavitation activity (Maxwell et al. Reference Maxwell, Cain, Hall, Fowlkes and Xu2013; Vlaisavljevich et al. Reference Vlaisavljevich, Maxwell, Warnez, Johnsen, Cain and Xu2014). Due to challenges with experiments, e.g. space/time resolution and optical access, numerical modelling shows promise towards understanding the basic mechanics. Originally, cavitation in tissue was modelled as occurring in a Maxwell-based viscoelastic fluid, e.g. with linear Maxwell (Allen & Roy Reference Allen and Roy2000a ) and nonlinear Oldroyd (Allen & Roy Reference Allen and Roy2000b ; Jimenez-Fernandez & Crespo Reference Jimenez-Fernandez and Crespo2006; Naude & Mendez Reference Naude and Mendez2008) models. Given that soft tissue typically relaxes to an original configuration after deformation, Kelvin–Voigt-based models have been explored more recently, e.g. linear Kelvin–Voigt (Yang & Church Reference Yang and Church2005) and Zener (Hua & Johnsen Reference Hua and Johnsen2013) models. Other studies have focused on the dynamics of bubbles encapsulated in a thin viscoelastic shell: Khismatullin & Nadim (Reference Khismatullin and Nadim2002) considered a linear Kelvin–Voigt shell in an Oldroyd fluid, while Liu et al. (Reference Liu, Sugiyama, Takagi and Matsumoto2012) investigated a neo-Hookean shell in water.
In biomedical applications, large and rapid bubble oscillations may be observed due to the small bubble sizes and high driving-pressure amplitudes, thus leading to large strains and high strain rates. As a result, two important effects absent in most previous studies may be important: compressibility of the surroundings and nonlinear (visco)elasticity. The former is addressed in Yang & Church (Reference Yang and Church2005) and Hua & Johnsen (Reference Hua and Johnsen2013), but the latter remains mostly unexplored. While Oldroyd-based models include nonlinearity of the viscoelasticity in the stress relaxation, they do not account for the presence of an original configuration and the finite-strain elasticity associated with large deformations (for the present purpose, hyperelasticity, finite-strain elasticity and nonlinear elasticity are synonymous). For inertial cavitation, characterized by large strains and high strain rates, linear elasticity is expected to fail.
The present work investigates the dynamics of a single spherical bubble in a homogeneous viscoelastic medium with nonlinear elasticity, ranging from small-amplitude oscillations to violent inertial collapse. An accurate description of elasticity is crucial because of its strong impact on bubble collapse (Fogler & Goddard Reference Fogler and Goddard1970), which is expected to be responsible for damage caused to the surroundings. To derive the equations of motion, we follow modern continuum mechanics finite-strain formalism (e.g. Mooney Reference Mooney1940; Ogden Reference Ogden1972; Holzapfel Reference Holzapfel2000). To the best of our knowledge, only one other study has accounted for nonlinear elasticity in the context of bubble dynamics (Liu et al. Reference Liu, Sugiyama, Takagi and Matsumoto2012), albeit to describe the infinitesimally thin shell encapsulating a bubble in water. In contrast, our work focuses on the more general problem of a bubble located directly in a viscoelastic medium whose elasticity is nonlinear. We consider several common models: linear elasticity, (one-parameter) neo-Hookean strain-energy function and (two-parameter) Mooney–Rivlin strain-energy function. Our approach can readily be generalized to other models. The focus of this article lies in purely mechanical effects; thermal effects are beyond the present scope but could be included. The present article is organized as follows. The governing equations for bubble dynamics in a neo-Hookean medium are first derived. Then, we analytically assess the role of (nonlinear) elasticity under various limits for different models. Next, we provide numerical results to illustrate the effect of (nonlinear) elasticity. Finally, the article summarizes the results and provides an outlook for future work.
2. Governing equations
The dynamics of a single spherical bubble in an infinite homogeneous isotropic medium is described following modern continuum mechanics formalism (e.g. Holzapfel Reference Holzapfel2000). Figure 1 shows the problem set-up. The bubble, whose centre lies at the origin, starts from an initial strain-free configuration with radius
$R_{o}$
. In the ‘current’ configuration at time
$t$
, its radius is denoted
$R(t)$
. Similarly, any other material point initially at
$\boldsymbol{x}_{o}=r_{o}\boldsymbol{e}_{r}$
is denoted
$\boldsymbol{x}=r(r_{o},t)\boldsymbol{e}_{r}$
in the current configuration. The deformation tensor is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D641}=\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{x}_{o}}=\left[\begin{array}{@{}ccc@{}}\displaystyle \frac{\partial r}{\partial r_{o}} & 0 & 0\\ 0 & \displaystyle \frac{r}{r_{o}} & 0\\ 0 & 0 & \displaystyle \frac{r}{r_{o}}\end{array}\right]\!,\end{eqnarray}$$
based on the gradient of a rank-two tensor in spherical coordinates (see, for example, Arfken, Weber & Harris Reference Arfken, Weber and Harris2013). The near field is assumed incompressible, with
$\text{det}(\unicode[STIX]{x1D641})=1$
. Thus,
$$\begin{eqnarray}\frac{\partial r}{\partial r_{o}}=\left(\frac{r_{o}}{r}\right)^{2}\Rightarrow r(r_{o},t)=[r_{o}^{3}+C_{1}(t)]^{1/3},\end{eqnarray}$$
where
$C_{1}(t)$
is an integration factor. Compressibility of the surrounding medium is included later by allowing for far-field wave propagation to provide a mechanism for energy dissipation by acoustic radiation (Herring 1941; Gilmore Reference Gilmore1952).
Figure 1. Problem setting using typical continuum mechanics convention.
Figure 2. Schematic of the present Kelvin–Voigt viscoelastic model consisting of a spring
${\it\eta}$
and a dashpot
${\it\mu}$
in parallel.
The surrounding medium is assumed to behave in a viscoelastic fashion. Given that soft tissue generally relaxes to an original configuration after loading, a Kelvin–Voigt viscoelastic medium with a spring and dashpot in parallel (Figure 2) is used, with the stress tensor given by
$$\begin{eqnarray}\unicode[STIX]{x1D64F}=\unicode[STIX]{x1D64F}_{\boldsymbol{f}}+\unicode[STIX]{x1D64F}_{\boldsymbol{e}}.\end{eqnarray}$$
The parallel spring allows for the medium to eventually recover its original configuration after deformation. In this formulation, a linear or nonlinear spring can readily be considered, corresponding to linear or nonlinear elasticity. Pressure is included in the elastic term, as explained later. The fluid component of the stress is purely deviatoric, given by Newton’s law of viscosity,
$$\begin{eqnarray}\unicode[STIX]{x1D64F}_{\boldsymbol{f}}=2{\it\mu}\unicode[STIX]{x1D64E}={\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\text{T}}),\end{eqnarray}$$
where
${\it\mu}$
is the viscosity,
$\unicode[STIX]{x1D61A}_{ij}$
is the strain-rate tensor and
$u_{i}$
is the velocity. The elastic component is the Cauchy stress tensor, which in its most general form is written
$$\begin{eqnarray}\unicode[STIX]{x1D64F}_{\boldsymbol{e}}=2\left[\left(\frac{\partial {\hat{W}}}{\partial I_{1}}+I_{1}\frac{\partial {\hat{W}}}{\partial I_{2}}\right)\unicode[STIX]{x1D63D}-\frac{\partial {\hat{W}}}{\partial I_{2}}\unicode[STIX]{x1D63D}^{2}\right]-\mathscr{P}{\bf\delta},\end{eqnarray}$$
where
$\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}^{\text{T}}$
is the left Cauchy–Green tensor,
${\hat{W}}$
the strain-energy function describing the surrounding medium,
${\bf\delta}$
is the identity tensor,
$\mathscr{P}(r,t)$
is a pseudo-pressure to be discussed later and
$I_{i}$
are the strain invariants, with
$i=1,2,3$
. For a prescribed strain-energy function, the dynamics is given by the momentum balance equation,
$$\begin{eqnarray}{\it\rho}\boldsymbol{a}=\text{div}(\unicode[STIX]{x1D64F}),\end{eqnarray}$$
where
${\it\rho}$
is the constant density of the surroundings and
$\boldsymbol{a}$
is the acceleration of any material point. By symmetry, only the radial component is retained in this equation, with an Eulerian acceleration in the current configuration given by
$$\begin{eqnarray}a_{r}=\left[\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}\right]\!,\end{eqnarray}$$
where the divergence-free velocity field is
$u(r,t)={\dot{R}}R^{2}/r^{2}$
.
Without loss of generality, the bubble dynamics equations are derived for a neo-Hookean model, the simplest one-parameter nonlinear elastic model. The corresponding strain-energy function is
${\hat{W}}={\it\eta}(I_{1}-3)/2$
(Treloar Reference Treloar1943a
,Reference Treloar
b
; Ogden Reference Ogden1972), where
${\it\eta}$
is the shear modulus from linear theory and
$I_{1}=\text{tr}(\unicode[STIX]{x1D63D})$
is the first strain invariant. Thus, the Cauchy stress tensor can be written
$$\begin{eqnarray}\unicode[STIX]{x1D64F}_{\boldsymbol{e}}={\it\eta}\unicode[STIX]{x1D63D}-\mathscr{P}\unicode[STIX]{x1D644}=\left[\begin{array}{@{}ccc@{}}\displaystyle {\it\eta}\left(\frac{r_{o}}{r}\right)^{4}-\mathscr{P} & 0 & 0\\ 0 & \displaystyle {\it\eta}\left(\frac{r}{r_{o}}\right)^{2}-\mathscr{P} & 0\\ 0 & 0 & \displaystyle {\it\eta}\left(\frac{r}{r_{o}}\right)^{2}-\mathscr{P}\end{array}\right]\!.\end{eqnarray}$$
Since
$\unicode[STIX]{x1D63D}$
is not purely deviatoric,
$\mathscr{P}(r,t)$
is not the hydrodynamic pressure,
$p(r,t)$
. The latter is obtained by writing the left Cauchy–Green tensor in terms of a deviatoric term,
$\text{dev}(\unicode[STIX]{x1D63D})$
, and an isotropic component (pressure),
$$\begin{eqnarray}\unicode[STIX]{x1D64F}_{\boldsymbol{e}}={\it\eta}\text{dev}(\unicode[STIX]{x1D63D})-\left[\mathscr{P}-\frac{{\it\eta}}{3}\text{tr}(\unicode[STIX]{x1D63D})\right]{\bf\delta}={\it\eta}\text{dev}(\unicode[STIX]{x1D63D})-p{\bf\delta}.\end{eqnarray}$$
For the neo-Hookean model,
$p=\mathscr{P}-{\it\eta}[(r_{o}/r)^{4}+2(r/r_{o})^{2}]/3$
. While not absolutely necessary, definition of
$\mathscr{P}$
leads to less complicated derivations and provides a clear connection to traditional continuum mechanics notation; the final results are written in terms of
$p$
. In light of these definitions, the Cauchy stress tensor can be rewritten:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64F}_{\boldsymbol{e}}=\left[\begin{array}{@{}ccc@{}}\displaystyle \frac{2{\it\eta}}{3}\left[\left(\frac{r_{o}}{r}\right)^{4}-\left(\frac{r}{r_{o}}\right)^{2}\right]-p & 0 & 0\\ 0 & \displaystyle \frac{{\it\eta}}{3}\left[\left(\frac{r}{r_{o}}\right)^{2}-\left(\frac{r_{o}}{r}\right)^{4}\right]-p & 0\\ 0 & 0 & \displaystyle \frac{{\it\eta}}{3}\left[\left(\frac{r}{r_{o}}\right)^{2}-\left(\frac{r_{o}}{r}\right)^{4}\right]-p\end{array}\right]\!. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Acknowledging that a point at infinity remains unperturbed, i.e.
$r_{o}/r\rightarrow 1$
at infinity, the pressure at infinity, which drives the bubble dynamics, is
$p_{\infty }(t)=\mathscr{P}_{\infty }(t)-{\it\eta}$
.
Now that
$\unicode[STIX]{x1D64F}$
is fully defined, (2.6) can be integrated from
$R$
to
$\infty$
to yield
$$\begin{eqnarray}R\ddot{R}+\frac{3}{2}{\dot{R}}^{2}=\mathscr{E}(t)+\frac{\mathscr{P}(R,t)-\mathscr{P}_{\infty }}{{\it\rho}}-\frac{4{\it\nu}_{L}{\dot{R}}}{R},\end{eqnarray}$$
where
${\it\nu}_{L}$
is the liquid viscosity and the elastic contribution is
$$\begin{eqnarray}\mathscr{E}(t)\equiv \frac{2{\it\eta}}{{\it\rho}}\int _{R(t)}^{\infty }\left[-\frac{r_{o}^{4}}{r^{5}}-\frac{r}{r_{o}^{2}}+\frac{2r_{o}}{r^{2}}\right]\!\text{d}r,\end{eqnarray}$$
where
$r_{o}=r_{o}(r,t)$
. For
${\it\eta}=0$
, (2.11) reduces to the Rayleigh–Plesset equation.
The bubble is constituted of vapour and contaminant gas, and its pressure is assumed to be homogeneous (Brennen Reference Brennen1995). Thus, the bubble pressure is
$$\begin{eqnarray}p_{B}(t)=p_{v}(T_{\infty })+p_{Go}\left(\frac{R_{o}}{R}\right)^{3k},\end{eqnarray}$$
where
$p_{v}$
is the vapour pressure,
$T_{\infty }$
is the ambient temperature,
$p_{Go}$
is the initial partial pressure of the contaminant gas and
$k$
is the polytropic index. For simplicity, heat and mass transfer are neglected; these effects could readily be incorporated into the model. The dynamic boundary condition at the bubble surface yields
$-p_{B}=\unicode[STIX]{x1D64F}_{rr}|_{r=R}-2S/R$
, where
$S$
is surface tension. Thus,
$$\begin{eqnarray}\mathscr{P}(R,t)=p_{v}+p_{Go}\left(\frac{R_{o}}{R}\right)^{3k}+{\it\eta}\left(\frac{R_{o}}{R}\right)^{4}-\frac{2S}{R}.\end{eqnarray}$$
An expression for
$r_{o}(r,t)$
can be obtained by evaluating (2.2) at the bubble wall:
$$\begin{eqnarray}r_{o}(r,t)=(r^{3}-R^{3}+R_{o}^{3})^{1/3}.\end{eqnarray}$$
Using integration by parts and rearranging, the elastic term yields
$$\begin{eqnarray}\mathscr{E}=-\frac{{\it\eta}}{2{\it\rho}}\left[3+\left(\frac{R_{o}}{R}\right)^{4}-4\left(\frac{R_{o}}{R}\right)\right]\!,\end{eqnarray}$$
thus leading to the equation governing spherical bubble dynamics in an incompressible and viscoelastic medium whose elasticity is described by a neo-Hookean model:
$$\begin{eqnarray}R\ddot{R}+\frac{3}{2}{\dot{R}}^{2}=\frac{p_{B}-p_{\infty }(t)}{{\it\rho}}-\frac{4{\it\nu}_{L}{\dot{R}}}{R}-\frac{2S}{{\it\rho}R}-\frac{E}{{\it\rho}},\end{eqnarray}$$
where the elastic stress for a neo-Hookean medium is given by
$$\begin{eqnarray}E_{NH}=\frac{{\it\eta}}{2}\left[5-4\left(\frac{R_{o}}{R}\right)-\left(\frac{R_{o}}{R}\right)^{4}\right]\!.\end{eqnarray}$$
Equation (2.18) is consistent with previous results for nucleation of voids in elastic solids (Gent & Tompkins Reference Gent and Tompkins1969; Horgan & Polignone Reference Horgan and Polignone1995), although those studies neglected inertia and viscosity and focused primarily on quasi-static behaviour, i.e.
${\dot{R}}=0$
and
$\ddot{R}=0$
in (2.17).
Compressible effects are crucial during violent bubble collapse (Prosperetti & Lezzi Reference Prosperetti and Lezzi1986; Fuster, Dopazo & Hauke Reference Fuster, Dopazo and Hauke2011). The approach of Keller & Miksis (Reference Keller and Miksis1980) incorporating acoustic radiation through liquid compressibility is implemented here, following Yang & Church (Reference Yang and Church2005) where it is assumed that, for the stress integral calculation, the divergence-free velocity is not modified by first-order compressible effects. The compressible extension to (2.17) yields
$$\begin{eqnarray}\displaystyle \left(1-\frac{{\dot{R}}}{c}\right)R\ddot{R}+\frac{3}{2}\left(1-\frac{{\dot{R}}}{3c}\right){\dot{R}}^{2} & = & \displaystyle \left(1+\frac{{\dot{R}}}{c}\right)\left[\frac{p_{B}-p_{\infty }(t)}{{\it\rho}}-\frac{4{\it\nu}_{L}{\dot{R}}}{R}-\frac{2S}{{\it\rho}R}-E\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{R}{{\it\rho}c}\frac{\text{d}}{\text{d}t}(\,p_{a}-p_{\infty }),\end{eqnarray}$$
where
$c$
is the constant sound speed of the undisturbed surrounding medium and
$$\begin{eqnarray}p_{a}(t)=p_{B}-\frac{2S}{R}-\frac{4{\it\nu}_{L}{\dot{R}}}{R}-E.\end{eqnarray}$$
While the focus in this section is on a neo-Hookean medium, similar expressions can be derived for other constitutive relations (cf. appendix A for elastic terms corresponding to a Mooney–Rivlin strain-energy function and ‘linear’ elastic models).
Figure 3. Absolute value of the dimensionless elastic force acting on the bubble (a, log–log) with a zoomed in view near
$R/R_{o}=1$
(b, linear–linear). For
$R/R_{o}<1$
, the force is negative.
Table 1. Summary of the elastic terms in the Rayleigh–Plesset/Keller–Miksis equations for the one-parameter constitutive models.
3. The role of nonlinear elasticity
3.1. The elastic force
The role of elasticity is investigated by determining the corresponding force under various limiting circumstances and for different constitutive models, to understand its main effects on the dynamics. The elastic terms correspond to a restoring force or ‘spring’; at any point along the interface, an infinitesimal force element goes as
${\sim}R^{2}E(R)$
. Table 1 summarizes the findings for the neo-Hookean, Yang & Church (Reference Yang and Church2005) and linear models, and Figure 3 provides a visual representation of the absolute value of the normalized elastic force. Since the force is dimensionless, the plots are valid for any
${\it\eta}$
. All the models reduce to zero for
${\it\eta}=0$
or zero strains, as expected. Unlike surface tension, the elastic terms can be positive (
$R<R_{o}$
, i.e. during collapse) or negative (
$R>R_{o}$
, i.e. for growth).
In general, the models with the largest elastic effects are the models that include a stronger nonlinearity, a behaviour that is most clear during collapse (
$R\ll R_{o}$
). By its sign and dependence on
$(R_{o}/R)^{n}$
, with
$n>0$
, the elastic force acts to reduce the acceleration of the bubble wall during the last stages of collapse; in other words, the elastic force (spring) is ‘pulling’ the bubble wall back towards the equilibrium condition, and increases as the bubble radius decreases for the nonlinear models. Given the asymptotic behaviour as
$R\ll R_{o}$
, this force is smallest in the linear elastic case, which in fact tends to a constant value as the radius decreases, with a slope given by the exponent of
$R_{o}/R$
in the fourth column in table 1. For other models, nonlinearity leads to an increase in this restoring force as the bubble size decreases. The neo-Hookean model exhibits the largest force. The rate at which the force of the nonlinear models increases depends on the exponent of
$R_{o}/R$
.
For small departures from the initial radius,
$R\approx R_{o}+{\it\epsilon}(t)$
with
${\it\epsilon}\ll R_{o}$
, perturbation analysis shows that the elastic terms reduce to
$-4{\it\eta}{\it\epsilon}$
to first order for all models, i.e. small-amplitude oscillations in a linear elastic medium. For larger departures from
$R_{o}$
, the behaviour depends on the model, though all are close. For
$R\lesssim R_{o}$
, the Yang & Church (Reference Yang and Church2005) model exhibits the largest force, followed by the neo-Hookean model. The behaviour for
$R\gtrsim R_{o}$
shows that the force is largest for the linear case, followed by the neo-Hookean model. In the limit
$R\gg R_{o}$
, elasticity hinders growth at the same rate for all models, although the coefficients are slightly different, with the neo-Hookean model being the largest.
As a comparison, the elastic term for the Mooney–Rivlin model goes as
${\it\eta}_{1}[5-4(R_{o}/R)-(R_{o}/R)^{4}]/2-{\it\eta}_{2}[2(R/R_{o})-(R_{o}/R)^{2}-1]$
, with
${\it\eta}={\it\eta}_{1}+{\it\eta}_{2}$
(see appendix A). Since
${\it\eta}_{2}>0$
, the total force is less. The behaviour as
$R\ll R_{o}$
is similar to that of the neo-Hookean model; for large growth however, the behaviour is governed by the term multiplied by
${\it\eta}_{2}$
. The Mooney–Rivlin model is not included in the comparison, because of the need to define a second parameter. For the viscoelastic shell of Liu et al. (Reference Liu, Sugiyama, Takagi and Matsumoto2012), the elastic term goes as
$2{\it\eta}[1-(R_{o}/R)^{6}]$
. For this case, the force at collapse is large (
${\sim}1/R^{4}$
); for large growth, the behaviour is similar to the linear case.
3.2. Natural and peak frequencies of the bubble
The frequency response of the bubble is of particular interest for ultrasound applications. To investigate this aspect, a bubble initially at equilibrium is subjected to a pressure waveform of frequency
${\it\omega}$
:
$$\begin{eqnarray}p_{\infty }(t)=p_{\infty }[1+g(t)],\end{eqnarray}$$
where
$0<|g|\ll 1$
. For simplicity, all forms of damping (viscous and compressible) are included in an effective viscosity. For this small-amplitude perturbation analysis, the same result is obtained with all models. The bubble is expected to oscillate about an equilibrium radius, with
$$\begin{eqnarray}R(t)=R_{eq}[1+f(t)],\end{eqnarray}$$
where
$0<|\,f|\ll 1$
and
$R_{eq}$
is not necessarily
$R_{o}$
. Substitution into (2.17) yields, to first order,
$$\begin{eqnarray}\displaystyle & & \displaystyle \ddot{f}+\frac{4{\it\nu}_{L}}{R_{eq}^{2}}{\dot{f}}+\frac{1}{R_{eq}^{2}}\left\{3k\frac{p_{Go}}{{\it\rho}}\left(\frac{R_{o}}{R_{eq}}\right)^{3k}-\frac{2S}{{\it\rho}R_{eq}}+\frac{2{\it\eta}}{{\it\rho}}\left(\frac{R_{o}}{R_{eq}}\right)\left[1+\left(\frac{R_{o}}{R_{eq}}\right)^{3}\right]\right\}f+\frac{p_{\infty }}{{\it\rho}R_{eq}^{2}}g\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{{\it\rho}R_{eq}^{2}}\left\{-(\,p_{\infty }-p_{v})+p_{Go}\left(\frac{R_{o}}{R_{eq}}\right)^{3k}-\frac{2S}{R_{eq}}+\frac{{\it\eta}}{2}\left[-5+4\frac{R_{o}}{R_{eq}}+\left(\frac{R_{o}}{R_{eq}}\right)^{4}\right]\right\}.\end{eqnarray}$$
The left-hand side of the equation contains all the perturbation terms, i.e. all terms that are time dependent and
$O(\,f,g)$
. The entire right-hand side is
$O(1)$
and corresponds to equilibrium conditions. Two main cases are of interest: finite
$g>0$
where
$R_{eq}\neq R_{o}$
, which corresponds to Rayleigh collapse, and
$R_{eq}=R_{o}$
. The former was investigated using perturbation methods by Hua & Johnsen (Reference Hua and Johnsen2013) and is thus not pursued further here. The latter simplifies greatly:
$$\begin{eqnarray}\ddot{f}+\frac{4{\it\nu}_{L}}{R_{eq}^{2}}{\dot{f}}+\frac{1}{R_{eq}^{2}}\left(3k\frac{p_{Go}}{{\it\rho}}-\frac{2S}{{\it\rho}R_{eq}}+\frac{4{\it\eta}}{{\it\rho}}\right)f=-\frac{p_{\infty }}{{\it\rho}R_{eq}^{2}}g.\end{eqnarray}$$
To first order, the damping is not affected by elasticity, but the oscillation frequency is. The peak frequency is given by
$$\begin{eqnarray}{\it\omega}_{p}=\sqrt{\frac{1}{{\it\rho}R_{o}^{2}}\left[3k(\,p_{\infty }-p_{v})+(3k-1)\frac{2S}{R_{o}}\right]-\frac{8{\it\nu}_{L}^{2}}{R_{o}^{4}}+\frac{4{\it\eta}}{{\it\rho}R_{o}^{2}}}\end{eqnarray}$$
and the natural frequency is obtained by setting
${\it\nu}_{L}=0$
. The frequency increases with increasing elasticity, and on setting
${\it\eta}=0$
reduces to the Minnaert frequency, consistent with Hua & Johnsen (Reference Hua and Johnsen2013).
Frequency response curves can be obtained following Lauterborn (Reference Lauterborn1976). The curves have the form
$$\begin{eqnarray}\frac{R_{max}-R_{eq}}{R_{eq}}=f\left(\frac{{\it\omega}}{{\it\omega}_{n}}\right),\end{eqnarray}$$
where
${\it\omega}_{n}$
is the natural frequency of the bubble with no elasticity,
$R_{eq}=R_{o}$
and
$R_{max}$
is the maximum radius of the bubble during its steady-state oscillations.
Figure 4 shows the frequency response contour map of a bubble in a neo-Hookean medium (with viscosity 1 cP) as the shear modulus and driving frequency are varied. A
$100\times 100$
point parameter space was used, and the amplitude was defined as the average amplitude in the latter half of fifteen acoustic cycles to remove the effect of the initial transient, with the bubble starting from rest. For low shear modulus, one main resonant peak is observed close to the bubble natural frequency, with many subharmonics occurring at lower frequencies. The spacing between these subharmonics decreases with decreasing frequency; at high frequencies, the amplitude of the oscillations is small. The resonance map does not significantly change until the shear modulus reaches
${\sim}1~\text{kPa}$
. At this point, the effect of elasticity is to shift the resonance peak and dominant subharmonics to the right to higher frequencies.
Figure 4. Frequency response contour map: dimensionless amplitude of a bubble in a neo-Hookean medium. The driving-pressure amplitude is
$0.8$
bar.
In Figure 5, the resonance curve is shown for
${\it\eta}=100~\text{kPa}$
, for the Newtonian, neo-Hookean, Yang & Church (Reference Yang and Church2005) and linear elastic media. For the neo-Hookean case, this plot corresponds to a horizontal cross-section of Figure 4. For elastic media, the frequency response curves are shifted to the right since
${\it\omega}_{n}$
is larger when
${\it\eta}\neq 0$
, as predicted by (3.5). Moreover, the resonance peaks are weakened because of elasticity. Since the peaks occur for radii
$R\gtrsim R_{o}$
, the elastic force is largest for the linear model, followed by the neo-Hookean model. The Yang & Church (Reference Yang and Church2005) model exhibits the weakest force in that regime and thus predicts larger resonance amplitudes, and at lower resonance frequencies, than the other models. Conversely, the linear elastic model predicts the lowest amplitude resonance responses, and at the highest frequencies. The Newtonian curve qualitatively corresponds to that obtained in Lauterborn (Reference Lauterborn1976). The present results are consistent with § 3.1.
Figure 5. Frequency response curves for
${\it\eta}=100~\text{kPa}$
and a pressure–wave amplitude of 0.8 bar.
This analysis can also be used to assess the presence of stable nuclei in an elastic medium. Elasticity modifies the Blake radius (Blake Reference Blake1949), e.g.
$$\begin{eqnarray}p_{v}-p_{\infty }>\frac{4\displaystyle \frac{S}{R_{eq}}+{\it\eta}}{3}\end{eqnarray}$$
for the isothermal case. The additional elastic term in this stability condition is always positive. Unsurprisingly, elasticity stabilizes the bubble equilibrium since the natural frequency exists in a broader range of values for
$p_{v}-p_{\infty }$
. Figure 6 shows the peak frequency as a function of
$R_{eq}$
for
${\it\eta}=0$
,
$10^{5}$
and
$10^{6}~\text{Pa}$
under isothermal conditions (
$k=1$
) and with a tension
$p_{v}-\bar{p}_{\infty }=1$
bar. The area under each curve corresponds to the stability region. As expected, the larger the shear modulus, the larger the stability region, such that smaller nuclei may exist. With a shear modulus of 1 MPa, the smallest equilibrium radius possible is approximately 50 nm, compared with 500 nm for zero shear modulus.
Figure 6. Peak frequency
${\it\omega}_{p}$
as a function of
$R_{eq}$
.
3.3. Minimum radius
The minimum radius and the maximum velocity during collapse are important quantities related to potential damage incurred to the surroundings. To estimate the role of elasticity on these quantities, analysis is performed for Rayleigh collapse following Brennen (Reference Brennen1995), in which viscous and compressible effects are neglected. Here, we include elasticity using the neo-Hookean model under isothermal conditions. From (2.13) and (2.17), the Rayleigh–Plesset equation can be rewritten
$$\begin{eqnarray}\displaystyle {\dot{R}}^{2} & = & \displaystyle \frac{2[p_{v}-p_{\infty }-{\rm\Delta}p]}{3{\it\rho}_{L}}\left[1-\left(\frac{R_{o}}{R}\right)^{3}\right]+\frac{2p_{Go}}{{\it\rho}_{L}}\left(\frac{R_{o}}{R}\right)^{3}\log \left(\frac{R}{R_{o}}\right)-\frac{2S}{{\it\rho}R}\left[1-\left(\frac{R_{o}}{R}\right)^{2}\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{{\it\eta}}{3{\it\rho}}\left[-5+6\left(\frac{R_{o}}{R}\right)+2\left(\frac{R_{o}}{R}\right)^{3}-3\left(\frac{R_{o}}{R}\right)^{4}\right]\!,\end{eqnarray}$$
where
${\rm\Delta}p$
is the instantaneous pressure increase. At collapse,
$R\ll R_{o}$
, such that
$$\begin{eqnarray}{\dot{R}}\rightarrow -\!\left(\frac{R_{o}}{R}\right)^{3/2}\sqrt{\frac{2[p_{\infty }+{\rm\Delta}p-p_{v}]}{3{\it\rho}_{L}}+\frac{2S}{{\it\rho}R_{o}}+\frac{2p_{Go}}{{\it\rho}_{L}}\log \left(\frac{R}{R_{o}}\right)+\frac{{\it\eta}}{{\it\rho}}\left(\frac{2}{3}-\frac{R_{o}}{R}\right)}.\end{eqnarray}$$
Since
$R<R_{o}$
for this problem, elasticity always reduces the velocity during collapse.
By defining
$A=2p_{Go}/{\it\rho}_{L}$
,
$B=-{\it\eta}/{\it\rho}$
and
$C=(2[p_{\infty }+{\rm\Delta}p-p_{v}]/3{\it\rho}_{L})+(2S/{\it\rho}R_{o})+(2{\it\eta}/3{\it\rho})$
, the minimum radius can be written in terms of the Wright function
$W(x)$
:
$$\begin{eqnarray}\frac{R_{min}}{R_{o}}=\frac{\exp [W(\log (-B/A)+(C/A))]}{\exp (C/A)}.\end{eqnarray}$$
It should be noted that if
${\it\eta}=0$
,
$R_{min}=R_{o}\exp (A/C)$
, which corresponds to the minimum radius in a Newtonian medium. The adiabatic case can be solved following the same approach given
$k$
.
4. Bubble response to various driving pressures
Comparisons between the different models are presented to illustrate the effect of elasticity on the bubble dynamics and verify the analysis in the previous section using three problems in which the driving pressures and initial/equilibrium radii are different, with relevance to biomedical applications. For all runs,
$c=1500~\text{m}~\text{s}^{-1}$
,
$p_{v}=2300~\text{Pa}$
,
$S=0.05~\text{N}~\text{m}^{-1}$
,
${\it\nu}_{L}=10^{-6}~\text{m}^{2}~\text{s}^{-1}$
and
${\it\rho}=1000~\text{kg}~\text{m}^{-3}$
. The ambient pressure is
$p_{\infty }(0)=1$
bar. The linear shear modulus is
${\it\eta}=100~\text{kPa}$
, representative of soft tissue (Madsen, Sathoff & Zagzebski Reference Madsen, Sathoff and Zagzebski1983; Wells & Liang Reference Wells and Liang2011). The Newtonian, neo-Hookean, Yang & Church (Reference Yang and Church2005) and linear models are considered. The initial radius
$R_{o}=10~{\rm\mu}\text{m}$
is the radius at time
$t=0$
, while the equilibrium radius
$R_{eq}$
is the radius achieved as
$t\rightarrow \infty$
, which need not be
$R_{o}$
. These initial-value problems are solved using a fifth-order accurate Runge–Kutta–Verner scheme (Verner Reference Verner1978).
4.1. Step increase in pressure (Rayleigh collapse)
First, the traditional Rayleigh collapse problem (Rayleigh Reference Rayleigh1917) is considered, in which the surrounding pressure
$p_{\infty }(0)$
is instantaneously increased to
$p_{\infty }(0)+{\rm\Delta}p$
for
$t>0^{+}$
, with
${\rm\Delta}p=100$
bar. This problem can be considered as an idealization of shock-induced bubble collapse in which the sound speed in the surroundings is infinite. In this problem,
$R_{o}=R_{eq}$
. The initial velocity is
${\dot{R}}_{o}=(\,p_{v}-p_{\infty }(0))/{\it\rho}c_{\infty }$
following Plesset & Prosperetti (Reference Plesset and Prosperetti1977).
Figures 7 and 8 show histories of the bubble radius and wall velocity. The equilibrium radius is different for each model since the pressure increase
${\rm\Delta}p$
changes the final elastic stress state (Hua & Johnsen Reference Hua and Johnsen2013). Considering the neo-Hookean model, the initial conditions prescribe
$p_{Go}=(\,p_{\infty }-p_{v})+2S/R_{o}$
. The equilibrium radius
$R_{eq}$
is given implicitly by
$$\begin{eqnarray}-\!(\,p_{\infty }+{\rm\Delta}p-p_{v})+p_{Go}\left(\frac{R_{o}}{R_{eq}}\right)^{3k}-\frac{2S}{R_{eq}}+\frac{{\it\eta}}{2}\left[4\frac{R_{o}}{R_{eq}}+\left(\frac{R_{o}}{R_{eq}}\right)^{4}-5\right]=0.\end{eqnarray}$$
Since the elastic term is different for each model, different equilibrium radii are achieved. As expected, the linear elastic model produces the smallest deviation from the Newtonian case, followed by the Yang & Church (Reference Yang and Church2005) and neo-Hookean models respectively.
Figure 7. History of the bubble radius following a step increase in pressure for different constitutive models (
${\rm\Delta}p=100$
bar). Time is non-dimensionalized by the Rayleigh collapse time
$t_{RC}$
.
Figure 8. Detailed view of the history of the bubble radius (a) and velocity (b) following a step increase in pressure at the first collapse. Time is non-dimensionalized by the Rayleigh collapse time
$t_{RC}$
.
The minimum radius and maximum velocity are important because the shocks emitted upon collapse, and thus the potential damage inflicted by cavitation, directly depend on these quantities. The smallest minimum radius, highest velocity and shortest collapse time are achieved in the Newtonian medium, followed in order by the linear elastic, Yang & Church (Reference Yang and Church2005) and neo-Hookean media, consistent with § 3.1. For the more nonlinear elastic models, the violence of the collapse is reduced, as exhibited by the increase in minimum radius and collapse time, and decrease in velocity. At this specific value of
${\rm\Delta}p$
, the velocity at collapse is more than five times smaller in the neo-Hookean medium than in a Newtonian liquid, thus reducing potential damage due to shocks.
After the first collapse, the largest relative rebound is achieved in the neo-Hookean material, followed by Yang & Church (Reference Yang and Church2005), linear elastic and Newtonian. This behaviour is inversely related to the minimum radius, as a collapse to a smaller radius dissipates more energy by acoustic radiation (Johnsen & Colonius Reference Johnsen and Colonius2009). Thus, the following rebounds are smaller and the oscillation period is shorter, consistent with § 3.2. This observation has important implications for the biomedical field, since larger oscillations lead to larger strains and higher strain rates.
4.2. Bubble out of equilibrium
Next, the growth and collapse of a bubble whose initial radius is set out of equilibrium keeping
$p_{\infty }$
constant are considered (Flynn Reference Flynn1975). This problem is a simpler case of cavitation-bubble growth/collapse and collapse/growth due to the passage of a transient tension/compression, set up in such a way that the transient does not affect the dynamics; the growth problem can be considered to model laser-induced cavitation (Lauterborn & Kurz Reference Lauterborn and Kurz2010). The initial state out of equilibrium corresponds to an initial internal bubble pressure that is lower or higher than the ambient pressure, i.e.
$R_{o}\neq R_{eq}$
, depending on the selected initial radius. For a prescribed equilibrium radius and in the absence of elasticity, the partial pressure (or number of moles) of non-condensible gas is given by
$p_{Geq}=(\,p_{\infty }-p_{v})+2S/R_{eq}$
. In practice, this problem is started by increasing or decreasing
$p_{Go}$
corresponding to a change in the initial radius
$R_{o}$
. Thus, collapse or growth follows.
Figures 9 and 10 show the history of the bubble radius for different ratios of the initial radius to the equilibrium radius, corresponding to
$R_{o}>R_{eq}^{N}$
(initial collapse) and
$R_{o}<R_{eq}^{N}$
(initial growth) respectively, where
$R_{eq}^{N}$
stands for the equilibrium radius obtained in the Newtonian case;
${\rm\Delta}p$
is the initial pressure difference between the bubble and the surroundings.
For
$R_{o}>R_{eq}^{N}$
(initial collapse), elasticity significantly inhibits collapse, more so than in the Rayleigh collapse of § 4.1. The linear model exhibits the highest force in this regime, such that the oscillation amplitude is smallest; the amplitude in the neo-Hookean medium is very close. For
$R_{o}<R_{eq}^{N}$
(initial growth), elasticity inhibits growth, with the Yang & Church (Reference Yang and Church2005) model exhibiting the strongest effect. The oscillations amplitude in the neo-Hookean medium are very close. These results are consistent with the analysis in § 3.1. Based on the inertial cavitation threshold criterion by Flynn (Reference Flynn1975), it is evident that, despite satisfying
$R_{o}/R_{eq}\gtrsim 2$
, the oscillations in viscoelastic media are not highly nonlinear.
4.3. Harmonic forcing
The response of a bubble to a sinusoidally varying pressure waveform is relevant to acoustic cavitation in biomedical applications. For this problem,
$R_{o}=R_{eq}$
. Figure 11 shows the response of a bubble to pressure waveforms of different driving frequencies
$f$
. The driving frequencies are in the range
$f=0.1{-}3~\text{MHz}$
with an amplitude of 1 MPa, for clinical relevance and comparisons with Yang & Church (Reference Yang and Church2005). Overall, as the frequency is increased, smaller maximum radii are achieved, consistent with figure 4. The results strongly depend on the frequency.
-
(i) For
$f<200~\text{kHz}$
, the bubble oscillates with the waveform at high amplitude. Elasticity restrains the bubble oscillations (e.g. maximum radius), and subharmonics are observed in the neo-Hookean medium. -
(ii) For
$200~\text{kHz}\lesssim f\lesssim 1000~\text{kHz}$
, different behaviours are observed depending on the model, particularly with increasing frequency. As explained in § 3.2, elasticity modifies the bubble natural frequency. For instance, the natural frequency with no elasticity is
$f_{n}=282~\text{kHz}$
, compared with
$f_{n}=425~\text{kHz}$
for a shear modulus of 100 kPa. Thus, a bubble subjected to harmonic forcing of the order of these frequencies may undergo large oscillations, with a behaviour different for each constitutive model. -
(iii) For
$f>1000~\text{kHz}$
, all models exhibit the same behaviour over the course of the first collapse, at which point the Newtonian solution departs from the elastic solutions due to the action of the elastic force. All elastic models reduce to the same solution, as oscillations about the time-varying steady state are small. This behaviour is expected, since the strains achieved remain small.
Figure 9. History of the bubble radius for the bubble initially out of equilibrium for different constitutive models (
$R_{o}/R_{eq}^{N}=2$
corresponds to
${\rm\Delta}p=-0.83$
bar). Time is non-dimensionalized by the Rayleigh collapse time
$t_{RC}$
.
Figure 10. History of the bubble radius for the bubble initially out of equilibrium for different constitutive models (
$R_{o}/R_{eq}^{N}=0.5$
corresponds to
${\rm\Delta}p=7.24$
bar). Time is non-dimensionalized by the Rayleigh collapse time
$t_{RC}$
.
Figure 11. History of the bubble radius under harmonic forcing with
$f=100$
, 282, 425, 1000 and 3000 kHz (a–e) and
${\it\eta}=100~\text{kPa}$
.
5. Conclusions
A theoretical framework to simulate spherical bubble dynamics in a viscoelastic medium with nonlinear elasticity (e.g. soft tissue) is presented. To account for the fact that tissue usually relaxes to an original configuration after loading, a Kelvin–Voigt viscoelastic model is chosen, with a nonlinear spring. Based on modern continuum mechanics formalism, different strain-energy functions are included to describe the nonlinear elastic component, which can then be inserted appropriately into a Rayleigh–Plesset/Keller–Miksis equation. One of the main contributions of this article lies in the derivation of these terms for several models: nonlinear neo-Hookean and Mooney–Rivlin, and linear elasticity. The model used previously by Yang & Church (Reference Yang and Church2005) was recovered using this approach. The response of each model is different for the problems of interest; these differences can be explained analytically by examining the elastic force. The following observations are made concerning the physics.
-
(i) In agreement with previous work on cavitation in viscoelastic media, elasticity overall acts to reduce the violence of collapse and growth, compared with the behaviour in a Newtonian medium. The elastic terms depend on the departure from the original configuration, with different coefficients and exponents depending on the constitutive model. These terms represent an elastic restoring force, i.e. a spring, as expected.
-
(ii) The bubble response may vary significantly depending on the constitutive model.
-
(iii) For large departures from the original configuration, models with stronger nonlinearity produce stronger forces, particularly at collapse. This force increases with decreasing radius. As a result, larger minimum radii and lower velocities are obtained at collapse. These quantities are important with regard to potential damage inflicted to the surrounding medium.
-
(iv) For small departures from the original configuration, all models reduce to the same response, as expected. For such cases, elasticity modifies the oscillation properties: the peak and natural frequencies increase, thus shifting the frequency response, and the amplitude is reduced, especially for the more nonlinear models. Furthermore, smaller stable nuclei are expected to be found in such media.
The proposed approach can readily be extended to other strain-energy functions. The findings of the present work have potentially important implications for cavitation research in biomedical applications, as several experimentally measurable quantities (frequency response, minimum radius) are reported. Overall, elasticity is expected to reduce typical cavitation damage due to shocks produced during violent collapse. Furthermore, the resulting bubble dynamics equation could be used for characterization of viscoelastic media (rheometry), e.g. determination of an appropriate model through fits, and measurement of viscoelastic properties, e.g. viscosity and shear modulus, even for complex constitutive models and high frequencies. At this time, several factors limit direct application of the model, e.g. constitutive modelling of tissue, inhomogeneities, compressibility, thermal and mass transport. To address these issues, development of numerical methods and models for cavitation and shock waves in tissue-like media, and single-bubble experiments in tissue phantoms are currently underway.
Acknowledgements
The authors thank A. Maxwell, E. Vlaisavljevich, Z. Xu, D. Miller, J. Freund and A. Wineman for helpful conversations. This work was supported in part by NSF grant number CBET 1253157 and NIH grant number 1R01HL110990-01A1.
Appendix A. Different constitutive models
The same approach as that used for the neo-Hookean model can be followed to derive the elastic terms for other constitutive models relevant for cavitation in tissue. Here, we show how to obtain the elastic term
$E$
in (2.17) and (2.19) for the Mooney–Rivlin strain-energy function and linear elastic models.
A.1. Nonlinear elasticity: Mooney–Rivlin model
The Mooney–Rivlin model (Mooney Reference Mooney1940; Rivlin Reference Rivlin1948), the simplest two-parameter nonlinear elastic model, is considered. The corresponding strain-energy function is
$$\begin{eqnarray}{\hat{W}}=\frac{{\it\eta}_{1}}{2}(I_{1}-3)+\frac{{\it\eta}_{2}}{2}(I_{2}-3),\quad {\it\eta}={\it\eta}_{1}+{\it\eta}_{2},\quad {\it\eta}_{1},{\it\eta}_{2}>0,\end{eqnarray}$$
where
${\it\eta}$
is the linear shear modulus and
$I_{2}$
is the second strain invariant. The Cauchy stress tensor is then given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64F}_{\boldsymbol{e}}=\left[\begin{array}{@{}ccc@{}}\displaystyle {\it\eta}_{1}\left(\frac{r_{o}}{r}\right)^{4}-2{\it\eta}_{2}\left(\frac{r_{o}}{r}\right)^{2}-\mathscr{P} & 0 & 0\\ 0 & {\it\eta}_{1}\displaystyle \left(\frac{r}{r_{o}}\right)^{2}-{\it\eta}_{2}\displaystyle \left[\left(\frac{r}{r_{o}}\right)^{4}-\left(\frac{r_{o}}{r}\right)^{2}\right]-\mathscr{P} & 0\\ 0 & 0 & {\it\eta}_{1}\displaystyle \left(\frac{r}{r_{o}}\right)^{2}-{\it\eta}_{2}\displaystyle \left[\left(\frac{r}{r_{o}}\right)^{4}-\left(\frac{r_{o}}{r}\right)^{2}\right]-\mathscr{P}\end{array}\right]\!, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
and the relation between the actual pressure and the pseudo-pressure is
$$\begin{eqnarray}\mathscr{P}=p+\frac{{\it\eta}_{1}}{3}\left[\left(\frac{r_{o}}{r}\right)^{4}+2\left(\frac{r}{r_{o}}\right)^{2}\right]-\frac{2{\it\eta}_{2}}{3}\left[\left(\frac{r}{r_{o}}\right)^{4}+2\left(\frac{r_{o}}{r}\right)^{2}\right]\!.\end{eqnarray}$$
The elastic terms for a Mooney–Rivlin strain-energy function are given by
$$\begin{eqnarray}E_{\mathit{MR}}(t)=\frac{{\it\eta}_{1}}{2}\left[5-4\left(\frac{R_{o}}{R}\right)-\left(\frac{R_{o}}{R}\right)^{4}\right]-{\it\eta}_{2}\left[2\left(\frac{R}{R_{o}}\right)-\left(\frac{R_{o}}{R}\right)^{2}-1\right]\!.\end{eqnarray}$$
If
${\it\eta}_{2}=0$
, the neo-Hookean model is recovered, thus verifying that the Mooney–Rivlin model is a refinement of the neo-Hookean model.
A.2. Linear elasticity: generalized Hooke’s law
The large strains expected in the problems of interest motivate the use of nonlinear elastic models. However, it is useful to consider the limiting case of infinitesimal strains, i.e. linear elasticity. The same assumptions as those made until now hold, except that the generalized Hooke’s law (Hoger & Johnson Reference Hoger and Johnson1995) is used for the elastic stress:
$$\begin{eqnarray}\unicode[STIX]{x1D64F}_{\boldsymbol{e}}=2{\it\eta}\unicode[STIX]{x1D640}-\mathscr{P}\unicode[STIX]{x1D644},\end{eqnarray}$$
where
$\unicode[STIX]{x1D640}$
is the linear strain tensor,
${\it\eta}$
is the shear modulus and
$\mathscr{P}(r,t)$
is a pseudo-pressure. The linear strain can be expressed in terms of the displacement gradient
$\unicode[STIX]{x1D643}=(\unicode[STIX]{x1D643}+\unicode[STIX]{x1D643}^{\text{T}})/2=\unicode[STIX]{x1D640}$
, since
$\unicode[STIX]{x1D643}$
is diagonal. Since
$\unicode[STIX]{x1D641}=\unicode[STIX]{x1D644}+\unicode[STIX]{x1D643}$
, the Cauchy stress tensor for a linear elastic medium is given by
$$\begin{eqnarray}\unicode[STIX]{x1D64F}=\left[\begin{array}{@{}ccc@{}}\displaystyle 2{\it\eta}\left[\left(\frac{r_{o}}{r}\right)^{2}-1\right]-\mathscr{P} & 0 & 0\\ 0 & 2{\it\eta}\displaystyle \left[\left(\frac{r}{r_{o}}\right)-1\right]-\mathscr{P} & 0\\ 0 & 0 & 2{\it\eta}\displaystyle \left[\left(\frac{r}{r_{o}}\right)-1\right]-\mathscr{P}\end{array}\right]\!,\end{eqnarray}$$
and the relation between the actual pressure and the pseudo-pressure is
$$\begin{eqnarray}\mathscr{P}=p+\frac{2{\it\eta}}{3}\text{tr}(\unicode[STIX]{x1D640})=p+\frac{2{\it\eta}}{3}\left[\left(\frac{r_{o}}{r}\right)^{2}+2\left(\frac{r}{r_{o}}\right)-3\right]\!.\end{eqnarray}$$
The elastic terms for the generalized Hooke’s law are given by
$$\begin{eqnarray}E_{L}=2{\it\eta}\left[1-\left(\frac{R_{o}}{R}\right)^{2}\right]\!.\end{eqnarray}$$
A.3. Linear elasticity: infinitesimal strains (Yang & Church Reference Yang and Church2005)
In Yang & Church (Reference Yang and Church2005), an equation similar to that obtained in § A.2 is derived, but with different coefficients and powers. We show here how to recover this result. By retaining the first-order terms in the Taylor series expansions of
$(r/r_{o})^{n}$
from (2.15), and assuming small strains (i.e.
$|(R^{3}-R_{o}^{3})/r^{3}|\ll 1$
), the
$rr$
component of the Cauchy stress is given by
$$\begin{eqnarray}-\!\frac{4{\it\eta}}{3}\frac{R^{3}-R_{o}^{3}}{r^{3}}-\mathscr{P}.\end{eqnarray}$$
The elastic contribution is then obtained by evaluating
$$\begin{eqnarray}E_{\mathit{YC}}=-3\int _{R}^{\infty }\frac{{\it\tau}_{rr}}{r}\text{d}r=\frac{4{\it\eta}}{3}\left[1-\left(\frac{R_{o}}{R}\right)^{3}\right]\!.\end{eqnarray}$$
