2.1. Magneto-elastic equations

The bubble surface is encapsulated by a lipid/polymer coating infused with magnetic nanoparticles. This facilitates the translation, radial and non-spherical oscillations of the bubble on the application of a magnetic field. In general, the strain energy density $W$ of magneto-elastomer is a function of $\boldsymbol{F}$ and $\boldsymbol m$ (Steigmann Reference Steigmann2004). Here, $\boldsymbol{F}$ is the gradient of the deformed position $\boldsymbol x$ with respect to reference position $\boldsymbol X$ , and $\boldsymbol m$ is magnetisation, which is the magnetic moment per unit volume. For ease of analysis and the range in which the magnetic field is applied $(0.5{-}1.5 \, \textrm { T})$ , the interface response is assumed to be weakly magnetic. This simplifies the expression for the free energy of the interface to a sum of elastic energy and magnetic energy as

(2.1) \begin{equation} W=W_e(\boldsymbol{F})-\frac {\mu _0\chi }{2}{\boldsymbol h}\boldsymbol{\cdot }{\boldsymbol h}. \end{equation}

Here, $\mu _0$ is the permeability of free space, $\chi$ is the interface susceptibility and $\boldsymbol{h}$ is the induced magnetic field at the interface. This simplifies the constitutive equations of interface stresses and magnetisation as

(2.2) \begin{equation} {\boldsymbol \sigma }=\frac {\partial W}{\partial \boldsymbol F}{\boldsymbol F}^{T}-q{\boldsymbol I},\quad {\boldsymbol m}=-\frac {1}{\mu _0}\frac {\partial W}{\partial \boldsymbol h}=\chi {\boldsymbol h}. \end{equation}

In the above equation, $q$ is a constraint pressure associated with the incompressibility condition of the membrane, ${\boldsymbol F}^{T}$ is the transformation of deformation gradient and $\boldsymbol I$ is the unit tensor. The induced field can be split as

(2.3) \begin{eqnarray} {\boldsymbol h}={\boldsymbol h}_a+{\boldsymbol h}_s, \end{eqnarray}

where ${\boldsymbol h}_a$ is the applied magnetic field and ${\boldsymbol h}_s$ is the self-field, which is the secondary field created due to the induced magnetic moment at the interface. Due to the weakly magnetic nature of the membrane, the self-field can be ignored. This simplifies the divergence of pondermotive stress $\boldsymbol \sigma _{\textit{pon}}$ , which is the sum of the magnetic stress $\boldsymbol{\sigma }_{\textit{mag}}$ and Maxwell stress $\boldsymbol{\sigma }_{\textit{max}}$ , as

(2.4) \begin{align} \boldsymbol{\sigma }_{\textit{max}}&=\mu _0(\boldsymbol{h\otimes h})-\frac {\mu _0(\boldsymbol{h}\boldsymbol{\cdot }\boldsymbol{h})}{2} \boldsymbol{I},\nonumber\\[4pt]\boldsymbol{\sigma }_{\textit{mag}}&=\mu _0\boldsymbol{(h\otimes m)},\nonumber\\[4pt]\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }_{\textit{pon}}&=\mu _0\boldsymbol{\nabla h}\boldsymbol{\cdot }\boldsymbol{m}\approx \mu _0\chi \boldsymbol{\nabla h}_a\boldsymbol{\cdot }\boldsymbol{h}_a, \end{align}

where $\boldsymbol{\nabla}$ is the del operator in the current configuration. The three-dimensional magneto-elastic equation in the absence of applied body forces in the current configuration is

(2.5) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{T}=0,\quad \boldsymbol{T}=\boldsymbol{\sigma }+\boldsymbol{\sigma }_{\textit{pon}}. \end{align}

Using the expression for the divergence of pondermotive stress from (2.4) in (2.5) and transforming the interface stress to the reference configuration modifies the above equation as follows:

(2.6) \begin{equation} \boldsymbol{\bar {\boldsymbol{\nabla }}}\boldsymbol{\cdot }\boldsymbol{P}+\mu _0\chi \boldsymbol{\nabla }\boldsymbol{h}_a\boldsymbol{\cdot }\boldsymbol{h}_a=0, \end{equation}

where $\boldsymbol{P}$ is the first Piola stress tensor and $\boldsymbol{\bar {\boldsymbol{\nabla }}}$ is del operator in reference configuration. The above equation is also valid at the mid-surface of the membrane

(2.7) \begin{equation} \boldsymbol{\bar {\boldsymbol{\nabla }}}\boldsymbol{\cdot }\boldsymbol{P}_{\kern-2pt s}+\mu _0\chi \boldsymbol{\nabla }\boldsymbol{h}_a\boldsymbol{\cdot }\boldsymbol{h}_a=0, \end{equation}

with $\boldsymbol{P}_{\kern-2pt s}$ as the Piola stress evaluated at the membrane mid-surface. Due to the thinness and the very low thickness-to-radius ratio of the bubble surface, of approximately $\mathcal{O}(10^{-3})$ , the interface can be modelled as a membrane (Barham et al. Reference Barham, Steigmann, McElfresh and Rudd2007, Reference Barham, Steigmann and White2012). Let the membrane thickness be $\epsilon$ , then divergence of the first Piola stress becomes

(2.8) \begin{equation} \boldsymbol{\bar {\boldsymbol{\nabla }}}\boldsymbol{\cdot }\boldsymbol{P}_{\kern-2pt s}=\boldsymbol{\bar {\boldsymbol{\nabla }}}_s\boldsymbol{\cdot }\boldsymbol{P}_{\kern-2pt s}+\boldsymbol{P}_{\kern-2pt s}^{\prime}\boldsymbol{k}. \end{equation}

Here, $\boldsymbol{\bar {\boldsymbol{\nabla }}}_s$ is the two-dimensional surface gradient operator in the reference/undeformed configuration, $\boldsymbol{P}_{\kern-2pt s}^{\prime}$ is the differentiation of first Piola stress in the normal direction and $\boldsymbol{k}$ is normal to the reference surface. Now, first expanding the Piola stress in the thickness direction

(2.9) \begin{align} \boldsymbol{P^{+}}&=\boldsymbol{P}_{s}+\frac {\epsilon }{2}\boldsymbol{P}^{\prime}_{s},\nonumber\\ \boldsymbol{P^{-}}&=\boldsymbol{P}_{s}-\frac {\epsilon }{2}\boldsymbol{P}^{\prime}_{s}, \end{align}

where $\boldsymbol{P^{\pm }}$ are the limiting values of Piola stress at the top and bottom surfaces, respectively, from inside, $\boldsymbol{k^{\pm }}$ are corresponding exterior normals at the top and bottom $\boldsymbol{\pm k}$ . Adding the traction on top $\boldsymbol t_{a}^{+}$ and on the bottom surface $\boldsymbol t_{a}^{-}$ the following equations are obtained:

(2.10) \begin{align} \boldsymbol{P^{+}k^{+}}&=\alpha \,\boldsymbol{t}_a^{+},\nonumber\\ \boldsymbol{P^{-}k^{-}}&=\alpha \,\boldsymbol{t}_a^{-},\nonumber\\ \boldsymbol{P^{+}k^{+}+P^{-}k^{-}}&=\epsilon \,\boldsymbol{P}^{\prime}_s\boldsymbol{k}. \end{align}

Here, $\alpha$ is the local area dilatation. Substituting the above equations in (2.7) and rewriting in the current configuration, the following membrane equation is obtained:

(2.11) \begin{align} \boldsymbol{\nabla }_s\boldsymbol{\cdot }\boldsymbol{\sigma }_s+\alpha \frac {\boldsymbol{t}_a^+ +\boldsymbol{t}_a^-}{\epsilon }+\mu _0\boldsymbol{\nabla h}_a\boldsymbol{\cdot }\boldsymbol{h}_a=0. \end{align}

Here, $\boldsymbol{\nabla }_s$ is the two-dimensional gradient in the current/deformed configuration and $\sigma _{s}$ is the Cauchy stress at mid-surface of the membrane. In (2.2) the expression for the stress can be represented in terms of principal stretches $\lambda _i, i=1,2,3$ as

(2.12) \begin{eqnarray} \boldsymbol{\sigma }_s=\sum _{i=1}^{3}\left (\lambda _{i}\frac {\partial W}{\partial \lambda _i}-q\right )\boldsymbol{v}_i\otimes \boldsymbol{v}_i, \end{eqnarray}

and the deformation gradient as

(2.13) \begin{eqnarray} \boldsymbol{F}=\sum _{i=1}^{3}\lambda _i\,\boldsymbol{v}_i\otimes \boldsymbol{u}_i, \end{eqnarray}

where $\boldsymbol{u}_i$ and $\boldsymbol{v}_i$ are the principal stress directions before and after the deformation and $\lambda _i$ is the stretch in the corresponding principal stress direction. The incompressibility assumption of the bubble implies $\textrm {det}({\boldsymbol F})=1$ , which leads to $\lambda _3=(\lambda _1\lambda _2)^{-1}$ . A generalised plane stress condition is assumed by which the stress vector component $\boldsymbol{\sigma }_s\boldsymbol{\cdot }\boldsymbol{\hat {n}}$ normal to the membrane mid-surface vanishes. This leads the mechanical part of the stress component in the normal direction to be zero and the expression for the constraint pressure can be written as

(2.14) \begin{align} q=\lambda _3\frac {\partial W}{\partial \lambda _3}. \end{align}

2.2. Kinematics of bubble surface

Figure 1.

(a) Deformed and (b) undeformed configurations of a bubble surface.

Figure 1. Long description

The diagram consists of two parts labeled (a) and (b), representing the deformed and undeformed configurations of a bubble surface, respectively. In part (a), the bubble is shown in a deformed state with labeled components including R0, Theta, E Theta, and E R0. Arrows indicate the direction of forces and movement. In part (b), the bubble is shown in an undeformed state with labeled components including r(Theta, t), Theta(Theta, t), e theta, and e r. Arrows again indicate the direction of forces and movement. Both parts include a rotational arrow indicating the direction of phi. The diagram illustrates the mechanical properties and dynamic behavior of encapsulated magnetic microbubbles, highlighting the influence of magnetic nanoparticles on their oscillations and interactions within a magnetic field.

