A nonlinear elasticity approach to modelling the collapse of a shelled microbubble

There is considerable interest in using shelled microbubbles as a transportation mechanism for localized drug delivery, specifically in the treatment of various cancers. In this paper a theoretical model is proposed which predicts the dynamics of an oscillating shelled microbubble. A neo-Hookean, compressible strain energy density function is used to model the potential energy per unit volume of the shell. The shell is stressed by applying a series of small radially directed stress steps to the inner surface of the shell whilst the outer surface is traction free. Once a certain radial deformation is reached, the stress load at the inner radius is switched off causing the shell to collapse and oscillate about its equilibrium (stress free) position. The inflated shell configuration is used as an initial condition to model the time evolving collapse phase of the shell. The collapse phase is modelled by applying the momentum balance law and mass conservation. The dynamical model which results is then used to predict the collapse time of the shelled microbubble as it oscillates about its equilibrium position. A linear approximation is used in order to gain analytical insight into both the quasistatic inflationary and the oscillating phases of the shelled microbubble. Results from the linearized model are then analysed which show the influence of the shell’s thickness, Poisson ratio and shear modulus on the rate of oscillation of the shelled microbubble. The nonlinear model for the quasistatic state is solved numerically and compared to the linearized quasistatic solution. At present, there is no solution to the nonlinear collapsed state. This is a future area of research for the current authors.


