2.1. Interface energy model

In the existing literature (Helfrich Reference Helfrich1973; Rangamani et al. Reference Rangamani, Benjamini, Agrawal, Smit, Steigmann and Oster2013; Argudo et al. Reference Argudo, Bethel, Marcoline and Grabe2016; Bian, Litvinov & Koumoutsakos Reference Bian, Litvinov and Koumoutsakos2020), modelling of cell membranes and lipid bi-layers is performed using Helfrich curvature energy which takes care of the bending energy through Gaussian curvature and spontaneous curvature. However, one cannot directly use Helfrich energy to model interfaces, which requires a generic framework to identify the structure of interface energy in terms of material symmetries. Also, it is important to note that these models are used to study the initial natural states of micro- or nanoscale structures. On the other hand, we are interested in modelling encapsulated microbubbles with the effect of residual fields developed due to the interface energy. Therefore, we follow the generic continuum framework developed in the previous work (Fung Reference Fung1977; Steigmann & Ogden Reference Steigmann and Ogden1999; Gao et al. Reference Gao, Huang, Qu and Fang2014) to incorporate the interface energy effects using the general notions of differential geometry (Itskov Reference Itskov2019; Steinmann Reference Steinmann2015).

We assume that the bubble interface is a deformable mathematical surface compatible with the deformation of the bulk. We assume that there exist a bulk stress-free fictitious configuration (FC) of the bubble as shown in figure 1(a) with inner and outer radii, $R_{e_1}$ and $R_{e_2}$, respectively. Even in the absence of external load this FC is unstable due to the presence of excess interface energy and the corresponding non-zero interface stress, which can force the bubble to attain a stable self-equilibrium state through expansion or shrinkage. This stable self-equilibrium state with bulk residual stress developed due to heterogeneous interface energy is shown in figure 1(b) with inner and outer radii $R_{10}$ and $R_{20}$, respectively. We also call this configuration a natural configuration (NC). Therefore, the bulk residual stress field in the bubble can be described as the state of self-equilibrium stress field due to the residual interface stresses. When excited by an external ultrasound pressure field, figure 1(c) represents the radial dynamic configuration (DC) of the same bubble with time-dependent inner and outer radii, $R_1(t)$ and $R_2(t)$, respectively.

Figure 1.

Schematic representation of a gas filled encapsulated microbubble suspended in a viscous medium at three different configurations. (a) Fictitious (bulk stress-free) configuration (FC) of the bubble with possible expansion or shrinkage to reach the natural (static equilibrium) configuration (NC) in (b), and (c) Dynamic configuration (DC) of the bubble at an instant of time $t$. (a) Fictitious configuration. (b) Natural configuration. (c) Dynamic configuration.

The approach presented above allows accounting for the influence of excess interfacial free energy and the corresponding non-zero residual interface stress on the mechanical behaviour of EBs. This FC is in agreement with the imaginary ‘fictitious stress-free configuration’ suggested by Gao et al. (Reference Gao, Huang, Qu and Fang2014). It is important to note that this FC may not exist physically. However, this assumption will help in calculating the bulk residual stress field induced by interfacial energy. Though many researchers have ignored the role of interfacial energy induced bulk residual stress, it should be noted that when we consider the effect of interface energy, the bulk residual stress field induced by the interfacial energy always becomes an important characteristic to account for and, hence, cannot be ignored. In very recent work, Chabouh et al. (Reference Chabouh, Dollet, Quilliet and Coupier2021) have also highlighted the importance of this stress-free configuration for developing dynamic models for EBs.

Let $(x^{1},x^{2},x^{3})$ be the three-dimensional Cartesian coordinate system with unit basis vector triad as $(\boldsymbol {e}_1,\boldsymbol {e}_2,\boldsymbol {e}_3)$. Let $(\rho,\phi,\theta )$ be the principal spherical polar coordinates with unit basis vector triad $(\boldsymbol {e}_\rho,\boldsymbol {e}_\phi,\boldsymbol {e}_\theta )$. Here, $\rho$ is the radial coordinate from the centre of the bubble, $\phi$ is the polar angle formed with the $x^{3}$-axis and $\theta$ is the azimuthal angle measured about the $x^3$-axis in the counterclockwise direction from the $x^1$-axis. The relation between the Cartesian and the spherical polar coordinates is given by

(2.1)\begin{equation} \left. \begin{aligned} x^{1} & =\rho\sin\phi\cos\theta,\\ x^{2} & =\rho\sin\phi\sin\theta,\\ x^{3} & =\rho\cos\phi. \end{aligned} \right\} \end{equation}

The unit vectors along the three spherical coordinate directions are given by

(2.2)\begin{equation} \left. \begin{aligned} \boldsymbol{e}_\rho & =\sin\phi\cos\theta\,\boldsymbol{e}_1+\sin\phi\sin\theta\,\boldsymbol{e}_2+\cos\phi\,\boldsymbol{e}_3,\\ \boldsymbol{e}_\phi & =\cos\phi\cos\theta\,\boldsymbol{e}_1+\cos\phi\sin\theta\,\boldsymbol{e}_2-\sin\phi\,\boldsymbol{e}_3,\\ \boldsymbol{e}_\theta & ={-}\sin\theta\,\boldsymbol{e}_1+\cos\theta\,\boldsymbol{e}_2. \end{aligned} \right\} \end{equation}

For describing bulk radial motion of the bubble, we assume that $\boldsymbol {x}_e=r_e \boldsymbol {e}_\rho$ as the position vector of a material point in the FC with radius $r_e$. Similarly, let $\boldsymbol {x}=r(r_e,t)\boldsymbol {e}_\rho$ be the position vector of a material point in the DC. Then the bulk deformation gradient tensor, $\boldsymbol{\mathsf{A}}$, and the corresponding left Cauchy–Green deformation tensor, $\boldsymbol{\mathsf{B}}$, can be expressed as (Itskov Reference Itskov2019)