In this work, only axisymmetric coils with uniform current distribution are considered, which creates axisymmetric magnetic fields, and the applied pressure is radial. Consequently, only the axisymmetric oscillations of the bubble are taken into account. Considering a spherical bubble of initial radius $R_0$ , $(R_0,\varTheta )$ in the reference configuration gets transformed to $(r,\theta )$ in the current configuration, as shown in figure 1. The position of the bubble surface after the deformation is expressed in terms of the spherical mode shapes as

(2.15) \begin{align} r(\varTheta ,t)&=R(t)+\!\sum _{k=2}^{\infty }a_k(t)P_k(\cos {\varTheta }),\nonumber\\[3pt]\theta (\varTheta ,t)&=\varTheta +\!\sum _{k=1}^{\infty }\frac {b_k(t)}{R(t)}P^1_k(\cos {\varTheta }), \end{align}

with the radial mode shape represented by the Legendre polynomial $P_k(\cos {\varTheta })$ and the tangential mode shape represented by the associated Legendre polynomial $P^1_k(\cos {\varTheta })$ . In (2.15) $R(t),a_k(t)$ and $b_k(t)$ are the deformed radius and the amplitude of the radial and tangential mode shapes, respectively. The deformed and undeformed mid-surfaces of encapsulation are

(2.16) \begin{align} {\boldsymbol x}&=r(\varTheta ,t) \cos (\phi )\sin (\theta ){\boldsymbol{e}}_1+r(\varTheta ,t) \sin (\phi )\sin (\theta ){\boldsymbol e}_2+r(\varTheta ,t) \cos (\theta ){\boldsymbol e}_3=r(\varTheta ,t)\boldsymbol{e}_r,\nonumber\\[5pt]{\boldsymbol X}&=R_0 \cos (\phi )\sin (\varTheta ){\boldsymbol E}_1+R_0 \sin (\phi )\sin (\varTheta ){\boldsymbol E}_2+R_0 \cos (\varTheta ){\boldsymbol E}_3=R_0\boldsymbol{E}_{R_{0}}, \end{align}

where $\boldsymbol E_i$ and $\boldsymbol e_i$ are the Cartesian unit bases in the undeformed and deformed interfaces, respectively. Additionally, $\boldsymbol E_{R_0}$ and $\boldsymbol e_{r}$ are the unit vectors in the radial direction on the undeformed and deformed interface, respectively. The curvilinear basis vectors on the deformed surface are

(2.17) \begin{align} \boldsymbol{g}_1&=\boldsymbol{x}_{,\varTheta }=r(\varTheta ,t)_{,\varTheta }\,\boldsymbol{e_r}+r(\varTheta ,t)\theta _{,\varTheta }\,\boldsymbol{e_{\theta }},\nonumber\\[3pt] \boldsymbol{g}_2&=\boldsymbol{x}_{,\phi }=r(\varTheta ,t)\sin {\theta }\,\boldsymbol{e_{\phi }}. \end{align}

The undeformed basis vectors are given by

(2.18) \begin{align} \boldsymbol{G}_1&=\boldsymbol{X}_{,\varTheta }=R_{0}\boldsymbol{E}_{\varTheta }, \quad \quad \quad \boldsymbol{G}^1=\frac {\boldsymbol{E}_{\varTheta }}{R_0},\nonumber\\[4pt] \boldsymbol{G}_2&=\boldsymbol{X}_{,\phi }=R_0\sin {\phi }\boldsymbol{E}_{\phi }, \quad \boldsymbol{G}^2=\frac {\boldsymbol{E}_{\phi }}{R_0\sin {\phi }}, \end{align}

where the lower and upper subscripts correspond to co-variant and contra-variant vectors. The gradient on the deformed mid-surface can be written as

(2.19) \begin{align} \boldsymbol{\bar {\boldsymbol{\nabla }}x}&=\boldsymbol{g}_{\alpha }\otimes \boldsymbol{G}^{\alpha }, \\[-10pt] \nonumber \end{align}

(2.20) \begin{align} \boldsymbol{\bar {\boldsymbol{\nabla }}x}&=\frac {r(\varTheta ,t)_{,\varTheta }}{R_0}(\boldsymbol{e_r}\otimes \boldsymbol{E_\theta })+\frac {r(\varTheta ,t)\theta _{,\varTheta }}{R_0}(\boldsymbol{e_\theta } \otimes \boldsymbol{E_\varTheta }) +\frac {r(\varTheta ,t)\sin {\theta }}{R_0\sin {\varTheta }}(\boldsymbol{e_{\phi }} \otimes \boldsymbol{E_{\phi }}),\nonumber\\[5pt]&=\frac {\sqrt {(r(\varTheta ,t)_{,\varTheta })^2+(r(\varTheta ,t)\theta _{,\varTheta })^2}}{R_0}(\boldsymbol{\hat {t}}\otimes \boldsymbol{E_{\theta }})+\frac {r(\varTheta ,t)\sin {\theta }}{R_0\sin {\varTheta }}(\boldsymbol{e_{\phi }} \otimes \boldsymbol{E_{\phi }}),\nonumber\\[5pt]\boldsymbol{F}&=\boldsymbol{\bar {\boldsymbol{\nabla }}x}+\frac {1}{\lambda _1\lambda _2}(\hat {\boldsymbol n}\otimes \boldsymbol{e}_r). \end{align}

Comparing the above equations with (2.13), the principal stretches in the meridional and circumferential directions become

(2.21) \begin{align} \lambda _1&=\frac {\sqrt {(r(\varTheta ,t)_{,\varTheta })^2+(r(\varTheta ,t)\theta _{,\varTheta })^2}}{R_0}\approx \frac {r\theta ^{\prime }}{R_0}\approx \frac {R(t)}{R_0}+\frac {a_k(t)}{R_0}P_k+\frac {b_k(t)}{R_0}\frac {\partial P^1_k}{\partial \varTheta },\nonumber\\[4pt]\lambda _2&=\frac {r(\varTheta ,t)\sin (\theta )}{R_0\sin (\varTheta )}\approx \frac {R(t)}{R_0}+\frac {a_k(t)}{R_0}P_k+\frac {b_k(t)}{R_0}\cot (\varTheta )P^1_k. \end{align}

Similarly, the normal $(\boldsymbol{\hat {n}})$ and the tangent $(\boldsymbol{\hat {t}})$ vectors on the deformed surface are

(2.22) \begin{align} \boldsymbol{\hat {n}}\approx \boldsymbol{e}_{r}-\frac {a_k(t)}{R(t)}P^1_k \boldsymbol{e}_{\theta },\nonumber\\[3pt]\boldsymbol{\hat {t}}\approx \boldsymbol{e}_{\theta }+\frac {a_k(t)}{R(t)}P^1_k \boldsymbol{e}_{r}. \end{align}

From the expression for the normal, the principal curvature of the deformed surface is obtained as

(2.23) \begin{align} \boldsymbol{\kappa }&=\boldsymbol{\nabla \hat {n}}=\kappa _1\,\boldsymbol{e}_{\theta }\otimes \boldsymbol{e}_{\theta }+\kappa _2\,\boldsymbol{e}_{\phi }\otimes \boldsymbol{e}_{\phi }, \end{align}

where the principal curvature components $\kappa _1$ and $\kappa _2$ are given by

(2.24) \begin{align} \kappa _1&\approx \frac {1}{R(t)}-\sum _{k=2}^{\infty }\frac {a_k(t)}{R(t)^2}\left (P_k+\frac {\textrm {d}^2P_k}{\textrm {d}\theta ^2}\right )\!,\nonumber\\[3pt]\kappa _2&\approx \frac {1}{R(t)}-\sum _{k=2}^{\infty }\frac {a_k(t)}{R(t)^2}\left (P_k+\frac {\cos (\theta )}{\sin (\theta )}P^1_k\right )\!. \end{align}

The membrane is assumed to follow Mooney–Rivlin constitutive law as follows:

(2.25) \begin{equation} W_e=C_1\left (\lambda _1^2+\lambda _2^2+\frac {1}{\lambda _1^2\lambda _2^2} -3\right ) +C_2\left (\frac {1}{\lambda _1^2}+\frac {1}{\lambda _2^2}+\lambda _1^2\lambda _2^2-3\right )\!. \end{equation}

The in-plane mechanical stress becomes

(2.26) \begin{align} \sigma _{1}&=\lambda _1\frac {\partial W_e}{\partial \lambda _1}=2C_1\bigg (\lambda _1^2-\frac {1}{\lambda _1^2\lambda _2^2}\bigg )+2C_2\bigg (\lambda _1^2\lambda _2^2-\frac {1}{\lambda _1^2}\bigg ),\nonumber\\[4pt]\sigma _{2}&=\lambda _2\frac {\partial W_e}{\partial \lambda _2}=2C_1\bigg (\lambda _2^2-\frac {1}{\lambda _2^2\lambda _1^2}\bigg )+2C_2\bigg (\lambda _2^2\lambda _1^2-\frac {1}{\lambda _2^2}\bigg ). \end{align}

The divergence of the in-plane stress becomes (Pozrikidis Reference Pozrikidis2001)

(2.27) \begin{align} (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma })_n&=-\kappa _1\sigma _1-\kappa _2\sigma _2,\nonumber\\[3pt](\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma })_t&=\frac {\partial \sigma _1}{\partial s}+\frac {1}{d}\frac {\partial d}{\partial s}(\sigma _1-\sigma _2), \end{align}

where $d$ is the horizontal distance of a point in the membrane from the axis of revolution of bubble

(2.28) \begin{equation} d=R(t)\sin (\theta )+a_k(t)\sin (\theta )P_k+b_k(t)\cos (\theta )P^1_k. \end{equation}

Substituting the stretches and $x$ in the divergence expression, we get