Introduction
Premanufactured shelled microbubbles are currently licensed in the UK as ultrasound imaging contrast agents (Narayan & Wheatley, 1999).Current research is focussing on using the microbubbles as a transportation mechanism for localized drug delivery specifically in the treatment of various cancers (Peregrino et al., 2013;Gourevich et al., 2015;Escoffre et al., 2013;Fan et al., 2013a;Yan et al., 2013;Niu et al., 2013).Ultrasound contrast agents (UCAs) are shelled microbubbles typically composed of a layer or several layers of a protein shell encapsulating a perfluoro gas that helps to stabilize the shelled microbubble when it is injected into the bloodstream (Stride & Saffari, 2003;Cosgrove, 2006;Stride & Coussios, 2010).The shelled microbubbles have a typical radius of between 1 and 4 μm allowing them to propagate through the capillaries in the human body and a shell thickness that varies between 4 and 100 nm depending on whether the UCA is a monolipid or polymer variant (McLaughlan et al., 2013).A typical shear modulus value for a monolipid UCA is 20 MPa with a Poisson ratio of ν = 0.48 (Falou et al., 2010;Sotiriou et al., 2014a).UCA's resonate with typical frequencies in the range of 1 to 10 MHz producing nonlinear, multiple harmonic signals that enhance the quality of the medical imaging process (Blomley et al., 2001).There has been a research momentum growing in recent years to use the UCAs as localized drug delivery agents (Lentacker et al., 2009).This procedure aims to minimize the pernicious side effects associated with current conventional chemotherapy treatments.It is worth emphasising that there are competitor technologies being proposed and, for example, some studies have focussed on exploiting the photothermal properties of gold and silver coated nanobubbles to kill the cancer cells Sotiriou et al. (2014b).Much progress has been made but much remains to be done before this can be deployed routinely in patients Dewitte et al. (2015).Hence there is a need to develop virtual simulation tools to better understand the challenges Fan et al. (2013b).This paper contributes to this effort by identifying how the material parameters of the shell influence the dynamics and collapse time of the shelled microbubble.In this paper we define the collapse time to be the time taken for a stressed shelled microbubble to collapse back to its stress free state.
Most current shelled microbubble models are based on the Rayleigh-Plesset equation for a free gas microbubble, which is derived by applying pressure balances to the inner surface of the shelled microbubble with those acting on the outside of the shelled microbubble's surface and the surrounding liquid.The Rayleigh-Plesset equation assumes that the shelled microbubble oscillations are purely radial and that the surrounding liquid is incompressible Rayleigh (1917); Plesset (1949).The gas in the shelled microbubble is assumed to behave adiabatically despite its polytropic index being relatively close to one which is associated with isothermal behaviour Doinikov & Bouakaz (2011).The damping contributors are the liquid viscosity that surrounds the shelled microbubble, thermal damping between the gas and the surrounding liquid, and damping associated with liquid compressibility via the external acoustic energy.The thermal damping can be accounted for by selecting the appropriate value for the polytropic index (Doinikov & Bouakaz, 2011).The Rayleigh-Plesset equation can be modified to take account of the compressible nature of the surrounding liquid (Doinikov & Bouakaz, 2011).Marmottant et al. (2005) modelled lipid shelled microbubbles using a Rayleigh-Plesset equation; this assumes that when the shelled microbubble oscillates, the elastic region holds only for a small range of radii.However, this model is capable of describing nonlinear effects particularly the compression only behaviour observed in the analysis of certain monolipid shelled microbubbles (Marmottant et al., 2005).This work was the first to propose a model for the surface tension of the shell which varies during various stages of its oscillatory motion.They achieved this by expressing the surface tension as a piecewise function dependent on the shell's area density with the shell experiencing a smaller surface tension during contraction.This is due to the area per molecule decreasing during contraction resulting in a smaller number of molecular interactions and therefore a smaller surface tension.Paul et al. also used empirical means to develop a Rayleigh-Plesset model using an effective surface tension that incorporated a dilatational elasticity term which was a function of the surface area change (Paul et al., 2010).The dilatational elasticity term was modelled using both quadratic and exponential terms with each approach giving similar results and predicting the subharmonic response with a reasonable level of accuracy.Both of these approaches are empirical and rely heavily on experimental observations for one particular type of UCA.Church developed a model for an encapsulated shelled microbubble incorporating both the inner and outer radii and the interfacial surface tension (Church, 1995).Church assumed that the incompressible encapsulated gas shelled microbubble was surrounded by an unbounded incompressible liquid and experienced radial oscillations when it was subjected to an external acoustic pressure.Church's model uses two surface tension expressions; a surface tension for the outer shell/liquid interface and an interfacial surface tension between the shell's inner radius and the gas layer.Church exploited the Rayleigh-Plesset equation and assumed that the shell material behaved as a viscoelastic solid obeying the Kelvin-Voigt constitutive equation (Church, 1995).Experimental evidence indicates that the Church model is more suited to Albumin shelled microbubbles with shell thicknesses of the order of 15-20 nm (Doinikov & Bouakaz, 2011).Doinikov & Dayton (2007) suggested a model for lipid encapsulated microbubbles where the shell is treated as a viscoelastic fluid that obeys the Maxwell constitutive equation.The Maxwell model developed by them helped to explain some of the experimental data that contradicted previous assumptions that the lipid shell behaves as a viscoelastic solid.The Maxwell model predicts that the resonance frequency of a shelled microbubble possessing a viscoelastic fluid shell can have both higher and lower resonance frequencies than that of a free microbubble.This depends on the choice of the shell parameters (Doinikov & Bouakaz, 2011;Doinikov & Dayton, 2007).Doinikov et al. (2009) proposed a modification to their model where the shell's viscous constant was replaced by a function of the shell's shear rate.They introduced this approach to attempt to solve the dependency of the shell viscosity on the initial shelled microbubble's radius (Doinikov & Dayton, 2007).Despite numerous Rayleigh-Plesset models existing for UCAs, experimental observations have been made that challenge the current, existing models.The experimental observations of compression only behaviour for monolipid coated UCAs is highlighted by the asymmetrical oscillation curves that are experimentally observed and that were subsequently modelled empirically by Marmottant et al. (2005) who separated the shell elastic behaviour into three different regions.According to current models, the shell's material parameters such as viscosity and elasticity, display a dependency on the initial radius of the shelled microbubble.However, large shelled microbubbles possess a greater mass and surface area yet should still have the same viscosity and elasticity for the same type of material.This highlights a current flaw in the modelling of the rheological properties of monolipid UCAs which clearly requires a more fundamental and mathematically rigorous treatment (Doinikov & Bouakaz, 2011).
There currently exists very little literature pertaining specifically to UCA modelling using nonlinear elasticity; which is a common approach for modelling large deformations of elastic materials and in particular soft materials such as in biological settings (Gorb & Walton, 2010;Fung, 1990Fung, , 1993Fung, , 1996;;Holzapfel, 2000;Holzapfel & Ogden, 2004).There are however, numerous publications related to the dynamics of spherical bodies using nonlinear elasticity (Lazopoulos & Ogden, 1999).A recent paper by Efthymium and Tsiglifis uses constitutive laws from nonlinear elasticity alongside the Kelvin-Voigt viscoelastic model to study the physical behaviour of various UCA types ranging from monolipids to polymers (Efthymiou et al., 2012).They reported that the polymer based UCAs were consistent with the neo-Hookean model whereas monolipid UCAs were consistent with the Mooney-Rivlin constitutive law due to the presence of strain softening behaviour.Strain softening behaviour occurs due to the area density of the monolipid decreasing as the material stretches radially outwards.This behaviour has been observed in monolipids typically used in UCA shells such as Sonovue Doinikov & Bouakaz (2011) and Tsiglifis & Pelekasis (2008).
This article will derive a new nonlinear elasticity model of a UCA to explore the role that the shell material properties and stress have on the UCA's dynamics.This article takes a novel approach to modelling shelled microbubbles and is quite distinct from the previous literature such as Marmottant et al. (2005) and Doinikov & Dayton (2007).A previous study by Moulton and Goriely has highlighted how difficult it is to model the full collapse of an elastic shell (Moulton & Goriely, 2011).This study focusses on calculating a collapse rate rather than trying to totally collapse a shelled microbubble.No previous study has focussed on the influence of the material parameters on the collapse time of the shelled microbubble specifically in relation to the shell's Poisson ratio and its shear modulus.
A neo-Hookean compressible strain energy density function is used to model the stress profiles within the thin shell.An initial series of small stress steps are applied symmetrically at the inner radius of the shell whilst the outer surface is traction free.Note that the applied stress at the inner radius is not attempting to model the gas in the shelled microbubble.In order to study the collapse of the shelled microbubble it has to be initiated with some non-equilibrium stress distribution.One way to achieve this is to inflate a shelled microbubble from a stress free reference configuration.The quasistatic phase describes this perturbation from its equilibrium state to a stressed state.This spherically symmetric deformation is quasistatic in nature and is thus considered to be independent of time.The spatial profiles for the Cauchy radial and hoop stresses within the shell will be used along with the conservation of mass and energy to model the subsequent collapse of the shelled microbubble.In order to collapse the shell, the applied stress at the inner radius is instantaneously switched off whilst the outer surface is traction free.This causes the shell to collapse and oscillate about its stress free configuration.This leads to a simple harmonic motion (SHM) description of the shelled microbubble dynamics where the collapse time can be determined.Note that the collapse time is defined in this paper to be the time taken for the collapsing shelled microbubble to first return to its reference (equilibrium) state.This is one quarter of the period of oscillation of the shelled microbubble with the SHM motion arising due to a lack of damping in the physical model.The nonlinear model for the quasistatic phase is solved numerically and compared and contrasted with the linearized analytical solution.This comparison is of course conducted over a small range of displacements, since this is what the linearized version is predicated upon, but future work will explore the nonlinear behaviour of the model derived in this paper.However, the collapse phase is solved for a linearized model and not the nonlinear model.The numerical solution of the nonlinear collapse phase is considered beyond the scope of this article.The material parameters of the shell are then adjusted to determine their influence on the time taken for the stressed shell to collapse back to its stress free equilibrium position.It is important to realize that the nature of the strain energy density function is only relevant to the nonlinear model for the quasistatic phase.