(2.3)\begin{equation} \left. \begin{aligned} \boldsymbol{\mathsf{A}} & ={\partial\boldsymbol{x}}/{\partial \boldsymbol{x}_e}=\frac{\partial{r}}{\partial r_e}\boldsymbol{e}_\rho\otimes\boldsymbol{e}_\rho+\frac{r}{r_e}\boldsymbol{e}_\phi\otimes\boldsymbol{e}_\phi+\frac{r}{r_e}\boldsymbol{e}_\theta\otimes\boldsymbol{e}_\theta,\\ \boldsymbol{\mathsf{B}} & =\boldsymbol{\mathsf{A}}\boldsymbol{\mathsf{A}}^\textrm{T}=\left(\frac{\partial{r}}{\partial r_e}\right)^2\boldsymbol{e}_\rho\otimes\boldsymbol{e}_\rho+\left(\frac{r}{r_e}\right)^2\boldsymbol{e}_\phi\otimes\boldsymbol{e}_\phi+\left(\frac{r}{r_e}\right)^2\boldsymbol{e}_\theta\otimes\boldsymbol{e}_\theta, \end{aligned} \right\} \end{equation}

where $\boldsymbol {c}\otimes \boldsymbol {d}$ represents the tensor product of two vectors $\boldsymbol {c}$ and $\boldsymbol {d}$.

Following the framework of Steigmann & Ogden (Reference Steigmann and Ogden1999), we idealize the gas-encapsulation and encapsulation-liquid interfaces as mathematical interfaces of zero thickness and call it Steigmann–Ogden interfaces (SOIs). Let $\theta ^1=\phi$ and $\theta ^2=\theta$ be the surface coordinates of the spherical interface of radius say $r_0$. In the convected coordinate system, let $r(r_0,t)$ be the deformed radius of the interface. For the bubble undergoing spherical deformation, the principal stretch ratio $\lambda =r(r_0,t)/r_0$. Let $\boldsymbol Z(\theta ^1,\theta ^2)$ and $\boldsymbol z(\theta ^1,\theta ^2)$ be the position vectors of the same point on the undeformed and deformed interfaces, respectively, given by

(2.4)\begin{equation} \left. \begin{aligned} \boldsymbol Z(\theta^1,\theta^2) & =r_0(\sin\theta^1\cos\theta^2\,\boldsymbol{e}_1+\sin\theta^1\sin\theta^2\,\boldsymbol{e}_2+\cos\theta^1\,\boldsymbol{e}_3),\\ \boldsymbol z(\theta^1,\theta^2) & =r(r_0,t)(\sin\theta^1\cos\theta^2\,\boldsymbol{e}_1+\sin\theta^1\sin\theta^2\,\boldsymbol{e}_2+\cos\theta^1\,\boldsymbol{e}_3). \end{aligned} \right\} \end{equation}

Due to the deformation of the bulk of the bubble, the interface can be assumed to be convected and the deformation mapping, $\boldsymbol \chi$, relates the same point before and after deformation as

(2.5)\begin{equation} \boldsymbol z(\theta^1,\theta^2)=\boldsymbol \chi(\boldsymbol Z(\theta^1,\theta^2)). \end{equation}

For the given position vectors of the interface, the tangent vectors $\boldsymbol G_\alpha$ and $\boldsymbol g_\alpha$ on the undeformed and deformed interfaces, respectively, are calculated using

(2.6a–d)\begin{equation} \boldsymbol{G}_\alpha=\boldsymbol{Z}_{,\alpha},\quad \boldsymbol{g}_\alpha=\boldsymbol{z}_{,\alpha},\quad \alpha\in\{1,2\},\quad ({\cdot})_{,\alpha}=\frac{\partial({\cdot})}{\partial\theta^\alpha}, \end{equation}

and the expressions for these vectors are given by

(2.7a,b)\begin{gather} \boldsymbol{G}_1=r_0\,\boldsymbol{e}_\phi,\quad \boldsymbol{G}_2=r_0\sin\phi\,\boldsymbol{e}_\theta, \end{gather}

(2.8a,b)\begin{gather}\boldsymbol{g}_1=r(r_0,t)\,\boldsymbol{e}_\phi,\quad \boldsymbol{g}_2=r(r_0,t)\sin\phi\,\boldsymbol{e}_\theta. \end{gather}

The components of the interface metric tensor for the undeformed and deformed interfaces, respectively, are given by

(2.9)\begin{gather} G_{\alpha\beta}=\begin{pmatrix} \boldsymbol{G}_1\boldsymbol{\cdot}\boldsymbol{G}_1 & \boldsymbol{G}_1\boldsymbol{\cdot}\boldsymbol{G}_2\\ \boldsymbol{G}_2\boldsymbol{\cdot}\boldsymbol{G}_1 & \boldsymbol{G}_2\boldsymbol{\cdot}\boldsymbol{G}_2 \end{pmatrix}=r_0^2\;{\textrm{{diag}}}(1,\sin^2\phi), \end{gather}

(2.10)\begin{gather}g_{\alpha\beta}=\begin{pmatrix} \boldsymbol{g}_1\boldsymbol{\cdot}\boldsymbol{g}_1 & \boldsymbol{g}_1\boldsymbol{\cdot}\boldsymbol{g}_2\\ \boldsymbol{g}_2\boldsymbol{\cdot}\boldsymbol{g}_1 & \boldsymbol{g}_2\boldsymbol{\cdot}\boldsymbol{g}_2 \end{pmatrix}=r(r_0,t)^2\;{\textrm{{diag}}}(1,\sin^2\phi) =\lambda^2 \it{G}_{\alpha\beta}. \end{gather}

For the given undeformed metric tensor components $G_{\alpha \beta }$, the Christoffel symbols are defined as

(2.11)\begin{equation} \tilde{\varGamma}_{\alpha\beta}^\gamma=\tfrac{1}{2}G^{\gamma\delta}(G_{\alpha\delta,\beta}+G_{\delta\beta,\alpha}-G_{\beta\alpha,\delta}). \end{equation}

Similarly, $\varGamma _{\alpha \beta }^\gamma$ are the Christoffel symbols constructed using the deformed interface metric tensor components $g_{\alpha \beta }$. The covariant derivative of a tensor generalizes the concept of partial derivative onto the curved surface. The following relation obtains the calculation of the components of the covariant derivative $(:)$ of a vector $(v^\alpha )$:

(2.12)\begin{equation} {v^\alpha}_{:\beta}={v^\alpha}_{,\beta}+\varGamma^\alpha_{\gamma\beta}\,v^\gamma. \end{equation}

The dual tangent vectors represented by $\boldsymbol {G}^\alpha$ and $\boldsymbol {g}^\alpha$ associated with the coordinates $\theta ^\alpha$ are given by

(2.13a,b)\begin{equation} \boldsymbol{G}^\alpha=G^{\alpha\beta}\boldsymbol{G}_\beta,\quad \boldsymbol{g}^\alpha=g^{\alpha\beta}\boldsymbol{g}_\beta. \end{equation}

The calculated values of the components of the dual metric tensors $G^{\alpha \beta }$ and $g^{\alpha \beta }$, respectively, are obtained as

(2.14)\begin{gather} G^{\alpha\beta}=r_0^2\;\textrm{{diag}}(1,\csc^2\phi), \end{gather}

(2.15)\begin{gather}g^{\alpha\beta}=r(r_0,t)^2\;\textrm{{diag}}(1,\csc^2\phi)= \lambda^{-2}\, \it{G}^{\alpha\beta}. \end{gather}

The second fundamental forms $Q_{\alpha \beta }$ and $q_{\alpha \beta }$ representing the normal curvatures associated with the deformed and the undeformed surfaces, respectively, are given by

(2.16a,b)\begin{equation} Q_{\alpha\beta}=\boldsymbol{N}\boldsymbol{\cdot}\boldsymbol{G}_{\alpha:\beta},\quad q_{\alpha\beta}=\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{g}_{\alpha:\beta}. \end{equation}

Here, $\boldsymbol {N}$ and $\boldsymbol {n}$ are the outward normal vectors to the undeformed and deformed interfaces, respectively, given by

(2.17)\begin{equation} \left. \begin{aligned} \boldsymbol{N} & =\tfrac{1}{2}\mu^{\alpha\beta}\boldsymbol{G}_\alpha\times\boldsymbol{G}_\beta,\\ \boldsymbol{n} & =\tfrac{1}{2}\varepsilon^{\alpha\beta}\boldsymbol{g}_\alpha\times\boldsymbol{g}_\beta,\quad\text{(summation convention)} \end{aligned} \right\} \end{equation}

where $\varepsilon ^{\alpha \beta }=e^{\alpha \beta }/\sqrt {g}, \mu ^{\alpha \beta }=e^{\alpha \beta }/\sqrt {G}$ with $G=\textrm {det}(G_{\alpha \beta })=r_0^4\sin ^2\phi, g=\textrm {det}(g_{\alpha \beta })=r(r_0,t)^4\sin ^2\phi$ and $e^{\alpha \beta }=e_{\alpha \beta }$ is the alternator symbol ($e^{12}=-e^{21}=1, e^{11}=e^{22}=0$). Furthermore, we calculate the relations

(2.18a–c)\begin{equation} \boldsymbol{G}_{1,1}={-}r_0\,\boldsymbol{e}_r,\quad \boldsymbol{G}_{1,2}=\boldsymbol{G}_{2,1}=r_0\cos\phi\,\boldsymbol{e}_\theta, \quad \boldsymbol{G}_{2,2}={-}r_0\sin\phi\,(\cos\theta\,\boldsymbol{e}_1+\sin\theta\,\boldsymbol{e}_2), \end{equation}

and obtain the expression for $Q_{\alpha \beta }$ using the relations in (2.16a,b) and (2.17) where $\boldsymbol {N}=\boldsymbol {e}_\phi \times \boldsymbol {e}_\theta =\boldsymbol {e}_\rho$ as

(2.19)\begin{equation} Q_{\alpha\beta}={-}r_0\;\textrm{{diag}}(1,\sin^2\phi) ={-}{{r}}_0^{{-}1}{G}_{\alpha\beta}. \end{equation}

From the kinematics of the interface deformation, the deformation gradient tensor $\boldsymbol{\mathsf{a}}$ can be written as (Itskov Reference Itskov2019)

(2.20)\begin{equation} \boldsymbol{\mathsf{a}}=\boldsymbol{g}_\beta\otimes\boldsymbol{G}^\beta. \end{equation}

The Jacobian $(J)$ for the deformation is defined as the ratio of deformed to underformed interface areas, given by

(2.21)\begin{equation} J=\textrm{det}\,\boldsymbol{\mathsf{a}}=\sqrt{g/G}. \end{equation}

The right Cauchy-Green deformation tensor ${\boldsymbol{\mathsf{C}}}$ is isotropic and can be expressed as

(2.22)\begin{equation} \boldsymbol{\mathsf{C}}=\boldsymbol{\mathsf{a}}^\textrm{T}\boldsymbol{\mathsf{a}}=g_{\alpha\beta}\boldsymbol{G}^\alpha\otimes\boldsymbol{G}^\beta. \end{equation}

The curvature tensor $(\boldsymbol{\mathsf{q}})$ and the relative curvature tensor $(\boldsymbol {\kappa })$, related to the deformed and the undeformed configurations, respectively, are defined as

(2.23)\begin{equation} \left. \begin{aligned} \boldsymbol{\mathsf{q}} & =q_{\alpha\beta}\,\boldsymbol{g}^\alpha\otimes\boldsymbol{g}^\beta,\\ \boldsymbol\kappa & =\kappa_{\alpha\beta}\,\boldsymbol{G}^\alpha\otimes\boldsymbol{G}^\beta,\\ \end{aligned} \right\} \end{equation}

with the relation $\boldsymbol \kappa =-\boldsymbol{\mathsf{a}}^\textrm {T}\boldsymbol{\mathsf{q}}\,\boldsymbol{\mathsf{a}}$. Then the relation between the components of the above two curvature tensors is $\kappa _{\alpha \beta }=-q_{\alpha \beta }$, where $q_{\alpha \beta }$ is calculated using (2.16a,b).