(2.29) \begin{align} F^{\textit{mem}}_n&=\frac {\left (\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol \sigma }_s\right )_n}{\lambda _1\lambda _2} \approx {\mathcal{S}}_0+{\mathcal{S}}_1a_k(t)+{\mathcal{S}}_2b_k(t),\nonumber\\[3pt]F^{\textit{mem}}_t&=\frac {\left (\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol \sigma }_s\right )_t}{\lambda _1\lambda _2}\approx {\mathcal{S}}^1_1a_k(t)+{\mathcal{S}}^1_2b_k(t), \\[-12pt] \nonumber \end{align}

(2.30) \begin{align} \mathcal{S}_0&=C_1\left (\frac {-4R_0^6}{R^{7}}+\frac {4}{R}\right )+C_2\left (-\frac {4R_0^{4}}{R^{5}}+\frac {4R}{R_0^2}\right )\!,\nonumber\\[3pt]\mathcal{S}_1&=C_1\left (-(2k^2+2k-28)\frac {R_0^6}{R^8}+\frac {2k^2+2k-4}{R^2}\right )\nonumber\\[3pt]&\quad +C_2\left (-(2k^2+2k-20)\frac {R_0^4}{R^6}+\frac {2k^2+2k+4}{R_0^2}\right )\!,\nonumber\\[3pt]\mathcal{S}_2&=C_1\left (-12(k^2+k)\frac {R_0^6}{R^8}\right )+C_2\left (-(8k^2+8k)\frac {R_0^4}{R^6}-\frac {4k^2+4k}{R_0^2}\right )\!,\nonumber\\[3pt]{\mathcal{S}}^{1}_{1}&=C_1\left (-\frac {8R_0^6}{R^8}-\frac {4}{R^2}\right )+C_2\left (-\frac {4R_0^4}{R^6}-\frac {8}{R_0^2}\right )\!,\nonumber\\[3pt]{\mathcal{S}}^{1}_{2}&=C_1\left ((4k^2+4k)\frac {R_0^6}{R^8}+\frac {4k^2+4k-4}{R^2}\right )+C_2\left ((4k^2+4k-4)\frac {R_0^4}{R^6}+\frac {4k^2+4k}{R_0^2}\right )\!. \end{align}

2.3. Magnetic field and forces

Figure 2.

Schematics for magnetic field. (a) Two coils placed symmetrically about the ${\boldsymbol e}_1-{\boldsymbol e}_2$ plane carrying currents in opposite directions. (b) Two dipoles placed symmetrically about the $\boldsymbol{e}_\rho -\boldsymbol{e}_\phi$ plane.

Figure 2. Long description

The image consists of two diagrams labeled (a) and (b). Diagram (a) shows two coils placed symmetrically about a plane, with each coil carrying currents in opposite directions. The coils are labeled with points O and A, and vectors e2, e3, and I are indicated. Diagram (b) depicts two dipoles placed symmetrically about the same plane. The dipoles are aligned along the ez axis, with vectors eρ, eφ, and M shown. The diagrams illustrate the arrangement and interaction of magnetic fields generated by the coils and dipoles.

2.3.1. Coil field

Two coils are symmetrically placed above and below a bubble, each carrying current $I$ in opposite directions and having $N_1$ turns, as shown in figure 2(a). While the Biot–Savart law offers a straightforward approach to calculating the resulting magnetic field, it often involves the use of complex elliptic integrals, making the computation cumbersome. To simplify the analysis, an alternative method proposed by Boridy (Reference Boridy1989) is used to calculate magnetic fields for axially symmetric systems. The generalised expression for an axisymmetric magnetic fields is

(2.31) \begin{align} \boldsymbol{h}_a&=h_{\textit{ar}}\boldsymbol{e}_r+h_{a\theta }\boldsymbol{e}_{\theta },\nonumber\\h_{\textit{ar}}&=\sum _{n=1}^{\infty }H_n\left (\frac {r}{a}\right )^{n-1}P_{n}(\cos {\theta }),\nonumber\\ h_{a\theta }&=\sum _{n=1}^{\infty }\frac {H_n}{n}\left (\frac {r}{a}\right )^{n-1}P_{n}^{1}(\cos {\theta }), \end{align}

where $r$ and $\theta$ are radial and angular position of any point in the $r{-}z$ plane, $a$ is the distance from the bubble $O$ to any point on the coil circumference $A$ and $\theta _c$ is the angle between the axis of symmetry $\boldsymbol{e}_3$ and the vector along $OA$ . For the case of coils carrying currents in opposite directions placed symmetrically above and below the equator, the odd terms vanish in (2.31) and only even terms are sustained. The expressions for even coefficients are given as follows:

(2.32) \begin{equation} H_{2n}=-\frac {I}{a}N_{1}\sin {\theta _c}\,P^{1}_{2n}(\cos {\theta _c}). \end{equation}

By choosing the value of $\theta _c$ as a positive root of $P^{1}_4$ , the coefficient $H_4$ can be made zero. As we are interested in fields and field gradients in the region of the bubble, the $r/a$ ratio is very small. Hence, the subsequent terms $H_{6}, H_{8},\ldots$ in (2.31) can be ignored and only leading terms can be considered, resulting in simplified expressions for the magnetic fields as

(2.33) \begin{align} h_{\textit{ar}}&=H_{2}\frac {r}{2a} \big(3\cos ^2{\theta }-1 \big),\nonumber\\[4pt]h_{a\theta }&=H_{2}\frac {r}{2a}(-3\sin {\theta }\cos {\theta }). \end{align}

In spherical coordinates the magnetic forces become

(2.34) \begin{align} F^{\textit{mag}}_{r}&=(\mu _0\chi \boldsymbol{\nabla }\boldsymbol{h}_{a}\boldsymbol{\cdot }\boldsymbol{h}_a)_r=\mu _0\chi \left (h_{\textit{ar}}\frac {\partial h_{\textit{ar}}}{\partial r}+\frac {h_{a\theta }}{r}\frac {\partial h_{\textit{ar}}}{\partial \theta }-\frac {h_{a\theta }^2}{r}\right )\!,\nonumber\\ &=\mu _0\chi \frac {H_{2}^{2}r}{4a^2}(3\cos ^2{\theta }+1),\nonumber\\ F^{\textit{mag}}_{\theta }&=(\mu _0\chi \boldsymbol{\nabla }\boldsymbol{h}_{a}\boldsymbol{\cdot }\boldsymbol{h}_a)_\theta =\mu _0\chi \left (h_{\textit{ar}}\frac {\partial h_{a\theta }}{\partial r}+\frac {h_{a\theta }}{r}\frac {\partial h_{a\theta }}{\partial \theta }+\frac {h_{\textit{ar}}h_{a\theta }}{r}\right )\!,\nonumber\\&=\mu _0\chi \frac {H_{2}^{2}r}{4a^2}(-3\sin {\theta }\cos {\theta }). \end{align}

2.3.2. Dipole field

The expression for the magnetic field generated due to a dipole is (Latulippe & Martel Reference Latulippe and Martel2014)

(2.35) \begin{align} \boldsymbol{h}_a=M\left [\frac {3(\boldsymbol{r}\boldsymbol{\cdot }{\boldsymbol{e}}_z)\boldsymbol{r}}{\|{\boldsymbol{r}}\|^5}-\frac {{\boldsymbol{e}}_z}{\|{\boldsymbol{r}}\|^3}\right ]\!. \end{align}

Here, $M$ is the magnetic dipole strength, $\boldsymbol{r}$ is the vector position of the bubble surface from the dipole and $\|{\boldsymbol{\cdot }}\|$ denotes the Euclidean norm. Assuming that the dipole is placed sufficiently far away from the bubble, $z\gg x$ and $z\gg y$ , the above expression for the magnetic field gets simplified as

(2.36) \begin{align} \boldsymbol{h}_a=\frac {3Mx}{z^4}\boldsymbol{e}_\rho +\frac {2M}{z^3}\boldsymbol{e}_z. \end{align}

As the two dipoles are place symmetrically about $\boldsymbol{e}_r$ , the field components along $\boldsymbol{e}_z$ cancel and the component along $e_r$ adds up, resulting in

(2.37) \begin{align} \boldsymbol{h}_a=\frac {6Mx}{z^4}\boldsymbol{e}_\rho . \end{align}

The forcing on the bubble surface due to the dipole field is

(2.38) \begin{align} \mu _0\chi _s\boldsymbol{\nabla h}_a\boldsymbol{\cdot }\boldsymbol{h}_a&=\frac {36\mu _0\chi _sM^2x}{z^8}\boldsymbol{e}_\rho , \end{align}

(2.39) \begin{align} \boldsymbol{F}^{mag}&=\left (\frac {1}{\lambda _1\lambda _2}\right )\frac {36\mu _0\chi M^2}{z^8}r\sin {\theta }\boldsymbol{e}_\rho . \end{align}

The calculations for the dipole field are made using cylindrical coordinates, as shown in figure 2(b). Substituting for $r$ and $\theta$ the deformed radial and meridional position from (2.15) in the above equations, the radial and meridional components of the magnetic force are obtained. Once these components are obtained in spherical coordinates, the dot product of the magnetic forces with the normal and tangential vectors is calculated. This gives us the corresponding force components in the normal and tangential directions as

(2.40) \begin{align} & F^{\textit{mag}}_{n} =\frac {1}{\lambda _1\lambda _2}\left (F^{\textit{mag}}_{r}-\frac {a_k(t)P^1_k}{R(t)}F^{\textit{mag}}_{\theta }\right )\!,\nonumber\\F^{\textit{mag}}_{t}&=\frac {1}{\lambda _1\lambda _2}\left (F^{\textit{mag}}_{\theta }+\frac {a_k(t)P^1_k}{R(t)}F^{\textit{mag}}_{r}\right )\!. \end{align}

The above equations are linearised in $a_k$ and $b_k$ to be consistent with the present formulation. The final equations, after using the recurrence relations of Legendre polynomials, are