Modelling the quasistatic stressing of a shelled microbubble
Let us start by considering the reference configuration of a stress free, fully formed spherical shell with an inner and outer radii of R I and R O , respectively.The current configuration is then an inflated spherical shell which will possess a stress.Let the current configuration basis vectors be represented by e r , e θ and e φ .To distinguish between the reference and current vector fields, we employ an uppercase (G) for the reference configuration and a lowercase (g) for the current configuration (Ogden, 1997).The stress is due to a stress load that is applied to the interior surface of the shelled microbubble shell.A mixed tensorial basis is used where the radial deformation is defined as χ = χ i g i and the gradient of the deformation F is given by Ogden (1997) In spherical polar coordinates, the current configuration can be transformed into physical components (Ogden (1997)) to give χ 1 = χ r , χ 2 = rχ θ , and χ 3 = r sin θχ φ .Using equation (1), we can determine the gradient of the deformation defined by χ 1 = r (R), and χ 2 = χ 3 = 0. Assuming θ (Θ) and φ (Φ), the nonzero terms in equation ( 1) are, (∇ ⊗χ) rR e r ⊗ e R = r e r ⊗ e R , (∇ ⊗χ) θΘ e θ ⊗ e Θ = (rθ /R)e θ ⊗ e Θ and (∇ ⊗ χ) φΦ e φ ⊗ e Φ = rφ sin θ/(R sin Θ)e φ ⊗ e φ (see Cowley et al., 2015).If we assume that the shell material is hyperelastic then there exists a strain energy density function expressing the potential energy per unit volume.We will use a neo-Hookean strain energy density function (Holzapfel, 2000;Gorb & Walton, 2010;Gou et al., 2012), W (F), which includes a compressible term that is used to model the change in volume of the shelled microbubble as it is inflated.The determinant of F, gives a measure of how the volume of the spherical shell changes as it maps from the stress free, reference configuration to the stressed, current configuration.This determinant is denoted by J and if we further assume that θ = Θ and φ = Φ then ( 2 ) A compressible, neo-Hookean (Gorb & Walton, 2010) strain energy density function is used to model the potential energy per unit volume of the shell.Previous studies by Tsiglifis & Pelekasis (2008) and Efthymiou et al. (2012) show that thick shelled ultrasound contrast agents (thicknesses of the order of 15nm and above) can be modelled by using a neo-Hookean strain energy density function whereas monolipid shells (thicknesses of the order of 4nm) are consistent with a Mooney-Rivlin strain energy density function.This study will focus on shell thicknesses between 20 and 50 nm which will be described via a compressible, neo-Hookean strain energy density function.Using a neo-Hookean strain energy density function has the added advantage of being among the simplest of nonlinear models to study.Nonlinear elasticity can account for large scale deformations and has been used in numerous studies pertaining to modelling rubber mediums and soft tissue (Gorb & Walton, 2010;Daniel et al., 2010).
Experimental studies (Falou et al., 2010;Tu et al., 2009) have shown that the Poisson ratio of a typical ultrasound contrast agent is between 0.48 and 0.49 which is similar to rubber.The common approach to modelling rubber for moderate strains is to utilize nonlinear strain energy density functions which are either neo-Hookean or Mooney-Rivlin (De Pascalis et al., 2013).The key reason for introducing a compressible, third invariant is to account for the material parameter described by Poisson's ratio and to identify how this parameter influences the collapse time of the shell whereas a neo-Hookean or Mooney-Rivlin strain energy density function only utilizes the shear modulus (De Pascalis et al., 2013).
The compressible neo-Hookean strain energy density function is (Gorb & Walton, 2010) where FF T is the left Cauchy-Green deformation tensor (Holzapfel, 2000), μ is the shear modulus, ν is Poisson's ratio and β = ν/(1 − 2ν).The stresses can be described using the first Piola-Kirchhoff stress tensor, relating the force in the current configuration to the area in the reference configuration (Gorb & Walton, 2010).On the other hand, the Cauchy stresses relate the force in the current configuration to the area in the current configuration.The first Piola-Kirchhoff stress tensor, S(F), is then Substituting equation (1) into equation ( 4) gives For the quasistatic solution of a stressed sphere the divergence of the first Piola-Kirchhoff stress tensor will be equal to zero.This is because the quasistatic deformation is independent of time which results in the Cauchy momentum equating to zero.To determine ∇ •S = 0, we need to be able to relate the physical coordinates for the mixed tensorial basis to the general basis vectors represented by the components g i and G i where i ∈ {1, 2, 3}.The first Piola-Kirchhoff stress tensor is represented by Ogden (1997), S = S j i g i ⊗ G j where S j i are the left-covariant components of S. Converting into physical coordinates gives and Calculating the divergence of S gives the following radial and angular equations and, If we further assume that the deformation is purely radial then θ = Θ, resulting in S 2 2 = S 3 3 .Equation ( 10) is the same as the model in Daniel et al. (2010) and the relationship S 2 2 = S 3 3 is because we have a spherically symmetric, purely radial deformation.This implies that both the polar and azimuthal stresses have the same dependency on the radial deformation.To formulate the Cauchy stresses, the radial differential equation (10) has to be rewritten in terms of the physical coordinates.Using equation ( 10) and rearranging, gives To determine the Cauchy radial and hoop stresses, equation ( 12) is solved for the inner surface of the shell being subjected to a load stress, p, and the outer shell surface is traction free.At R = R I (the inner radius), this inner boundary condition is S rR e R = −pJF −T • e R = −pr 2 /R 2 e R , and so, from equations ( 2) and ( 7) Calculating the boundary condition for the outer shell's radius R = R O gives, S rR = 0, and so . ( 14) Equation ( 12) can now be solved subject to the boundary conditions given by equations ( 13) and ( 14).The Cauchy stresses represented by τ Holzapfel (2000) are then obtained from the first Piola-Kirchhoff stress tensor via τ = SF T /J.The radial and hoop stresses are then evaluated using equation ( 6), to give and, τ φφ = μ