In a detailed study accounting for the material symmetry of the interface, Steigmann (Reference Steigmann2001) concluded that the interfacial energy is not an isotropic scalar-valued tensor function corresponding to the reference configuration, rather it is a function of the right Cauchy–Green interface deformation tensor $\boldsymbol{\mathsf{C}}$ and the relative curvature tensor $\boldsymbol {\kappa }$. Here, the interface is assumed to be a micropolar, which is not isotropic with respect to inversion, also known as hemitropic. In other words, an interface whose symmetry group consists of all rotations but no reflections are said to be hemitropic. For such an hemitropic interface, its energy density is expressed as (Fung Reference Fung1977; Steigmann & Ogden Reference Steigmann and Ogden1999; Gao et al. Reference Gao, Huang, Qu and Fang2014)

(2.24)\begin{equation} \gamma=\gamma(\boldsymbol{\mathsf{C}},{\boldsymbol \kappa}), \end{equation}

and satisfies the relation

(2.25)\begin{equation} \gamma(\boldsymbol{\mathsf{C}},{\boldsymbol \kappa})=\gamma(\boldsymbol{\mathsf{Q}}\boldsymbol{\mathsf{C}}\boldsymbol{\mathsf{Q}}^\textrm{T},\boldsymbol{\mathsf{Q}}{\boldsymbol \kappa}\boldsymbol{\mathsf{Q}}^\textrm{T}), \end{equation}

where $\boldsymbol{\mathsf{Q}}$ is a proper-orthogonal second-order tensor.

Hence, according to standard representation theory, the interface energy density can be expressed in terms of six basis invariants of the right Cauchy–Green interface deformation tensor $\boldsymbol{\mathsf{C}}$, the relative curvature tensor $\boldsymbol {\kappa }$ and the permutation tensor density $\boldsymbol {\mu }$ on the undeformed interface as

(2.26)\begin{equation} \gamma=\gamma(I_1, I_2, I_3, I_4, I_5, I_6), \end{equation}

and the invariants are given by (see Appendix A)

(2.27)\begin{equation} \left. \begin{aligned} I_1 & =\textrm{tr}\,\boldsymbol{\mathsf{C}}=G^{\alpha\beta}C_{\alpha\beta},\\ I_2 & =\textrm{det}\,\boldsymbol{\mathsf{C}}=J^2=g/G,\\ I_3 & =\textrm{tr}{\,\boldsymbol \kappa}=G^{\alpha\beta}\kappa_{\alpha\beta},\\ I_4 & =\textrm{det}{\,\boldsymbol \kappa}=\tfrac{1}{2}\mu^{\alpha\beta}\mu^{\gamma\delta}\kappa_{\alpha\gamma}\kappa_{\beta\delta},\\ I_5 & =\textrm{tr}(\boldsymbol{\mathsf{C}}\boldsymbol{\kappa})=C_{\alpha\beta}\kappa^{\alpha\beta}=C^{\alpha\beta}\kappa_{\alpha\beta},\\ I_6 & =\textrm{tr}(\boldsymbol{\mathsf{C}}\boldsymbol{\kappa}\boldsymbol{\mu})=G_{\alpha\beta}C_{\gamma\delta}\kappa^{\alpha\gamma}\mu^{\beta\delta}=G_{\alpha\beta}\kappa_{\gamma\delta}C^{\alpha\gamma}\mu^{\beta\delta}. \end{aligned} \right\} \end{equation}

The choice of these basis invariants are unique to hemitropic material property. These invariants mentioned above were first listed by Zheng (Reference Zheng1993) and further, the two-dimensional Cayley–Hamilton theorem can be used to obtain the equivalence of these invariants to that provided by Zheng. The expressions for the Cauchy interface stress $(\boldsymbol \sigma )$ and the moment $({\boldsymbol{\mathsf{m}}})$ tensors are calculated using the relations (for detailed derivation see Steigmann & Ogden Reference Steigmann and Ogden1999; Gao et al. Reference Gao, Huang, Qu and Fang2014)

(2.28a,b)\begin{equation} {\boldsymbol \sigma}=2\,\frac{\partial \gamma}{\partial{\boldsymbol{\mathsf{C}}}},\quad {{\boldsymbol{\mathsf{m}}}}=\frac{\partial\gamma}{\partial \boldsymbol \kappa}. \end{equation}

Using the expression for interface energy in terms of invariants, one can obtain the expressions for the components of interface stress tensor and the Eulerian bending moment tensor as

(2.29)\begin{equation} \left. \begin{aligned} \frac{1}{2}J\sigma^{\alpha\beta} & =\frac{\partial \gamma}{\partial I_1}G^{\alpha\beta}+ \frac{\partial \gamma}{\partial I_2}\tilde{C}^{\alpha\beta}+\frac{\partial \gamma}{\partial I_5}\kappa^{\alpha\beta}+\frac{1}{2}\frac{\partial \gamma}{\partial I_6}\left(D^{\alpha\beta}+D^{\beta\alpha}\right),\\ Jm^{\alpha\beta} & =\frac{\partial \gamma}{\partial I_3}G^{\alpha\beta}+ \frac{\partial \gamma}{\partial I_4}\tilde{\kappa}^{\alpha\beta}+\frac{\partial \gamma}{\partial I_5}C^{\alpha\beta}+\frac{1}{2}\frac{\partial \gamma}{\partial I_6}\left(E^{\alpha\beta}+E^{\beta\alpha}\right), \end{aligned} \right\} \end{equation}

where

(2.30a,b)\begin{equation} D^{\alpha\beta}=G_{\gamma\delta}\mu^{\alpha\gamma}\kappa^{\beta\delta},\quad E^{\alpha\beta}=G_{\gamma\delta}\mu^{\alpha\gamma}C^{\beta\delta}. \end{equation}

Here, $\tilde {\boldsymbol{\mathsf{C}}}$ and $\boldsymbol {\tilde {\kappa }}$ are the adjugate of $\boldsymbol{\mathsf{C}}$ and $\boldsymbol {\kappa }$, respectively, with components given by

(2.31)\begin{gather} \tilde{C}^{\alpha\beta}=(\textrm{tr}\,\boldsymbol{\mathsf{C}})G^{\alpha\beta}-C^{\alpha\beta}, \end{gather}

(2.32)\begin{gather}\tilde{\kappa}^{\alpha\beta}=(\textrm{tr}\,\boldsymbol{\kappa})G^{\alpha\beta}-\kappa^{\alpha\beta}. \end{gather}