(2.41) \begin{align} & F^{\textit{mag}}_{n} =\mathcal{H}_0+f_1b_1P_1+\mathcal{H}_1P_2+\mathcal{H}\left [\sum _{k=2}^{\infty }\mathcal{F}_{1}(k)a_kP_k +\!\sum _{k=2}^{\infty }\mathcal{F}_{2}(k)b_kP_k\right .\nonumber\\[4pt]\quad &+\left .\sum _{k=0}^{\infty }\mathcal{F}_{3}(k)a_{k+2}P_k +\!\sum _{k=0}^{\infty }\mathcal{F}_{4}(k)b_{k+2}P_k+\!\sum _{k=4}^{\infty }\mathcal{F}_{5}(k)a_{k-2}P_k+\!\sum _{k=4}^{\infty }\mathcal{F}_{6}(k)b_{k-2}P_k\right ]\!,\nonumber\\[4pt]& F^{\textit{mag}}_{t}=f^1_1b_1P^1_1+\mathcal{H}^1_2P^1_2+\mathcal{H}\left [\sum _{k=2}^{\infty }\mathcal{F}^{1}_{1}(k)a_kP^1_k +\!\sum _{k=2}^{\infty }\mathcal{F}^{1}_{2}(k)b_kP^1_k\right .\nonumber\\[4pt]&\quad +\left .\sum _{k=0}^{\infty }\mathcal{F}^{1}_{3}(k)a_{k+2}P^1_k+\!\sum _{k=0}^{\infty }\mathcal{F}^{1}_{4}(k)b_{k+2}P^1_k+\!\sum _{k=4}^{\infty }\mathcal{F}^{1}_{5}(k)a_{k-2}P^1_k+\!\sum _{k=4}^{\infty }\mathcal{F}^{1}_{6}(k)b_{k-2}P^1_k\!\right ]\!, \end{align}

where the expressions for $\mathcal{H}$ , $\mathcal{H}_i$ , $\mathcal{H}^1_i$ , $\mathcal{F}_i$ and $\mathcal{F}^1_i$ are given in Appendix A for both coil and dipole fields. It is evident from Liu et al. (Reference Liu, Sugiyama, Takagi and Matsumoto2012, Reference Liu, He, Gong and Huang2018) that an applied acoustic pressure primarily excites the radial mode. If we disregard any initial disturbances to the non-spherical modes, the applied magnetic field excites the radial mode as well as the $a_2$ and $b_2$ modes. Additionally, the forcing on the $k^{\textrm{}}$ th mode depends on the neighbouring modes, specifically the $k-2$ and $k+2$ modes. As a result, under the influence of the magnetic field generated by both the coil and dipole, only the even-numbered modes are triggered. This allows us to ignore the translations and odd-numbered modes in the equations for mode shapes.

2.4. Fluid equations

The assumption is made that the flow field inside the bubble is negligible and fluid flow is considered only outside the bubble, governed by the incompressible Navier–Stokes equations

(2.42) \begin{align} {\boldsymbol{\nabla }}\boldsymbol{\cdot }{\boldsymbol u}&=0,\nonumber\\ \rho \frac {\partial \boldsymbol{u}}{\partial t}+\rho [{\boldsymbol u}\boldsymbol{\cdot }{\boldsymbol{\nabla }}]\boldsymbol{u}&=-\boldsymbol{\nabla } p+\eta \boldsymbol{\nabla }\boldsymbol{\cdot }{(\boldsymbol{\nabla }u+\boldsymbol{\nabla }u^T)}, \end{align}

where $\boldsymbol{u}$ and $p$ are the fluid velocity and pressure, respectively. Here, $\rho$ and $\eta$ represent the density and viscosity of the fluid, respectively. The non-spherical oscillations are assumed to be small compared with the radial oscillations, and the velocity is split into potential and viscous flow (Prosperetti Reference Prosperetti1977) as follows:

(2.43) \begin{equation} \boldsymbol{u}=\boldsymbol{u}_p+\boldsymbol{u}_v. \end{equation}

The velocity of the potential flow is calculated by substituting $\eta =0$ . The generalised solution for the axisymmetric Laplace equation for incompressible flow is

(2.44) \begin{align} \psi &=\frac {C_0P_0}{r}+\!\sum _{k=2}^{\infty }\frac {C_kP_k}{r^{k+1}}, \\[-12pt] \nonumber \end{align}

(2.45) \begin{align} \boldsymbol{u}_p&=\frac {\partial \psi }{\partial r}\boldsymbol{e}_r+\frac {1}{r}\frac {\partial \psi }{\partial \theta }\boldsymbol{e}_\theta . \end{align}

The fluid velocity continuity at the interface in the normal direction is given by

(2.46) \begin{align} S&=r-R-a_kP_k,\nonumber\\ 0&=\frac {\partial S}{\partial t}+(\boldsymbol{\nabla} \psi )\boldsymbol{\cdot }\boldsymbol{\nabla } S. \end{align}

Expanding the terms in (2.46) gives

(2.47) \begin{align} -\dot {R}+\!\sum _{k=2}^{\infty } -\dot {a}_{k}P_k-\frac {C_0}{R^2}+\frac {2C_0a_kP_k}{R^3}-\frac {(k+1)C_kP_k}{R^{k+2}}=0. \end{align}

The coefficients in (2.47) are determined by using the orthogonality property of Legendre polynomials. These conditions are given by collecting

zeroth-order coefficient

(2.48) \begin{equation} \dot {R}+\frac {C_0}{R^2}=0,\nonumber\\ C_0=-R^2\dot {R}, \end{equation}

the kth-order coefficient

(2.49) \begin{align} -\dot {a}_k&+\frac {2C_0a_k}{R^3}-\frac {C_k(k+1)}{R^{k+2}}=0,\nonumber\\[5pt]C_k&=-\frac {R^{k+2}}{k+1}\left (\dot {a}_k+\frac {2a_k\dot {R}}{R}\right )\!. \end{align}

In the Navier–Stokes equation viscosity is substituted as zero and the equation integrated by substituting the potential flow velocity, which gives pressure developed due to potential flow as

(2.50) \begin{equation} p_p=p_\infty -\rho \left (\frac {\partial \psi }{\partial t}+\frac {1}{2}|\boldsymbol{\nabla }\psi |^2\right )\!. \end{equation}

The viscous correction to the potential flow is obtained by solving the linearised Navier–Stokes equation in velocity–vorticity form

(2.51) \begin{align} \frac {\partial \boldsymbol{\omega }}{\partial t}+\boldsymbol{(\boldsymbol{\nabla }\times \omega \times u)}&=-\nu \boldsymbol{\nabla \times (\boldsymbol{\nabla }\times \omega )},\nonumber\\ \boldsymbol{\omega }&=\boldsymbol{\nabla \times u}. \end{align}

Due to axisymmetric motion of the bubble, only the toroidal component of vortex is considered as

(2.52) \begin{align} \boldsymbol{\omega }=\boldsymbol{\nabla \times }(T_kP_k\boldsymbol{e}_r). \end{align}

Here, $T_k$ is the toroidal component of vorticity (Prosperetti Reference Prosperetti1977), substituting the above in the Navier–Stokes equation gives (Liu et al. Reference Liu, Sugiyama, Takagi and Matsumoto2012, Reference Liu, He, Gong and Huang2018)

(2.53) \begin{align} \rho \frac {\partial T_k}{\partial t}+\rho \frac {\partial }{\partial r}\left (\frac {\dot {R}R^2T_k}{r^2}\right )-\eta \frac {\partial ^2T_k}{\partial r^2}+\eta k(k+1)r^{-2}T_k=0. \end{align}

The viscous velocity and viscous pressure are obtained as

(2.54) \begin{align} \boldsymbol{u}_v&=\left (\sum _{k=2}^{\infty }T_kP_k-\frac {\partial \varPhi }{\partial r}\right ){\boldsymbol e_r}-\frac {1}{r}\frac {\partial \varPhi }{\partial \theta }\boldsymbol{e}_\theta , \\[-12pt] \nonumber \end{align}

(2.55) \begin{align} p_v&=\sum _{k=1}^{\infty }k\left [\mu \frac {T_k}{R}+\frac {\rho \dot {R}}{R}\int _{R}^{\infty }\left [\frac {R^3}{s^3}-1\right ]\frac {R^k}{s^k}T_k \,\textrm {d}s\right ]P_k, \\[-12pt] \nonumber \end{align}

(2.56) \begin{align} \varPhi &=\sum _{k=2}^{\infty }P_k\bigg [\bigg (-\frac {k+1}{2k+1}\int _{R}^{\infty }s^{-k}T_k \,\textrm {d}s+\frac {k+1}{2k+1}\int _{R}^{r}s^{-k}T_k \,\textrm {d}s\bigg )r^{k}\nonumber \\& \quad -\bigg (\frac {k}{2k+1}R^{2k+1}\int _{R}^{\infty }s^{-k}T_k\,\textrm {d}s -\frac {k}{2k+1}\int _{R}^{r}s^{k+1}T_k\,\textrm { d}s\bigg )r^{-(k+1)}\bigg ]. \end{align}

The viscous stress due to this flow on the bubble surface becomes

(2.57) \begin{equation} \boldsymbol{T^{\textit{fluid}}_{\textit{ex}}}=\eta (\boldsymbol{\nabla u+\boldsymbol{\nabla }u^T}). \end{equation}

The components of this flow in the normal and tangential directions become

(2.58) \begin{align} \boldsymbol{\hat {n}}\boldsymbol{\cdot }\boldsymbol{T}^{\textit{fluid}}_{\textit{ex}}\boldsymbol{\cdot }\boldsymbol{\hat {n}}&=\frac {-4\eta \dot {R}}{R}+\!\sum _{k=2}^{\infty }2\eta \bigg [(k+2)\frac {\dot {a}_k}{R}-2(k-1)\frac {\dot {R}}{R^2}a_k\bigg .\nonumber \\&\quad \bigg .+k(k+1)R^{k-2}\int _{R}^{\infty }s^{-k}T_k\,\textrm { d}s\bigg ]P_k, \\[-12pt] \nonumber \end{align}