Nondimensionalization and linearization of the inflationary model
Nondimensionalization is used to assist in solving the quasistatic and collapse phases of the shelled microbubble.This is achieved by using the reference configuration's inner radius, R I .This leads to Nondimensionalizing the boundary condition at the inner radius given by equation ( 13) and rearranging gives where p = p/μ.Similarly, nondimensionalizing the boundary condition at the outer radius represented by equation ( 14) leads to Nondimensionalising the Cauchy stresses given by equations ( 15)-( 17) results in and, τφφ = To afford some analytical insight, equation ( 18) is linearized using where 0 < p << 1. Equation ( 24) assumes that the displacement y(Y ) in the current configuration for the quasistatic phase, is perturbed by a very small amount which is denoted by pf (Y ) where p is a nondimensionalized parameter and is expressed in terms of the applied stress load p and the shear modulus μ. f (Y ) denotes a function of Y , the independent variable in the reference configuration, and is related to the difference between the displacement in the current and reference configurations.The use of this is justified if all the stress components developed throughout the body are small in comparison with the shear modulus of the shell.Typical imaging pressures are of the order of p = 200 kPa whereas the shear modulus of a monolipid is μ = 20MPa resulting in p/μ = 0.01 (Marmottant et al., 2005).Substituting equation ( 24) into equation ( 18) results in the leading order equation This can be solved to give where f (Y ) = A/Y 2 + BY and A and B are constants determined using the boundary conditions represented by equations ( 13) and ( 14).The linearized boundary conditions give expressions for the parameters A and B, and the radial deformation y(Y ) is then approximated by Using this solution and linearising the expressions for the Cauchy stresses gives the approximations and