We assume that the bubble undergoes deformation from a FC such that the bulk deformation tensor $\boldsymbol{\mathsf{A}}$ is compatible with the interface deformation tensor $\boldsymbol{\mathsf{a}}$. From the balance of linear momentum, the general governing equation for the radial dynamics of the shell material and the outer liquid can be expressed in terms of the bulk stress field ${\boldsymbol{\mathsf{T}}}$ and velocity field $\boldsymbol {u}(r,t)$ as

(2.33)\begin{equation} \varrho\left[{\boldsymbol u}_{,t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}){\boldsymbol u}\right]=\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}, \end{equation}

where $\varrho$ is the density of the medium. Due to the existence of interface energy, the generalized Young–Laplace interface jump condition can be expressed in terms of interface stress and bending moment as (for detailed derivations, see Steigmann & Ogden Reference Steigmann and Ogden1999; Gao et al. Reference Gao, Huang, Qu and Fang2014)

(2.34)\begin{equation} [\![{\boldsymbol{\mathsf{T}}}]\!]\boldsymbol{\cdot n}=\left[(\sigma^{\alpha\beta}-m^{\lambda\alpha}q^{\beta}_{\lambda})\boldsymbol{g}_\beta+m^{\beta\alpha}_{\quad:\beta}\boldsymbol{n}\right]_{:\alpha}, \end{equation}

where the surface covariant derivative of a second-order tensor $m^{\beta \alpha }_{:\beta }$ is given by the expression

(2.35)\begin{equation} m^{\beta\alpha}_{:\beta}=m^{\beta\alpha}_{,\beta}+m^{\delta\alpha}\,\varGamma^\beta_{\;\,\delta\beta}+m^{\beta\delta}\,\varGamma^\alpha_{\delta\beta}, \end{equation}

and $[\![{\boldsymbol{\mathsf{T}}}]\!]$ is the jump in the stress due to interface energy. Similar to the calculations of $Q_{\alpha \beta }$ in (2.19), using the relations (2.10), (2.15) and (2.17), the expression for the curvature tensor are obtained for the inner (minus sign) and outer (plus sign) interfaces as

(2.36)\begin{equation} \left. \begin{aligned} q_{\alpha\beta} & ={-}({\mp} r(r_0,t)^{{-}1}g_{\alpha\beta})={-}({\mp} r_0^{{-}1}\lambda G_{\alpha\beta}),\\ q^\alpha_{\beta} & =g^{\alpha\gamma}q_{\beta\gamma}={-}({\mp} r_0^{{-}1}\lambda^{{-}1} \delta^\alpha_{\beta}), \end{aligned} \right\} \end{equation}

and the expressions for the components of the relative curvature tensor $\boldsymbol {\kappa }$ are given by

(2.37)\begin{equation} \left. \begin{aligned} \kappa_{\alpha\beta} & ={\mp} r_0^{{-}1}\lambda G_{\alpha\beta},\\ \kappa^{\alpha\beta} & ={\mp} r_0^{{-}1}\lambda^3 g^{\alpha\beta}. \end{aligned} \right\} \end{equation}

Therefore, the contra- and co-variant components of the interface deformation tensor components can be written as

(2.38)\begin{equation} \left. \begin{aligned} C_{\alpha\beta} & =\lambda^2 G_{\alpha\beta},\\ C^{\alpha\beta} & =G^{\alpha\gamma}G^{\beta\delta}C_{\gamma\delta}=\lambda^4g^{\alpha\beta}. \end{aligned} \right\} \end{equation}

The invariants in (2.27) are calculated as

(2.39a–f)\begin{equation} I_1=2\,\lambda^2, \quad I_2=\lambda^4, \quad I_3=\mp2\,\lambda/r_0, \quad I_4=\lambda/r_0^2, \quad I_5=\mp2\,\lambda^3/r_0, \quad I_6=0. \end{equation}

For the undeformed interface configuration, the identity tensor (unit tensor) $\boldsymbol{\mathsf{1}}$ is given by

(2.40)\begin{equation} \boldsymbol{\mathsf{1}}=G_{\alpha\beta}\,\boldsymbol{G}^\alpha\otimes\boldsymbol{G}^\beta. \end{equation}

Using the identity tensor one can write $\boldsymbol{\mathsf{C}}=\lambda ^2\boldsymbol{\mathsf{1}}$ and the components of the adjugate of $\boldsymbol{\mathsf{C}}$ are $\tilde {C}^{\alpha \beta }=C^{\alpha \beta }.$ From the relations in (2.29), components of interface stress and Eulerian bending moment tensors, respectively, can be simplified to yield the relations

(2.41a,b)\begin{equation} \sigma^{\alpha\beta}=\sigma g^{\alpha\beta},\quad m^{\alpha\beta}=m g^{\alpha\beta}, \end{equation}

where

(2.42)\begin{gather} \sigma=\gamma_0+2\left(\gamma_1+\gamma_2\lambda^2\mp\frac{\gamma_5\lambda}{r_0}\right), \end{gather}

(2.43)\begin{gather}m=\gamma_3\mp\frac{\gamma_4\lambda}{r_0}+\gamma_5\lambda^2, \end{gather}

with ‘$-$’ and ‘$+$’ signs are supposed to be used for the inner and outer interfaces, respectively. Here, deformation independent interface tension is introduced as $\gamma _0$, and $\gamma _i=\partial \gamma /\partial I_i, i=1,2,$ are considered as interface material parameters. The generalized Young–Laplace equation (2.34) for the radial dynamics can be expressed as

(2.44)\begin{equation} [\![{\boldsymbol{\mathsf{T}}}]\!]\boldsymbol{\cdot} {\boldsymbol n}={-}({\mp} 2 r_0^{{-}1}\lambda^{{-}1}\varSigma)\;\boldsymbol{n}, \end{equation}

where $\varSigma =\sigma \mp r_0^{-1}\lambda ^{-1}\,\textrm {m}$. The residual interface stress and bending moment tensors at static equilibrium are calculated by the substitution $\lambda =1$, which yields $\sigma ^{*\alpha \beta }=\sigma ^*g^{\alpha \beta }$ and $m^{*\alpha \beta }=m^*g^{\alpha \beta }$, where $\sigma ^*$ and $m^*$ are given by