(2.59) \begin{align} \boldsymbol{\hat {t}}\boldsymbol{\cdot }\boldsymbol{T}^{\textit{fluid}}_{\textit{ex}}\boldsymbol{\cdot }\boldsymbol{\hat {n}} &=\sum _{k=2}^{\infty }2\eta \bigg [\frac {k+2}{k+1}\frac {\dot {a}_k}{R}+\frac {1-k}{k+1}\frac {\dot {R}a_k}{R^2}-\frac {T_k}{2R}-R^{k-2}\int _{R}^{\infty }s^{-k}T_k \,\textrm {d}s\bigg ]P^1_k. \end{align}

2.5. Final equations

The traction on the lower and upper surfaces of the membrane are

(2.60) \begin{align} \boldsymbol{t}_a^+&=-p_{\textit{ex}}\boldsymbol{\hat {n}}+\boldsymbol{T}^{\textit{fluid}}_{\textit{ex}}\boldsymbol{\hat {n}},\nonumber\\ \boldsymbol{t}_a^-&=\,\,\,p_{\textit{in}}\boldsymbol{\hat {n}}. \end{align}

Dividing (2.11) by $-\alpha$ and substituting $F^{\textit{mem}}$ , $F^{\textit{mag}}$ , $\boldsymbol{t}^{+}_a$ and $\boldsymbol{t}^{-}_{a}$ in the governing equation stress balance at the interface, the following equations are obtained:

(2.61) \begin{align} F^{\textit{mem}}_{n}+\frac {p_{\textit{ex}}}{\epsilon }-\frac {p_{\textit{in}}}{\epsilon }-\frac {\boldsymbol{\hat {n}}\boldsymbol{\cdot }\boldsymbol{T^{\textit{fluid}}_{\textit{ ex}}}\boldsymbol{\hat {n}}}{\epsilon }-F^{\textit{mag}}_{\textit{rad}}=0,\nonumber\\ F^{\textit{mem}}_{t}-\frac {\boldsymbol{\hat {t}}\boldsymbol{\cdot }\boldsymbol{T^{\textit{fluid}}_{\textit{ex}}}\boldsymbol{\hat {n}}}{\epsilon }-F^{\textit{mag}}_{\textit{tan}}=0. \end{align}

From the expression for the divergence of stresses at the interface, tractions from above and below it can be observed to be well separated in terms of Legendre or associated Legendre polynomials. The orthogonality property of these polynomials can be utilised to decompose the equation into separate modes

(2.62) \begin{equation} \left . \begin{array}{ll} \displaystyle \int _{-1}^{1}P_{n}^{k}P_m^k\,\text{d}x=0 \quad w\textit{hen}\quad [m\not = n],\nonumber\\[8pt] \displaystyle \int _{-1}^{1}P_{n}^{k}P_m^k\,\text{d}x=\frac {2}{2n+1}\frac {(n+m)!}{(n-m)!} \quad w\textit{hen}\quad [m=n]. \end{array}\!\!\right \} \end{equation}

Multiplying (2.61) $_1$ by $P_k$ and integrating the equations for radial, and kth mode shape yields the following equation:

(2.63) \begin{align} & \rho \left (\!R\ddot {R}+\frac {3}{2}\dot {R}^2\!\right )+p_\infty -p_{in}+4\eta \frac {\dot {R}}{R}+F^{\textit{ mem}}_{n0}-\mathcal{H}_0-\mathcal{H}[\mathcal{F}_3(0)a_2(t)+\mathcal{F}_4(0)b_2(t)]=0, \\[-28pt] \nonumber \end{align}

(2.64) \begin{align} &\rho \left (\frac {R\ddot {a}_2}{3}+\dot {R}\dot {a}_2-\frac {\ddot {R}}{3}a_2+\frac {2\dot {R}}{R}\int _{R}^{\infty }\left [\frac {R^3}{s^3}-1\right ]\frac {R^2}{s^2}T_2 \,\textrm {d}s\right ) +\frac {2\eta T_2}{R}\nonumber \\[4pt]&\qquad + {2\eta }\left [\frac {4\dot {a}_2}{R}+\frac {2\dot {R}a_2}{R^2}-6\int _{R}^{\infty }s^{-2}T_2 \,\textrm {d}s\right ]+F^{\textit{ mem}}_{n2}-\mathcal{H}_2\nonumber \\[4pt] &\qquad -\mathcal{H}[\mathcal{F}_1(2)a_2(t)+\mathcal{F}_2(2)b_2(t) +\mathcal{F}_3(2)a_4(t)+\mathcal{F}_4(2)b_4(t)]=0, \\[-28pt] \nonumber \end{align}

(2.65) \begin{align} \begin{split} &\rho \left (\frac {R\ddot {a}_k}{k+1}+\left [\frac {3\dot {R}}{k+1}+\frac {2(k+2)\eta }{\rho R}\right ]\dot {a}_k+\left [\frac {-(k-1)\ddot {R}}{k+1}+\frac {4(k-1)\eta }{\rho }\frac {\dot {R}}{R^2}\right ]a_k +\frac {k\eta T_k}{\rho R}\right .\\[4pt]&\left .+k\frac {\dot {R}}{R}\int _{R}^{\infty }\left [\frac {R^3}{s^3}-1\right ]\frac {R^k}{s^k}T_k \,\textrm { d}s -\frac {2\eta }{\rho }k(k+1)\int _{R}^{\infty }s^{-k}T_k \,\textrm {d}s\right ) +F^{\textit{ mem}}_{nk}-\mathcal{H}\left [\mathcal{F}_{1}(k)a_k(t)\right .\\[4pt]&\left .+\mathcal{F}_{2}(k)b_k(t)+\mathcal{F}_{3}(k)a_{k+2}(t)+\mathcal{F}_{4}(k)b_{k+2}(t)+\mathcal{F}_{5}(k)a_{k-2}(t)+\mathcal{F}_{6}(k) b_{k-2}(t)\right ]=0. \end{split} \\[-28pt] \nonumber \end{align}

Similarly, orthogonality conditions are used to simplify (2.61) $_2$ to obtain the equations for $b_k$ as

(2.66) \begin{align} &2\eta \left [\frac {k+2}{k+1}\frac {\dot {a}_k}{R}+\frac {1-k}{k+1}\frac {\dot {R}a_k}{R^2}-\frac {T_k}{2R}-R^{k-2}\int _{R}^{\infty }s^{-k}T_k \,\textrm {d}s\right ]+F^{\textit{ mem}}_{tk}-\mathcal{H}^1_2\delta _{k2}\nonumber\\[4pt]&\quad -\mathcal{H}[\mathcal{F}^{1}_{1}(k)a_k(t)+\mathcal{F}^{1}_{2}(k)b_k(t)+\mathcal{F}^{1}_{3}(k)a_{k+2}(t) +\mathcal{F}^{1}_{4}(k)b_{k+2}(t)\nonumber\\[4pt]&\quad +\mathcal{F}^{1}_{5}(k)a_{k-2}(t)+\mathcal{F}^{1}_{6}(k) b_{k-2}(t)]=0. \end{align}

Along with this, a no-slip boundary condition at the surface of the bubble in the tangential direction is used

(2.67) \begin{equation} u_{\theta }=r\frac {\partial \theta }{\partial t}. \end{equation}

Substituting the deformed radial and angular positions of the bubble surface from (2.15) in the above equation, it gets simplified as

(2.68) \begin{equation} -\frac {\dot {a}_k}{k+1}-\frac {2\dot {R}a_k}{(k+1)R}+R^{k-1}\int _{R}^{\infty }s^{-k}T_k \,\textrm {d}s=\dot {b}_k-\frac {\dot {R}b_k}{R}. \end{equation}

It can be observed from (2.63) to (2.65) that the radial oscillations are directly influenced by both the applied pressure and the magnetic field. In contrast, the $a_2$ shape mode is directly driven by the magnetic field, while its response to pressure is indirect, arising from the coupling between the radial and shape mode equations. All other even shape modes are similarly excited through the coupling between equations, rather than by direct external forcing.

5.1. Solution verification

To prove the convergence of the numerical scheme as no exact solution to the system of equations is available a priori, a solution obtained with a very fine time step $\bar {a}_i(t^{\prime})$ is used as the reference (LeVeque Reference LeVeque2007; Birken et al. Reference Birken, Quint, Hartmann and Meister2010). The solutions computed with coarser time steps are then compared against these reference data to demonstrate convergence. As the magnitudes of the mode shapes lie in the range of $\mathcal{O}(10^{-3} {-} 10^{-30})$ , only the relative error is considered. The system of equations has a fluid subsystem and a structural subsystem. Although the numerical algorithms for the structure and fluid are order four and order one, the resulting loosely coupled method is neither fourth-order nor first-order accurate. To calculate the order of accuracy of the stagger scheme, the following error measures are introduced:

(5.4) \begin{align} e_s=\,\underset {i}{\textrm {max}}\,\frac {\|{a_i(t^{\prime})-\bar {a}_i(t^{\prime})}\|_2}{\|{\bar {a}_i(t^{\prime})}\|_2},\quad \quad e_f=\,\underset {i}{\textrm {max}}\,\frac {\|{T^{100}_i-\bar {T}^{100}_i}\|_2}{\|{\bar {T}^{100}_i}\|_2}. \end{align}