Linearization of the time evolving collapse phase of the shell
Applying the momentum balance law where ρ o denotes the density in the reference configuration, v represents the velocity and t denotes the time, gives where the radial component of the material derivative is given by Dv r /Dt = ∂v r /∂t (Ogden, 1997, pp. 143-145).Writing equation ( 30) in terms of the Cauchy stress leads to where v = v r e r , v φ = 0 and v θ = 0 (radial dependency only) (Acheson, 1990, pp. 354-355).To collapse the shell a change in the boundary conditions has to be applied at the inner radius of the shell.
In the inflationary picture there is a stress applied at the inner radius, directed radially outwards, but to collapse the shell the stress at the inner radius is traction free.The right hand side of equation ( 31) is given by where τ rr and τ θθ are given by equations ( 15) and ( 16).For v r = ∂r/∂t, the term r denotes the position of a material particle which can be written as r = r(t; R).Using the chain rule Equation ( 31) can be rewritten as To nondimensionalize time we use t = γ t and set γ = ρ o R 2 I /μ which leads to the nondimensionalized momentum balance Applying the boundary conditions at the inner and outer radii which is traction free such that τ yy (Y and where the nondimensionalized initial conditions are given by the inflated solution and Linearization is performed once again where Linearizing ( 35) gives The initial conditions are given by and with boundary conditions that are stress free at both the inner and outer radius given by Applying a Laplace transform to equation ( 41) and applying the initial conditions leads to where Consider first the case when s = 0.Here equation ( 47) becomes which has solution G(Y , 0) = E/Y 2 + FY , where E, F ∈ R. The Laplace transformed boundary conditions arising from equations ( 44) and ( 45) result in E = F = 0. Hence the solution pertaining to the case s = 0 is the trivial solution G(Y , 0) = 0. From now on we consider the case when s = 0.The homogeneous equation has the solution where i 1 (αsY ) is the modified spherical Bessel function of the first kind of order one and k 1 (αsY ) is the modified spherical Bessel function of the second kind of order one (Abramovitz & Stegun, 1972, pp. 443-445) and (Arfken, 1985, pp. 633-634).The modified spherical Bessel functions can be written in terms of more familiar functions as The particular integral solution is since f (y) satisfies equation ( 25).The general solution is then given by Taking the Laplace transform of the stress free boundary conditions given by equations ( 44) and ( 45) and substituting equation ( 51) into the Laplace transforms of ( 44) and (45) gives We can write equations ( 52) and ( 53) as the matrix equation where Equation ( 54) leads to and where the determinant is given by The physical system becomes singular when the determinant given by equation ( 57) is equal to zero i.e. det(M) = 0.This corresponds to the poles of the integrand in the inverse Laplace transform formulation.
The location of the poles in the complex plane can be empirically observed by producing numerical plots of det(M(s)) where s = x + iω, over a range of material parameter sets.The poles ω j were identified numerically using FindRoot in Mathematica (Wolfram Research Inc., 2016) which uses the Newton method to identify roots.The residuals were also evaluated numerically using NResidue in Mathematica (Wolfram Research Inc., 2016) for a range of β and Y O values.The general form of the solution is given by which simplifies to where A n (Y ) denotes the amplitude of each residual whose pole is given by s = ±iω n where equation ( 59) is a Fourier series.The amplitudes for higher harmonics are very small in comparison to A 1 .Hence the collapse time is dominated by the first angular frequency ω 1 .Every component in expression ( 59) independently satisfies the boundary conditions given by equations ( 44) and ( 45).The collapse time t * , which is the time taken for the stressed shell to collapse to its stress free state, is therefore given by t * = π/(2ω 1 ).

Results for the inflationary phase of the shelled microbubble
By numerically solving the nonlinear inflating shelled microbubble equations given by equations ( 18)-( 20) for different stress loads (p) applied to the inner surface of the shelled microbubble, the effect on the shell deformation can be determined.Equations ( 18)-(20) were solved numerically using an in-built solver in Mathematica using the shooting method (Wolfram Research Inc., 2016). Figure 1 shows the effect on the radial deformation of the shelled microbubble for a series of stress loads with shelled  18) subject to the boundary conditions given by equations ( 19) and ( 20).The linearized model is calculated using equation ( 27).18) and ( 21) subject to the boundary conditions given by equations ( 19) and ( 20).The linearized model is calculated using equation ( 28).
microbubbles of initial thickness Y O − Y I = 0.02. Figure 1 shows that both the nonlinear and linearized models are approximately linear with respect to Y with the nonlinear model displaying a slightly larger radial deformation.For this choice of parameters, the linearized model is valid up to p ≈ 0.003.It transpires that decreasing the shell's thickness from Y O − 1 = 0.2 to Y O − 1 = 0.02 results in a greater radial deformation for the shelled microbubbles.This greater radial deformation is due to thinner shells requiring a greater strain in order to achieve the stresses required to balance with the applied stress load.Note that when no stress load is applied (p = 0) then y(Y ) = Y which is the trivial solution satisfying equations ( 18) to (20) as y(1) = 1, y (Y ) = 1 and y(Y O ) = Y O . Figure 2 illustrates the behaviour of the nondimensionalized Cauchy radial stresses for the spherical inflation of a shelled microbubble of initial thickness Y O − 1 = 0.02.Each curve arises from a specific stress load applied to the inner surface of the shell.The negative values imply that the shell is being compressed as the stress load increases.The compressive nature is a consequence of the shell becoming thinner.It can also be seen that the largest compression occurs at the inner surface of the spherical shell.The gradients of the curves are positive and increase in magnitude as the stress load increases.Note that the Cauchy radial stress at the inner radius is equal in magnitude to the applied stress steps that are applied to the inner radius of the shelled microbubble.Figure 2 shows that both the nonlinear and linearized models are approximately linear with respect to y with the nonlinear model experiencing a slightly larger radial deformation.Figure 3 shows the nondimensionalized Cauchy hoop stress τθθ for an initial shell thickness of Y O − Y I = 0.02; the positive values indicating a stretching behaviour.As the thickness of the shell reduces the size of the Cauchy hoop stress increases along with the radial deformation of the shelled microbubble.The Cauchy hoop stress τ θθ is equal in magnitude to τ φφ due to the spherically directed nature of the radial deformation.There is a significant increase in the magnitude of the Cauchy hoop stresses in comparison to the Cauchy radial stresses, indicating that the hoop stresses play a dominant role in the collapse of the shelled microbubble.The thinning down of the shell normal to the surface due to a compressive stress and the stretching out of the shell tangential to the surface as a result of the positive hoop stresses, result in the material particles that make up the shell at the inner radius experiencing greater tensions.Figure 3  18) and ( 22) subject to the boundary conditions given by equations ( 19) and ( 20).The linearized model is calculated using equation ( 29).
shows that both the nonlinear and linearized models are approximately linear with respect to y with the nonlinear model experiencing a slightly larger radial deformation and hoop stress.