(2.45)\begin{gather} \sigma^*=\gamma_0+2\left(\gamma_1+\gamma_2\mp\frac{\gamma_5}{r_0}\right), \end{gather}

(2.46)\begin{gather}m^*=\gamma_3+\gamma_5\mp\frac{\gamma_4}{r_0}. \end{gather}

Here, the assumption of $\gamma _i=\partial \gamma /\partial I_i$ as constant gives a simple linear structure to the interface energy density function with a linear combination of the basis invariants defined in (2.27). However, in a more general case, the interface energy can be a nonlinear function of these basis invariants. Therefore, there is always scope to improve and modify the mathematical structure of the Young–Laplace equation and the residual stress field.

Comparing the residual stress $\sigma ^*$ at the interface with two leading-order terms of well-known Tolman's formula (Tolman Reference Tolman1949) $\sigma _s\sim \sigma _\infty (1-2\delta _\infty /r_0)$ for a curved interface indicates that $\gamma _5$ has the role of Tolman's length scale $\delta _\infty$. In $\varSigma$, $\gamma _3$ appears to have a similar interpretation to $\gamma _5$ for $\lambda =1$; however, it is important to identify the distinct effects of $\gamma _3$ and $\gamma _5$ on the interface tension, as will be shown in (2.53).

As mentioned earlier, the present model naturally introduces residual stress in the bulk due to the expansion or shrinkage of the bubble after preparation, which was also observed in the experiments (de Jong et al. Reference de Jong, Emmer, Chin, Bouakaz, Mastik, Lohse and Versluis2007). This also indicates the important feature of equivalent NC bubbles with the same radii and internal gas pressure at the NC, which can show completely different mechanical behaviour. This behaviour is strongly connected to the process of attaining NC, which in-turn is related to the interface energy induced residual stress. On the other hand, one would also expect that different size bubbles made of the same shell material should demonstrate the behaviour based on same material constants. Here, we demonstrate these features by considering the above interface energy and corresponding interface parameters in the mathematical model of EB dynamics.

2.2. Constitutive relation for the shell and the fluid

We analyse the EB dynamics by assuming that the shell material and the fluid surrounding the bubble are homogeneous, isotropic and incompressible materials. Furthermore, the shell is assumed to be viscoelastic. A Kelvin–Voigt type constitutive model with elastic part described by the Mooney–Rivlin (MR) material model with elastic constants $C_1$ and $C_2$, and a viscous part by the Newton's law of viscosity is considered. In order to capture the material nonlinearity of the shell material, we have chosen this basic nonlinear material model. The ratio $C_2/C_1$ in the MR material model is usually interpreted as a strain-softening or hardening parameter depending on its value. For an incompressible bulk MR material, the strain energy density $W=C_1(\tilde {I}_1-3)+C_2(\tilde {I}_2-3)$, where $\tilde {I}_1,\;\tilde {I}_2,\;\tilde {I}_3$ are the three invariants of ${\boldsymbol{\mathsf{B}}}$ given by

(2.47)\begin{equation} \left. \begin{aligned} \tilde{I}_1 & =\textrm{tr}{\,{\boldsymbol{\mathsf{B}}}}=2\varLambda^2+\varLambda^{{-}4},\\ \tilde{I}_2 & =\textrm{tr}{\,{\boldsymbol{\mathsf{B}}}^2}=2\varLambda^{{-}2}+\varLambda^4,\\ \tilde{I}_3 & =\textrm{det}\,{\boldsymbol{\mathsf{B}}}=1, \end{aligned} \right\} \end{equation}

where $\varLambda =r/r_e$ is the principal stretch ratio in both $\theta ^1$ and $\theta ^2$ directions. Since the bulk is assumed to be incompressible, ${\tilde {I}}_3=1$ gives

(2.48)\begin{equation} \frac{\partial{r}}{\partial r_e}=\left(\frac{r_e}{r}\right)^2=\varLambda^{{-}2}. \end{equation}

Using superscripts $S$ and $L$ for the shell and liquid parameters, respectively, the Cauchy stress tensor for the shell is given by

(2.49) \begin{equation} {\boldsymbol{\mathsf{T}}}^{S}={-}p^{S}{\boldsymbol{\mathsf{I}}}+2(C_1+\tilde{I}_1C_2){\boldsymbol{\mathsf{B}}}-2C_2{\boldsymbol{\mathsf{B}}}^2+\eta^{S}[\boldsymbol{\nabla}\boldsymbol{u}^{S}+(\boldsymbol{\nabla}\boldsymbol{u}^{S})^\textrm{T}], \end{equation}

where $\eta ^{S}$ is the shell viscosity coefficient, $\boldsymbol{\mathsf{I}}$ is the identity tensor and $p^{S}$ is commonly referred to as the arbitrary hydrostatic pressure due to incompressibility. The liquid viscous stress tensor with viscosity coefficient $\eta ^{L}$ is given by the constitutive equation

(2.50)\begin{equation} {\boldsymbol{\mathsf{T}}}^{L}={-}p^{L}{{\boldsymbol{\mathsf{I}}}}+\eta^{L}[\boldsymbol{\nabla}\boldsymbol{u}^{L}+(\boldsymbol{\nabla}\boldsymbol{u}^{L})^\textrm{T}]. \end{equation}

2.3. Governing equation

For the radial dynamics of the shell exerted by an ultrasound field and the surrounding liquid, the governing equation in spherical polar coordinates with $(r,\phi,\theta )$ as principal coordinate directions can be reduced to the radial equation

(2.51)\begin{equation} \varrho\left(\frac{\partial u}{\partial t}+{u}\frac{\partial u}{\partial r}\right)=\frac{\partial T_{rr}}{\partial r}+\frac{2 T_{rr}-T_{\phi\phi}-T_{\theta\theta}}{r}. \end{equation}

The stress jump conditions at the inner and outer interfaces are