Here, $\|{.}\|_2$ denotes the $2$ -norm, $T_i^{100}$ is the ith toroidal vortex at $100\,\unicode{x03BC} \textrm {s}$ and $a_i$ is ith mode amplitude. To test the implementation, simulations are carried out for a total duration of $100\,\unicode{x03BC} \text{s}$ , with the reference solution $\bar {a}_i$ obtained using a non-dimensional time step of $\Delta t^{\prime} = 1.7 \times 10^{-4}$ . The convergence behaviour is then assessed through the log–log plots shown in figure 3, where $\log$ refers to the common logarithm. It is observed that all the shape modes have qualitatively similar log–log plots. As expected, the maximum error is observed at higher mode numbers considered in the numerical analysis. This convergence in error proves that the time series converges to some solution. To ensure the correctness of the solution, the full-scale simulation is compared with the solution obtained by using the boundary-layer approximation, as shown in figure 4. It is observed that the oscillation patterns and mean mode amplitudes are in close agreement, with only a small deviation in the maximum mode amplitude between the boundary-layer approximation and the full-scale simulation. Based on these two results, it can be reasonably concluded that the numerical solution method provides correct results. Since the upwind scheme is used for discretising the convection term, it must satisfy the Courant–Friedrichs–Lewy condition $(u_{\textit{conv}}\Delta t/\Delta y\leqslant 1)$ . Because the scheme is only conditionally stable, the solution becomes numerically unstable and eventually diverges when larger time steps $\Delta t^{\prime}$ are used. The measured slope of the $\log {e_s}$ plot is $1.995$ , and for the $\log {e_f}$ plot is $2.065$ , based on which it can be concluded that the loosely coupled scheme is approximately second-order accurate. From the error plots shown in figure 3, $\Delta t^{\prime} = 1 \times 10^{-3}$ is chosen for our simulations.

Figure 3.

Error convergence plot for the (a) structural subsystem and (b) fluid subsystem.

Figure 3. Long description

Two line graphs showing error convergence for structural and fluid subsystems. The x-axis represents the logarithm of delta t prime, and the y-axis represents the logarithm of error. Each graph contains a blue line with red data points, indicating a linear relationship between the variables. All values are approximated.

Figure 4.

Comparison of full-scale simulation with solution obtained by boundary-layer approximation.

Figure 4. Long description

The image contains four line graphs comparing full-scale simulation with boundary-layer approximation for different coefficients over time. The x-axis represents time in seconds, ranging from 7.90 times 10 to the power of -4 to 8.00 times 10 to the power of -4. The y-axes represent the ratios of different coefficients to R0. The top left graph shows the ratio a2/R0, with values ranging from 0.01468 to 0.01482. The top right graph shows the ratio a4/R0, with values ranging from 6.50 times 10 to the power of -5 to 6.75 times 10 to the power of -5. The bottom left graph shows the ratio a6/R0, with values ranging from 3.45 times 10 to the power of -7 to 3.75 times 10 to the power of -7. The bottom right graph shows the ratio a8/R0, with values ranging from 1.60 times 10 to the power of -9 to 1.90 times 10 to the power of -9. The blue line represents the full-scale simulation, and the orange dotted line represents the approximation. All values are approximated.

5.2. Stability diagrams

In the study of the bubble dynamics under an applied magnetic field, all even-numbered shape modes are excited, making it challenging to determine the number of modes that need to be considered for analysis. To address this challenge, the pressure–frequency stability diagram (see figure 5(a)) is analysed as a function of the number of modes, while applying the pressure at infinity as $p_{\infty }=p_0(1+\epsilon _p\cos {\omega _dt}).$ Here, $\epsilon _p$ and $\omega _d$ represent the dimensionless acoustic pressure with ambient pressure and the circular frequency, respectively. The stability curve represents the critical pressure above which the amplitude response of at least one non-spherical mode $a_k$ changes its behaviour from bounded oscillations to exponential blowup. In the simulations, the coils are assumed to carry a current of $I = 100$ kA, with the number of coil turns set to $ N_1 = 1000$ , with a thickness of $\epsilon = 20$ nm and an initial bubble radius of $R_0 = 10$ $\unicode{x03BC}$ m.

Figure 5.

Convergence of the stability diagram with number of modes and results of linear stability analysis. (a) Stability diagram with total modes considered in the analysis and (b) mode number unstable at a particular frequency.

Figure 5. Long description

The line graph consists of two subplots labeled (a) and (b). Subplot (a) displays the stability diagram with total modes considered in the analysis. The x-axis is labeled with omega sub d over 2 pi times 10 to the power of 6, ranging from 0 to 6. The y-axis is labeled with epsilon sub p, ranging from 0 to 50. Five lines represent different values of k: 12, 14, 16, 24, and 32. Subplot (b) shows the mode number unstable at a particular frequency. The x-axis is labeled with omega sub d over 2 pi times 10 to the power of 6, ranging from 0 to 7. The y-axis is labeled with k, ranging from 0 to 25. A red line represents the mode number unstable at different frequencies. All values are approximated.

Figure 5(a) shows that, as the number of modes increases from 12 to 32, the critical pressure curves in the $\epsilon _p - \omega _d$ plane stabilise, defining a well-converged stable region where the linearised non-spherical oscillation equations remain valid. This convergence occurs because, in the absence of initial disturbances in the shape modes, the $a_2$ mode is initially dominant. The amplitude of each successive even-numbered mode gradually decreases, approaching nearly zero at higher mode numbers.

Furthermore, the stability analysis results discussed in § 4 indicate that the first mode to become unstable during a frequency sweep shifts to higher modes as the frequency increases, as shown in figure 5(b). This explains the convergence of the critical pressure curve when more than 16 modes are considered for driving frequencies up to $6$ MHz.

Figure 6.

(a) Variation of natural frequency with mode number plot, $\omega _k\,\,\textrm {vs}\,\,k$ , at various values of applied static current. (b) The critical pressure boundary in the pressure–frequency plane, $\epsilon _p\,\,\textrm {vs}\,\,\omega _d$ , shown for different values of applied static current. (c) Maximum eigenvalue of fundamental matrix variation with driving frequency at various values of applied current.

Figure 6. Long description

The line graph consists of three subplots labeled (a), (b), and (c). Subplot (a) shows the variation of natural frequency with mode number at various values of applied static current. The x-axis represents the mode number, and the y-axis represents the natural frequency. Different symbols and colors represent different values of applied static current: 50 kiloamperes, 100 kiloamperes, and 150 kiloamperes. Subplot (b) depicts the critical pressure boundary in the pressure-frequency plane for different values of applied static current. The x-axis represents the driving frequency divided by 2 pi, and the y-axis represents the critical pressure. Different line styles and colors represent different values of applied static current: 50 kiloamperes, 75 kiloamperes, 100 kiloamperes, 125 kiloamperes, and 150 kiloamperes. Subplot (c) shows the maximum eigenvalue of the fundamental matrix variation with driving frequency at various values of applied current. The x-axis represents the driving frequency divided by 2 pi, and the y-axis represents the maximum eigenvalue. Different line styles and colors represent different values of applied static current: 50 kiloamperes, 100 kiloamperes, and 150 kiloamperes. All values are approximated.

5.2.1. Effect of static current and material parameters on stability diagram for coil field

Figure 6(b) shows the influence of an applied static current on the critical pressure above which the mode amplitude response becomes unstable. From this, it is evident that the applied magnetic field does not affect the pressure–frequency stability diagram of surface instability. Similarly, the $\omega _k - k$ plot in figure 6(a), which shows the relation between natural frequency and mode number, suggests that the applied current does not noticeably influence the natural frequency. Additionally, the same plot reveals that several harmonics, sub-harmonics and super-harmonics of different modes are closely spaced around the fundamental frequency of a particular mode. As a result, when the critical pressure is exceeded, multiple modes become unstable simultaneously. This complexity makes stability analysis feasible only through computational methods. Figure 6(c) is obtained by calculating the maximum eigenvalue of the fundamental matrix $X(T)$ by varying the pressure amplitude until the eigenvalue becomes greater than one. The $\text{max}(\lambda _i)$ plot of figure 4(c) indicates that the maximum eigenvalue remains unaffected by increasing the current in the system of equations, which explains the invariance of the stability diagram with respect to variations in current. Furthermore, the frequency is non-dimensionalised as $\omega ^{\prime}=\sqrt {\rho _f R_0^3/2c_1}$ and ambient pressure as $P=p_0R_0/c_1$ . As shown in figure 7(a), this scaling causes the frequency curves to collapse onto a single curve when the undeformed radius $R_0$ is varied. On the other hand, the frequency curves do not collapse when $c_1$ is varied, as seen in figure 8(a). This is due to the presence of the $c_2$ term in the Mooney–Rivlin model, which introduces an additional non-dimensional parameter $\alpha ={c_2}/{c_1}$ that affects the scaling. As $c_1$ increases while $c_2$ remains fixed, $\alpha$ decreases, reducing its influence and causing the frequency curves to converge. The pressure–frequency stability curves do not collapse when either $R_0$ or $c_1$ is varied. This is because the instability with increasing pressure arises from the presence of the $\ddot {R}$ term in (2.65). Substituting $\ddot {R}$ from (2.63) and simplifying leads to a form resembling Mathieu’s equation, $\ddot {q}+(\delta +\eta \cos {\omega _d t})q=0$ , with some extra terms. The stability of this modified equation is highly sensitive to the parameters $\delta$ and $\eta$ , which depend on $R_0$ , $c_1$ and $\omega _d$ . So, the stability curves do not collapse onto a single curve as the undeformed radius or surface shear modulus is varied.

Figure 7.

(a) Variation of natural frequency with mode number for shear modulus $c_1=0.1\,\mathrm{N\,m}^{-1}$ and (b) critical pressure versus driving frequency curve for shear modulus $c_1=0.1\,\mathrm{N\,m}^{-1}$ in the $\epsilon _p\,\,\textrm {vs}\,\,\omega _d^{\prime}$ plane shown for different values of initial radius.

Figure 7. Long description

The first line graph displays the variation of natural frequency with mode number for shear modulus, with three data lines representing initial radius values of ten micrometers, fifteen micrometers, and twenty micrometers. The x axis represents the mode number ranging from zero to thirty, while the y axis represents the natural frequency scaled by a factor involving density, initial radius, and other constants. The second line graph shows the critical pressure versus driving frequency curve for shear modulus, with three data lines for the same initial radius values. The x axis represents the driving frequency divided by two pi, ranging from zero to eight, while the y axis represents the critical pressure. All values are approximated.

Figure 8.