Results for the collapse phase of the shelled microbubble
Having produced results for the inflationary phase of the microbubble's evolution, this can be used as an initial condition to model the collapsing shell.The collapse is achieved by setting the boundary condition on the inner surface of the shell to be stress free.Figure 4 plots g(Y , t) at the nondimensionalized inner and outer radius of a collapsing shell versus the nondimensionalized time for an initial condition based on an inflationary stress load that is 1% of the shear modulus of the shell (p = 0.01).The resulting nonlinear trend is a sinusoidal function with a collapse time, t * ≈ 0.47.There are no dissipative, damping terms in the momentum balance law such as viscoelastic terms connected to the viscosity of the shell and the surrounding fluid.Neglecting the fluid resistance results in a simple harmonic behaviour with no resulting energy dissipation to the surrounding medium.Both the inner and outer radius of the shell experience the same collapse time despite their radial deformations being different with the inner radius experiencing a slightly larger radial deformation.Figure 5 highlights the nonlinear relationship between the redimensionalized collapse time of the shell and the redimensionalized shear modulus.The nondimensionalized angular collapse frequency is evaluated using equation ( 54) where the nondimensionalized collapse time is given by t * = π/(2ω).We redimensionalize the collapse time using t = γ t where γ = ρ o R 2 I /μ. Figure 5 shows that as the shear modulus increases the microbubble shells experience faster collapse times.Figure 6 shows the relationship between the collapse time of the shell and the nondimensionalized original (stress free) thickness of the shell.Thinner shells will strain more to balance the tensions that they are subjected to which will result in larger radial deformations (this is for a fixed initial stress of p = 0.01).The hoop stresses play a key role in the collapse phase of the shell with thinner shells experiencing a larger hoop stress.Hence, despite thin shells experiencing larger radial deformations, they also experience a significant increase in their hoop stresses.(Falou et al., 2010) and its inner radius is R I = 1μm.Typical shear modulus values are within the range of μ = 20 → 100MPa (Falou et al., 2010;Tu et al., 2009).The collapse time is calculated from the zeros of the determinant of the matrix M(s) in equation ( 54) and is redimensionalized using t = γ t where γ = ρ o R 2 I /μ.
A NONLINEAR ELASTICITY APPROACH 797 Fig. 6.Graph of the nondimensionalized collapse time of the shell versus the nondimensionalized original thickness of the shell where β = 12, ν = 0.48, and μ = 1.This is calculated from the zeros of the determinant of the matrix M(s) in equation ( 54).
Thinner shells also possess a smaller mass and it is the combination of the mass and the hoop stress that must be considered when determining the collapse time.Figure 6 illustrates that thicker shells have slightly longer collapse times despite their smaller radial deformations.Figure 7 highlights the approximately linear relationship between the nondimensionalized collapse time and the Poisson ratio for an initial stress of p = 0.01.The Poisson ratio is the negative of the ratio of the transverse to axial strain (Gorb & Walton, 2010) and so shells with smaller Poisson ratios will experience larger radial deformations for a given fixed initial stress.This will result in a longer collapse time for smaller Poisson ratios due to their larger radial deformation.Figure 7 illustrates how the collapse times are lower for a higher Poisson ratio.