(2.52)\begin{equation} \left. \begin{aligned} \left. T_{rr}^{S}\right\lvert_{r=R_1}+p_{g} & =\frac{2\sigma_1}{R_1},\\ \left[T_{rr}^{S}-T_{rr}^{L}\right]_{r=R_2} & ={-}\frac{2\sigma_2}{R_2}, \end{aligned} \right\} \end{equation}

where $\sigma _1$ and $\sigma _2$ are the interface tension parameters given by

(2.53)\begin{equation} \left. \begin{aligned} \sigma_1 & =\gamma_{10}+2\left(\gamma_{11}+\gamma_{12}\frac{R_1^2}{R_{10}^2}-\gamma_{15}\frac{R_1}{R_{10}^2}\right)-\frac{1}{R_1}\left(\gamma_{13}-\gamma_{14}\frac{R_1}{R_{10}^2}+\gamma_{15}\frac{R_1^2}{R_{10}^2}\right),\\ \sigma_2 & =\gamma_{20}+2\left(\gamma_{21}+\gamma_{22}\frac{R_2^2}{R_{20}^2}+\gamma_{25}\frac{R_2}{R_{20}^2}\right)+\frac{1}{R_2}\left(\gamma_{23}+\gamma_{24}\frac{R_2}{R_{20}^2}+\gamma_{25}\frac{R_2^2}{R_{20}^2}\right), \end{aligned} \right\} \end{equation}

with $\gamma _{i(\cdot )}, i=1,2$ correspond to the inner and outer interfaces, respectively. Unlike existing shell models, interface tensions are considered on either side of the bubble shell and dependent on the interface radius.

Substituting $R_1=R_{10}$ and $R_2=R_{20}$ in the stress jump conditions in (2.52), one can rewrite the interface conditions at the NC of the bubble as

(2.54)\begin{gather} \left. T_{rr}^{S}\right\lvert_{R_1=R_{10}}+p_{g_0}=\frac{2}{R_{10}}\left[\gamma_{10}+2\left(\gamma_{11}+\gamma_{12}-\frac{\gamma_{15}}{R_{10}}\right)-\frac{1}{R_{10}}\left(\gamma_{13}-\frac{\gamma_{14}}{R_{10}}+\gamma_{15}\right)\right], \end{gather}

(2.55)\begin{gather}\left. T_{rr}^{S}\right\lvert_{R_2= R_{20}}\left.-T_{rr}^{L}\right\lvert_{R_2=R_{20}}\,{=}\,\frac{-2}{R_{20}}\left[\gamma_{20}\,{+}\,2\left(\gamma_{21}\,{+}\,\gamma_{22}\,{+}\,\frac{\gamma_{25}}{R_{20}}\right)\,{+}\,\frac{1}{R_{20}}\left(\gamma_{23}+\frac{\gamma_{24}}{R_{20}}\,{+}\,\gamma_{25}\right)\right]. \end{gather}

To obtain the governing equation, integrating (2.51) with respect to the radial coordinate $r$ for the shell and liquid regions separately, we have

(2.56)$$\begin{gather} \int_{R_1}^{R_2} \varrho^{S}\left(\frac{\partial u^{S}}{\partial t}+{u^{S}}\frac{\partial u^{S}}{\partial r}\right) \,\textrm{d}r+\int_{R_2}^{\infty} \varrho^{L}\left(\frac{\partial u^{L}}{\partial t}+{u^{L}}\frac{\partial u^{L}}{\partial r}\right) \,\textrm{d}r \nonumber\\ = \int_{R_1}^{R_2} \left(\frac{\partial T^{S}_{rr}}{\partial r}+\frac{2T_{rr}^{S}-T^{S}_{\phi\phi}-T^{S}_{\theta\theta}}{r}\right) \,\textrm{d}r+\int_{R_2}^\infty \left(\frac{\partial T^{L}_{rr}}{\partial r}+\frac{2T^{L}_{rr}-T^{L}_{\phi\phi}-T^{L}_{\theta\theta}}{r}\right) \,\textrm{d}r, \end{gather}$$

where the superscripts $S$ and $L$ represent the physical terms corresponding to the shell and liquid, respectively. From (2.49) we calculate the expressions for the components of the Cauchy stress tensor for the shell as

(2.57)\begin{equation} \left. \begin{aligned} T^{S}_{rr} & ={-}p^{S}+2\lambda^{{-}4}C_1+4\lambda^{{-}2}C_2+2\eta^{S}\left(\frac{\partial u^{S}}{\partial r}\right),\\ T^{S}_{\phi\phi} & =T^{S}_{\theta\theta}={-}p^{S}+2\lambda^{2}C_1+2(\lambda^4+\lambda^{{-}2})C_2+2\eta^{S}\left(\frac{u^{S}}{r}\right),\\ T^{S}_{r\phi} & =T^{S}_{\phi\theta}=T^{S}_{\theta r}=0. \end{aligned} \right\} \end{equation}

Similarly from (2.50), we calculate the components of the liquid viscous stress tensor as

(2.58)\begin{equation} \left. \begin{aligned} T^{L}_{rr} & ={-}p^{L}+2\eta^{L}\left(\frac{\partial u^{L}}{\partial r}\right),\\ T^{L}_{\phi\phi} & =T^{L}_{\theta\theta}={-}p^{L}+2\eta^{L}\left(\frac{u^{L}}{r}\right),\\ T^{L}_{r\phi} & =T^{L}_{\phi\theta}=T^{L}_{\theta r}=0. \end{aligned} \right\} \end{equation}

Using the above expressions for the stress components normal to the three principal coordinate planes and the expression for the velocity field in Appendix B, the governing equation for the radial dynamics of a bubble exerted by an ultrasound field is given by

(2.59)$$\begin{gather} \varrho^{S}\left[R_1\ddot{R}_1+\frac{3}{2}\dot{R}_1^2\right]-\left(\varrho^{S}-\varrho^{L}\right)\left[R_2\ddot{R}_2+\frac{3}{2}\dot{R}_2^2\right]=p_g-p_0-p_d+p_e-\frac{2\sigma_1}{R_1}-\frac{2\sigma_2}{R_2}\nonumber\\ -4\eta^{L}\frac{\dot{R}_2}{R_2}-4\eta^{S}\left(\frac{\dot{R}_1}{R_1}-\frac{\dot{R}_2}{R_2}\right), \end{gather}$$