Plots for an undeformed bubble radius $R_0=10\,\unicode{x03BC} \textrm{m}$ with different values of surface encapsulation shear modulus $c_1$ . (a) Variation of natural frequency with mode number. (b) Critical excitation pressure ( $\epsilon _p$ ) versus driving frequency ( $\omega _d$ ) curves representing the stability boundaries.

Figure 8. Long description

The first line graph shows the variation of natural frequency with mode number for an undeformed bubble radius with different values of surface encapsulation shear modulus. The x-axis represents the mode number ranging from 0 to 30, and the y-axis represents the natural frequency scaled by a factor involving fluid density, initial bubble radius, and shear modulus. Four data series are plotted, each corresponding to different shear modulus values: 0.05 newton per meter, 0.1 newton per meter, 0.15 newton per meter, and 0.2 newton per meter. The second line graph shows the critical excitation pressure versus driving frequency curves representing the stability boundaries for the same conditions. The x-axis represents the driving frequency scaled by a factor involving 2 pi, and the y-axis represents the critical excitation pressure. Four data series are plotted, each corresponding to different shear modulus values: 0.05 newton per meter, 0.1 newton per meter, 0.15 newton per meter, and 0.2 newton per meter. All values are approximated.

5.2.2. Effect of time-varying current in the coil on stability

Time-varying current and pressure are employed as

(5.5) \begin{align} p_\infty =p_0(1+\epsilon _p\cos {\omega _1t}),\quad I=I_0\cos {\omega _2t}. \end{align}

Here, $I_0$ denotes the dimensional current. When only a time-varying current is applied, with $\epsilon _p = 0$ and $\omega _2 = 2\pi f_2$ , the resulting current–frequency stability diagram is obtained for various radii, as shown in figure 9. In this study the magnetic field of range $0.5-1.5\,\textrm {T}$ is used as it is the range used in medical applications and the initial radii of the bubbles considered are $10-20\,\unicode{x03BC} \textrm {m}$ , from this the maximum current amplitude lies in the range $133-266\,\textrm {kA}$ . From the plots, it can be seen that instability occurs mainly at higher current amplitudes, except at a few frequencies. At higher frequencies, no dips are observed, which appear only at lower frequencies. Since the magnetic force is proportional to $H_2R_0^2$ , increasing the radius reduces the critical current amplitude.

Figure 9.

Critical current amplitude ( $I_0$ ) versus frequency $(f_2)$ with $\epsilon _p=0$ , $c_1=0.1\,\mathrm{N\,m}^{-1}$ for different values of initial radius.

Figure 9. Long description

The line graph presents the relationship between critical current amplitude and frequency for different initial radii. The x-axis represents frequency in Hertz, ranging from 1.0 to 4.0 times 10 to the power of 5. The y-axis represents critical current amplitude in Amperes, ranging from 0 to 7 times 10 to the power of 5. Three data lines are shown: one for an initial radius of 12 micrometers in red, one for 15 micrometers in blue, and one for 20 micrometers in black. Each line shows how the critical current amplitude varies with frequency for the respective initial radius. All values are approximated.

Figure 10.

Stability diagram for time-varying current. (a) Critical pressure amplitude versus frequency curves in $\epsilon _p$ vs $\omega ^{\prime}_d$ plane for $f_2=0.5\,\textrm {MHz}$ , initial radius $R_0=10\,\unicode{x03BC} \textrm {m}$ and $c_1=0.1\,\mathrm{N\,m}^{-1}$ . (b) Critical pressure amplitude versus frequency curves in $\epsilon _p$ vs $f_d$ plane for $I_0=150\,\textrm {kA}$ , initial radius $R_0=10\,\unicode{x03BC}\textrm{m}$ and $c_1=0.1\,\textrm {N m}^{-1}$ for different frequencies of current $f_2$ .

Figure 10. Long description

The image contains two line graphs. The first graph (a) shows critical pressure amplitude versus frequency curves for initial currents of fifty kiloamperes and one hundred fifty kiloamperes. The x-axis represents the frequency in megahertz, ranging from zero to three, while the y-axis represents the critical pressure amplitude, ranging from zero to fifteen. The second graph (b) shows critical pressure amplitude versus frequency curves for different frequencies, specifically zero megahertz and point five megahertz. The x-axis and y-axis are labeled similarly to the first graph. All values are approximated.

From the pressure–frequency stability diagram shown in figure 10, it can be observed that the frequency and amplitude of the applied current have little effect on the stability plot, with noticeable variation appearing only at higher pressure frequencies $f_d$ . It is also observed from other simulations that, for other frequencies and amplitudes of the current, the critical pressure curves lie within these two bounds. Interpreting the simulation results becomes challenging, since the governing equations now take the form of a dual-frequency Mathieu equation

(5.6) \begin{equation} \ddot {x}+(\delta +\epsilon _1\cos {\omega _1t}+\epsilon _2\cos {\omega _2t})x=0, \end{equation}

with nonlinear coefficients. A separate, detailed analysis using perturbation methods is needed to gain deeper insights into bubble behaviour.

Another way of generating a time-varying magnetic field is to use oscillating coils with distance $a$ (shown in figure 2) varying as $a=\tilde {a}(1+\beta \cos {2\pi f_3 t})$ . The amplitude of coil oscillation is assumed to be small, with $\beta \ll 1.$ However, from the simulations, it is observed that neither $\beta$ nor $f_3$ influences the pressure–frequency stability plots (not presented in the results section). On the other hand, the periodicity of the modal oscillations depends on the relation between $f_1$ and $f_3$ , which is also observed in the case of a time-varying current. Furthermore, the analysis presented here pertains to specific sources, considering that a more general applied magnetic field could offer a broader understanding of the magnetic bubble dynamics.

5.2.3. Stability diagrams for dipole field

Figure 11.

Critical pressure vs natural frequency curves for a bubble of shear modulus $c_1=0.1\,\textrm {Nm}^{-1}$ , radius $R_0=10\,\unicode{x03BC} \textrm {m}$ for various values of (a) dipole strength at $z=0.01\,\textrm {m}$ and (b) distance of dipole from centre of bubble at $M=1\,\textrm {Am}^2$ .

Figure 11. Long description

The image contains two line graphs labeled (a) and (b). Graph (a) shows critical pressure versus natural frequency curves for a bubble of shear modulus and radius for various values of dipole strength at a specific distance. The graph includes three data lines representing different dipole strengths: one ampere squared, five amperes squared, and ten amperes squared. Graph (b) shows critical pressure versus natural frequency curves for a bubble of shear modulus and radius for various distances of the dipole from the center of the bubble at a specific dipole strength. The graph includes two data lines representing different distances: zero point zero one meters and zero point zero zero six meters. The x-axis represents the natural frequency divided by two pi, and the y-axis represents the critical pressure. All values are approximated.

The computation is performed for bubble of radius $R_0=10\,\unicode{x03BC} \textrm {m}$ , shear modulus $c_1=0.1\,\textrm {Nm}^{-1}$ and susceptibility $\chi =5$ . From figure 11 it is observed that stability decreases as dipole strength $M$ is increased and the distance of the dipole from the centre of the bubble is decreased. The non-dimensional term for magnetic forcing is

(5.7) \begin{equation} \mathcal{M}_n=\frac {12\mu _0\chi M^2R_0^2}{z^8c_1}\sim \mathcal{O}(10^{-4}{-}10^{-2}). \end{equation}

The above order of magnitude is for $z=0.01\,\textrm {m}$ . Although the order of magnitude of forcing is the same as that of the coil field, the dipole field slightly decreases $\delta$ compared with that of the coil field in Mathieu’s equation $\ddot {q}+(\delta +\eta \cos {\omega _dt})q=0$ , which makes it enter the unstable zone at a lower pressure amplitude, decreasing the stability. Here also, at lower values of dipole strength, the stability diagrams remain nearly the same. Decreasing $z$ by half increases the order of the forcing 200 times and significantly reduces the stability zone.

5.3. Time-series analysis of mode shapes

Radial and shape oscillations are analysed only within the stable region of the $\epsilon _p-\omega _d$ plane. In this stable zone, the dynamics is predominantly governed by $R(t)$ and $a_2(t)$ . The behaviour of these two modes is studied as material properties vary. To investigate the influence of different material properties on the time series of radial and shape modes, the governing equations are non-dimensionalised. This is done by introducing the transformations

(5.8) \begin{align} R(t)=R_0\bar {R}(t),\,a_k=R_0\bar {a}_k(t), \,t=\sqrt {\frac {\rho _fR_0^3}{2c_1}}\,\bar {t}, \end{align}

where $c_1=\epsilon \,C_1$ and $c_2=\epsilon \,C_2$ . These transformations lead to the introduction of key non-dimensional quantities: the Reynolds number $Re=\sqrt {{2\rho _fc_1R_0}/{\eta ^2}}$ , which compares the elastic and viscous forces in a liquid, the non-dimensional magnetic force parameter $Rn=({\mu _o\chi _s H_2^2R_0^2}/{2c_1a^2})$ for a coil field and $\mathcal{M}_n$ defined in (5.7) for a dipole field and the strain stiffening parameter $\alpha ={c_2}/{c_1}$ . These non-dimensional parameters facilitate a systematic exploration of how material properties influence the dynamics of the system.

5.3.1. Variation with magnetic susceptibility for coil field

Simulations are conducted with the following parameters: $R_0=10\,\unicode{x03BC} \textrm{m}$ , $c_1=0.1\,\mathrm{N\,m}^{-1}$ , $c_2=0.12\,\mathrm{N\,m}^{-1}$ , $N_1=1000$ turns, $\epsilon =20\, \textrm {nm}$ , $\omega _d=1\,\textrm {MHz}$ , $\epsilon _p=0.5$ and $I=100\,\textrm {kA}$ over a duration of $2000\,\unicode{x03BC} \textrm {s}$ . The time-dependent behaviour of the radial mode and the $a_2$ shape mode is analysed. The results show that the radial mode remains largely unaffected, which can be explained by the relative magnitudes of the non-dimensional forces involved. This is because the orders of magnitude of the various terms are $R_0\sim \mathcal{O}(10^{-5})$ , $c_1\sim \mathcal{O}(0.1)$ , $\chi _s\sim \mathcal{O}(10^{-8})$ , $H_2\sim \mathcal{O}({IN_1}/{a})\sim \mathcal{O}(10^{9})$ and $p_\infty \sim \mathcal{O}(10^5)$ , which gives the order of magnitude of the pressure force and magnetic force as