Conclusion
This article described an analytical and numerical approach to modelling the inflationary process of a shelled microbubble via a quasistatic radially directed stress load applied to its inner surface.The full nonlinear quasistatic equation was solved numerically using the shooting method and this solution was compared to the linearized model.The stress load was switched off (a traction free boundary condition was applied) and the time for the microbubble's shell to collapse back down to its equilibrium position was determined by applying the momentum balance law and the quasistatic phase's inflated radial deformation as an initial condition.It is important to recognize that the microbubble collapse does not end once it has reached its stress free state but continues to oscillate about its stress free state indefinitely due to the lack of viscous damping.Key material parameters such as the thickness of the shell, its Poisson ratio and the shell's shear modulus were varied to determine their influence on the collapse rate of the shell.A typical redimensionalized collapse time for a shelled microbubble of interior radius 1μm, a shell density of ρ = 1100kgm −3 , a redimensionalized shear modulus of μ = 20 MPa and a Poisson ratio of ν = 0.48, subjected to a redimensionalized stress load of p = 200 kPa, is of the order t * = 0.47 × (ρ o R 2 I /μ) ≈ 3.5ns where t * = 0.47 is the nondimensionalized collapse time.Shells with a larger shear modulus possessed faster collapse times.As the thickness of the shell increased the collapse time of the shell increased in an approximately linear manner.Shells with a larger Poisson ratio had smaller initial radial deformations and therefore exhibited faster collapse times.It would be extremely advantageous to have a numerical solution for the nonlinear model to not only validate the linearized model's solution but also to enable us to extend this study into the regimes of higher, non-perturbative stress loads.However, the development of a numerical solution for the nonlinear collapse model is a formidable challenge and is a possible area for future study.Current literature pertaining to shelled microbubble dynamics focusses on the forcing of gas loaded shells via an external pressure load such as an ultrasound signal.Studies such as Church (1995), Marmottant et al. (2005), Paul et al. (2010), Doinikov & Bouakaz (2011), Doinikov & Dayton (2007) and Doinikov et al. (2009) do not consider the effects of varying material parameters of the thin shell specifically its Poisson ratio and its shear modulus.Our model highlights the very small dependency between the collapse time and the thickness of the shell.Unfortunately, there is no relevant experimental data at present to support this model.Müller & Stannarius (2009) consider a collapsing and unfolding shelled smectic A liquid crystal millibubble which is collapsing as a consequence of the shell's surface tension whereas the mechanism responsible for the collapse of the shell in our model is the change in boundary conditions at the inner surface of the shell.Such a change in boundary conditions results in unloading at the inner shell's surface and changes the stress profiles within the shell.The main contributions of this study are (a) the nonlinear framework which can be a template for investigating other shell materials by choosing appropriate strain energy density functions (b) the nonlinear differential equation which models the quasistatic phase (no doubt this has a rich structure and so can be studied for its large deformation behaviour), and (c) the qualitative insight that it gives us into how the material parameters such as the shear modulus, the thickness of the shell and the Poisson ratio influence the collapse time of the shell.Hopefully such models will aid soft matter scientists' understanding of UCA localized drug delivery and gene therapy.To improve on this study, future work will consider the numerical solution for the nonlinear collapse model, the influence of viscoelasticity and an external stress load applied to the outer surface of the shell.

Fig. 4 .
Fig. 4. Graph of g(Y , t) at the nondimensionalized inner and outer radius of a collapsing shell versus the nondimensionalized time where ν = 0.48, β = 12, Y I = 1 and μ = 1 and the initial thickness is Y O − Y I = 0.02.This is calculated using equations (54) and (59).

Fig. 5 .
Fig. 5. Graph of the redimensionalized collapse time of the shell versus the redimensionalized shear modulus where β = 12, ν = 0.48, and the initial stress free thickness is R O − R I = 0.02nm.The density of the shell ρ o = 1100kgm −3(Falou et al., 2010) and its inner radius is R I = 1μm.Typical shear modulus values are within the range of μ = 20 → 100MPa(Falou et al., 2010;Tu et al., 2009).The collapse time is calculated from the zeros of the determinant of the matrix M(s) in equation (54) and is

Fig. 7 .
Fig. 7. Graph of the nondimensionalized collapse time of the shell versus the Poisson ratio of the shell's material for an initial nondimensionalized stress load of p = 0.01 where μ = 1 and the initial stress thickness is Y O − Y I = 0.02.This is calculated from the zeros of the determinant of the matrix M(s) in equation (54).