where $p_g$ is the gas pressure, $p_0$ is the ambient pressure, $p_d=p_a\sin (2{\rm \pi} f t)$ is the external acoustic field with pressure $p_a$ and frequency $f$, and $p_e$ is the elastic restoring force (see Appendix C) given by

(2.60) \begin{align} p_e&={-}4\left[C_1\left(\frac{1}{4\varLambda_2^4}+\frac{1}{\varLambda_2}\right)+C_2\left(\frac{1}{2\varLambda_2^2}-\varLambda_2\right)\right]\nonumber\\&\quad +4\left[C_1\left(\frac{1}{4\varLambda_1^4}+\frac{1}{\varLambda_1}\right)+C_2\left(\frac{1}{2\varLambda_1^2}-\varLambda_1\right)\right]. \end{align}

Here, $\varLambda _1=R_1/R_{e_1}$ and $\varLambda _2=R_2/R_{e_2}$ are the stretch ratios. The interface tension parameters $(\sigma _1,\sigma _2)$ and the elastic restoring force $(p_e)$ are the additional terms which makes the governing equation (2.59) different from that of Church's model (Church Reference Church1995). When the bubble is in its NC the time-dependent terms in the governing equation vanish and the static equilibrium equation is given by

(2.61)\begin{align} p_{g_0}&=p_0+\frac{2}{R_{10}}\left[\gamma_{10}+2\left(\gamma_{11}+\gamma_{12}-\frac{\gamma_{15}}{R_{10}}\right)-\frac{1}{R_{10}}\left(\gamma_{13}-\frac{\gamma_{14}}{R_{10}}+\gamma_{15}\right)\right]\nonumber\\ &\quad +\frac{2}{R_{20}}\left[\gamma_{20}+2\left(\gamma_{21}+\gamma_{22}+\frac{\gamma_{25}}{R_{20}}\right)+\frac{1}{R_{20}}\left(\gamma_{23}+\frac{\gamma_{24}}{R_{20}}+\gamma_{25}\right)\right]\nonumber\\ &\quad +4\left[C_1\left(\frac{1}{4\varLambda_{20}^4}+\frac{1}{\varLambda_{20}}\right)+C_2\left(\frac{1}{2\varLambda_{20}^2}-\varLambda_{20}\right)\right]\nonumber\\ &\quad -4\left[C_1\left(\frac{1}{4\varLambda_{10}^4}+\frac{1}{\varLambda_{10}}\right)+C_2\left(\frac{1}{2\varLambda_{10}^2}-\varLambda_{10}\right)\right]. \end{align}

Let $p_{g_e}$ be the uniform gas pressure inside the bubble in the FC. Assuming that the gas inside the bubble is undergoing expansion or compression with polytropic expansion index $k=1.4$, the relation between $p_g,p_{g_0}$ and $p_{g_e}$ can be expressed as

(2.62)\begin{equation} p_g=p_{g_0}\left(\frac{R_{10}}{R_1}\right)^{3k}=p_{g_e}\left(\frac{R_{e_1}}{R_1}\right)^{3k}. \end{equation}

With the inner and outer radii of the bubble at NC being known, the bubble radii at FC can be found by solving the nonlinear algebraic equation (2.61).

As we are interested in bubbles of radius ${O}(10^{-6})$ m and thickness ${O}(10^{-9})$ m, inspecting the terms in (2.53), it is reasonable to assume the order of parameters $\gamma _{10},\gamma _{20},\gamma _{11},\gamma _{21},\gamma _{12},\gamma _{22}$ as ${O}(1)\,\textrm {N}\,\textrm {m}^{-1}$, $\gamma _{13},\gamma _{23},\gamma _{15},\gamma _{25}$ as ${O}(10^{-6}$ )N and $\gamma _{14},\gamma _{24}$ as ${O}(10^{-12})$ N m. Primarily this assumption ensures that the effective interface tension parameters $\sigma _1$ and $\sigma _2$ possess reasonable values of ${O}(1)\,\textrm {N}\,\textrm {m}^{-1}$. Also, it is worth mentioning that the order of these interface parameters are chosen such that it will not lead to large deviations in the bubble radii while undergoing deformation from the FC to the NC. In the rest of the paper, $\varrho ^{S}=1100\,\textrm {kg}\,\textrm {m}^{-3}$, $\varrho ^{L}=1000\,\textrm {kg}\,\textrm {m}^{-3}$, $\gamma _{10}=0.04\,\textrm {N}\,\textrm {m}^{-1}$, $\gamma _{20}=0.005\,\textrm {N}\,\textrm {m}^{-1}$ (Krasovitski & Kimmel Reference Krasovitski and Kimmel2006; Shao & Chen Reference Shao and Chen2015) and a numerical value of $\gamma _{ij}$ by default carries these orders apart from any other explicitly mentioned factors.

In (2.53), $\gamma _{i2}$ term retains the effect of relative interface area. However, in the present model, the dynamic interface tension cannot become zero due to the change in the interface area, unlike the other models (Paul et al. Reference Paul, Katiyar, Sarkar, Chatterjee, Shi and Forsberg2010). On the other hand, both $\gamma _{13}$ and $\gamma _{15}$ can make the dynamic interface tension negative at the inner interface, which is purely an effect of interface bending rigidity. These relations also show distinct effects of $\gamma _3$ and $\gamma _5$ terms when the interface starts deforming. It is interesting to observe that the contribution of $\gamma _5$ terms in (2.59) is only through static equilibrium pressure and depends only on the initial size of the bubble. On the other hand, the dynamic interface tension has terms proportional to the interface radius and curvature through $\gamma _5$ and $\gamma _3$, respectively, which shows distinct effects.

Here, we study the influence of three important interface parameters $\gamma _2,\gamma _3$ and $\gamma _5$ on the EB's NC and radial dynamics. These three interface parameters are related to the deformation dependent interface area, curvature and the coupling between the two. To understand the influence of these aspects further, we have constructed a multistart optimization problem by defining a compression or expansion ratio to estimate the interface and material parameters.