(5.9) \begin{align} P&=\frac {p_0 R_0}{c_1}\sim \mathcal{O}(10^1),\nonumber\\ Rn&=\frac {\mu _o\chi _s H_2^2R_0^2}{2c_1a^2}\sim \mathcal{O}(10^{-2}-10^{-3}). \end{align}

Therefore, the magnetic forces are significantly smaller compared with the applied pressure, rendering their influence negligible on increasing $\chi$ . However, for the $a_2$ shape mode, the magnetic force serves as the primary driving term, significantly influencing $a_2$ maximum oscillation amplitude. This dependence on susceptibility is clearly observed, as illustrated in figure 12.

5.3.2. Variation with undeformed radius for coil field

Simulations are conducted with $c_1=0.1\,\mathrm{N\,m}^{-1}$ , $c_2=0.12\,\mathrm{N\,m}^{-1}$ , $\chi =1$ , $\epsilon =20\,\textrm {nm}$ , $\omega _d=1\,\textrm {MHz}$ , $\epsilon _p=0.5$ and $I=100\,\textrm {kA}$ over a duration of $100\,\unicode{x03BC} \textrm{s}$ . From the expressions for the non-dimensional magnetic force, it is evident that the forcing on the shape modes increases with an increase in radius. Additionally, as shown in (2.63), the pressure acting on the radial mode opposes the radial motion, which leads to a reduction in the amplitude of radial oscillations with increasing radius. This behaviour is clearly depicted in figure 13. Also, the pressure force is proportional to $R_0$ and the magnetic forcing is proportional to $R_0^2$ . This scaling is reflected in the $a_2$ response, where its amplitude increases, respectively, by approximately $2.25$ and $4$ times as $R_0$ varies from 10 to $20\,\unicode{x03BC}$ m.

Figure 12.

Variation of radial mode amplitude $R(t)$ and second mode amplitude $a_2(t)$ with time for a bubble of initial radius $R_0=10\,\unicode{x03BC} \textrm {m}$ at various values of encapsulation susceptibility. The simulations are conducted at $\omega _d=1\,\textrm {MHz}$ and $\epsilon _p=0.5$ .

Figure 12. Long description

The line graph consists of four subplots. The top left subplot shows the variation of radial mode amplitude with time, with the x-axis labeled as time in seconds and the y-axis labeled as R(t)/R0. The top right subplot shows the variation of second mode amplitude with time for a value of chi equal to 0.5, with the x-axis labeled as time in seconds and the y-axis labeled as a2(t)/R0. The bottom left subplot shows the variation of second mode amplitude with time for a value of chi equal to 1, with the x-axis labeled as time in seconds and the y-axis labeled as a2(t)/R0. The bottom right subplot shows the variation of second mode amplitude with time for a value of chi equal to 1.5, with the x-axis labeled as time in seconds and the y-axis labeled as a2(t)/R0. All values are approximated.

Figure 13.

Variation of radial mode amplitude $R(t)$ and second mode amplitude $a_2(t)$ with time for a bubble encapsulation of susceptibility $\chi =1$ at different values of undeformed bubble radius $R_0$ . The simulations are conducted at $\omega _d=1\,\textrm {MHz}$ and $\epsilon _p=0.5$ .

Figure 13. Long description

The image contains four line graphs that illustrate the variation of radial mode amplitude and second mode amplitude with time for a bubble encapsulation of susceptibility at different values of undeformed bubble radius. The top left graph shows the radial mode amplitude over time for a bubble radius of 10 micrometers, with the blue line representing the data. The top right graph displays the second mode amplitude over time for the same bubble radius, also with a blue line. The bottom left graph presents the second mode amplitude over time for a bubble radius of 15 micrometers, depicted by a red line. The bottom right graph shows the second mode amplitude over time for a bubble radius of 20 micrometers, represented by a green line. Each graph has time on the x-axis and the respective amplitude ratio on the y-axis. The simulations are conducted at specific conditions, and the graphs show how the amplitudes vary over a short time interval from 1.995 to 2.000 milliseconds. The data indicates oscillatory behavior in the amplitudes, with different frequencies and magnitudes depending on the bubble radius.

Figure 14.

Variation of second mode amplitude $a_2(t)$ with time for a bubble encapsulation of susceptibility $\chi =1$ , initial radius $R_0=10\,\unicode{x03BC} \textrm{m}$ , current amplitude $I_0=150\,\textrm {kA}$ , pressure amplitude $\epsilon _p=0.5$ and frequency of pressure $f_1=1\,\textrm {MHz}$ for different frequencies of current $f_2$ .

Figure 14. Long description

The image contains four line graphs, each representing the variation of second mode amplitude with time for a bubble encapsulation of susceptibility. The graphs are arranged in a 2x2 grid. Each graph shows a different frequency of current. The x-axis represents time in seconds, ranging from 1.990 to 2.000 times 10^-3 seconds. The y-axis represents the normalized amplitude a2/R0, with different scales for each graph. The blue line graph in the top left shows a frequency of 0.5 times the initial frequency. The red line graph in the top right shows a frequency of 0.75 times the initial frequency. The green line graph in the bottom left shows a frequency equal to the initial frequency. The purple line graph in the bottom right shows a frequency of 0.5 times the initial frequency. Each graph displays oscillatory behavior with varying amplitudes and frequencies. All values are approximated.

5.3.3. Time-series analysis of modes with applied frequency of current in the coil

Figure 15.

Variation of second mode amplitude $a_2(t)$ and phase plot with time for a bubble encapsulation of susceptibility $\chi =1$ , initial radius $R_0=10\,\unicode{x03BC} \textrm{m}$ , current amplitude $I_0=150\,\textrm {kA}$ , pressure amplitude $\epsilon _p=0.5$ and frequency of pressure $f_1=1\,\textrm {MHz}$ and $f_2=1.777f_1$ .

Figure 15. Long description

The line graph on the left shows the variation of the second mode amplitude over time, with the x-axis representing time in seconds and the y-axis representing the amplitude ratio. The amplitude fluctuates between approximately 0.275 and 0.290. The phase plot on the right illustrates the relationship between the amplitude ratio and the phase, with the x-axis representing the amplitude ratio and the y-axis representing the phase. The phase plot forms a complex, intertwined pattern. All values are approximated.

It can be seen from figure 14 that, as the frequency of the applied current is increased, proportional to the frequency of the pressure field, the periodicity of the oscillation pattern of $a_2$ mode continuously varies with pressure and current varying as

(5.10) \begin{align} p_\infty =p_0(1+\epsilon _p\cos {\omega _1t}), \quad I=I_0\cos {\omega _2t}. \end{align}

Here, $\omega _1=2\pi f_1$ is the angular frequency of the pressure field and $\omega _2=2\pi f_2$ is the angular frequency of the applied current. The pattern of only the $a_2$ mode is studied, as it is the dominant mode within the stable zone. The frequency of the current is found to have only a slight influence on the amplitude of the $a_2$ mode, with $f_2 = 0.5f_1$ producing the maximum effect. This occurs because $\cos ^2{\omega _2 t}$ expands to $1 + \cos {2\omega _2 t}$ , which directly adds to the pressure force, thereby increasing the amplitude. Further simulations reveal that the oscillations are predominantly chaotic, except when $f_2$ is a fractional multiple of $f_1$ , as shown in figure 15. This observation is noteworthy, since operating at frequencies that induce a strong chaotic behaviour makes controlling the bubble dynamics increasingly challenging. A more comprehensive analysis of bifurcation patterns under varying material parameters may be conducted separately to gain deeper insights into the oscillatory characteristics of magnetic bubbles.

5.3.4. Time-series variation with dipole strength for dipole field

Simulations are conducted with the following parameters: $R_0=10\,\unicode{x03BC} \textrm{m}$ , $c_1=0.1\,\mathrm{N\,m}^{-1}$ , $c_2=0.12\,\mathrm{N\,m}^{-1}$ , $\omega _d=1\,\textrm {MHz}$ , $\chi =5$ , $z=0.01\,\textrm {m}$ and $\epsilon _p=0.5$ over a duration of $2000\,\unicode{x03BC} \textrm{s}$ . It can be observed from figure 16 that increasing the dipole strength leads to a slight increase in the amplitude of radial oscillations, while significantly enhancing the magnitude of the second mode oscillation. From the expression for $\mathcal{M}_n$ in (5.7), it is evident that the magnetic forcing scales with $M^2$ . Consequently, increasing the dipole strength by a factor of ten amplifies the magnetic forcing by a factor of one hundred, leading to a noticeable impact on radial oscillations and a pronounced effect on second mode.

Figure 16.

Time-series response of radial and second mode for a bubble of radius $R_0=10\,\unicode{x03BC} \textrm{m}$ , surface shear modulus $c_1=0.1\,\textrm {Nm}^{-1},$ susceptibility $\chi =5$ for different values of dipole strength $M$ .

Figure 16. Long description

The image contains four line graphs showing the time-series response of radial and second mode for a bubble of radius surface shear modulus susceptibility for different values of dipole strength. The top left graph shows the response for M equals 1 ampere squared, the top right graph for M equals 5 ampere squared, and the bottom right graph for M equals 10 ampere squared. The x-axis represents time in seconds, and the y-axis represents the normalized response. Each graph displays oscillatory behavior with varying amplitudes and frequencies. The top left graph shows a higher frequency oscillation with smaller amplitude, while the bottom right graph shows a lower frequency oscillation with larger amplitude. The top right graph shows a moderate frequency oscillation with medium amplitude. The bottom left graph shows the radial second mode response for M equals 5 ampere squared. All values are approximated.