Likely cavitation and radial motion of stochastic elastic spheres

TL;DR

This paper investigates cavitation and radial motion in stochastic elastic spheres, revealing how probabilistic material properties influence bifurcation types and dynamic behaviors under various loading conditions.

Contribution

It introduces the analysis of cavitation in stochastic neo-Hookean spheres, including composites and inhomogeneous cases, highlighting the probabilistic effects on bifurcation and dynamic responses.

Findings

Center material determines critical cavitation load.

Supercritical bifurcation with stable cavitation occurs under dead-load in static spheres.

Subcritical bifurcation with snap cavitation is possible in composites and inhomogeneous spheres.

Abstract

The cavitation of solid elastic spheres is a classical problem of continuum mechanics. Here, we study this problem within the context of "stochastic elasticity" where the constitutive parameters are characterised by probability density functions. We consider homogeneous spheres of stochastic neo-Hookean material, composites with two concentric stochastic neo-Hookean phases, and inhomogeneous spheres of locally neo-Hookean material with a radially varying parameter. In all cases, we show that the material at the centre determines the critical load at which a spherical cavity forms there. However, while under dead-load traction, a supercritical bifurcation, with stable cavitation, is obtained in a static sphere of stochastic neo-Hookean material, for the composite and radially inhomogeneous spheres, a subcritical bifurcation, with snap cavitation, is also possible. For the dynamic spheres, oscillatory motions are produced under suitable dead-load traction, such that a cavity forms and expands to a maximum radius, then collapses again to zero periodically, but not under impulse traction. Under a surface impulse, a subcritical bifurcation is found in a static sphere of stochastic neo-Hookean material and also in an inhomogeneous sphere, whereas in composite spheres, supercritical bifurcations can occur as well. Given the non-deterministic material parameters, the results can be characterised in terms of probability distributions.

Equations

$$\left({\lambda }_i\frac{\partial \mathcal{W}}{\partial {\lambda }_i}-{\lambda }_j\frac{\partial \mathcal{W}}{\partial {\lambda }_j}\right)\left({\lambda }_i-{\lambda }_j\right)>0\text{mbox}if{\lambda }_i\ne {\lambda }_j,i,j=1,2,3.$$

$$\mathcal{W}({\lambda }_1,{\lambda }_2,{\lambda }_3)=\frac{\mu }{2}\left({\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2-3\right),$$

$$\displaystyle\begin{cases}E\left[\mu\right]=\underline{\mu}>0,&\\ E\left[\log\ \mu\right]=\nu,&\mbox{such that $|\nu|<+\infty$}.\end{cases}$$

$$\underline{\mu }={\rho }_1{\rho }_2,\text{Var}[\mu ]={\rho }_1{\rho }_2^2,$$

$$g(\mu ;{\rho }_1,{\rho }_2)=\frac{{\mu }^{{\rho }_1-1}{e}^{-\mu /{\rho }_2}}{{\rho }_2^{{\rho }_1}\Gamma ({\rho }_1)},\text{mbox}for\mu >0\text{mbox}and{\rho }_1,{\rho }_2>0,$$

$$\Gamma (z)={\int }_0^{+\infty }{t}^{z-1}{e}^{-t}\text{d}t.$$

$$\rho \ddot{\text{x}}=\mathrm{div}\text{T}+\rho \text{b},$$

$$\text{T}={\text{T}}^T,$$

$$\ddot{\text{x}}=-\mathrm{grad}\xi .$$

$$\mathrm{curl}\ddot{\text{x}}=\text{0}.$$

$$\text{T}=-\rho \xi \text{I}+{\text{T}}^{(0)},$$

$$r^3=R^3+c^3,\theta =\Theta ,\varphi =\Phi ,$$

$$\textbf{0}=\mathrm{curl}\ \ddot{\textbf{x}}=\left[\begin{array}[]{c}(\partial\ddot{\phi}/\partial\theta)/r-(\partial\ddot{\theta}/\partial\phi)/(r\sin\theta)\\ (\partial\ddot{r}/\partial\phi)/(r\sin\theta)-\partial\ddot{\phi}/\partial r\\ \partial\ddot{\theta}/\partial r-(\partial\ddot{r}/\partial\theta)/r\end{array}\right],$$

$$-\frac{\partial \xi }{\partial r}=\ddot{r}=\frac{2c{\dot{c}}^2+c^2\ddot{c}}{r^2}-\frac{2c^4{\dot{c}}^2}{r^5},$$

$$-\xi =-\frac{2c{\dot{c}}^2+a^2\ddot{c}}{r}+\frac{c^4{\dot{c}}^2}{2r^4}=-r\ddot{r}-\frac{3}{2}{\dot{r}}^2.$$

$$\text{F}=\mathrm{diag}\left(\frac{R^2}{r^2},\frac{r}{R},\frac{r}{R}\right),$$

$${\lambda }_1=\frac{R^2}{r^2},{\lambda }_2={\lambda }_3=\frac{r}{R}.$$

$$\text{B}={\text{F}}^2=\mathrm{diag}\left(\frac{R^4}{r^4},\frac{r^2}{R^2},\frac{r^2}{R^2}\right),$$

$$\begin{array}{rl}{I}_1=& \mathrm{tr}(\text{B})=\frac{R^4}{r^4}+2\frac{r^2}{R^2},\\ I_2=& \frac{1}{2}\left[{\left(\mathrm{tr}\text{B}\right)}^2-\mathrm{tr}\left({\text{B}}^2\right)\right]=\frac{r^4}{R^4}+2\frac{R^2}{r^2},\\ I_3=& \det \text{B}=1.\end{array}$$

$${\text{T}}^{(0)}=-p^{(0)}\text{I}+{\beta }_1\text{B}+{\beta }_{-1}{\text{B}}^{-1},$$

$${\beta }_1=2\frac{\partial W}{\partial I_1},{\beta }_{-1}=-2\frac{\partial W}{\partial I_2},$$

$$\begin{array}{rl}{T}_{rr}^{(0)}& =-p^{(0)}+{\beta }_1\frac{R^4}{r^4}+{\beta }_{-1}\frac{r^4}{R^4},\\ T_{\theta \theta }^{(0)}& =T_{rr}^{(0)}+\left({\beta }_1-{\beta }_{-1}\frac{r^2}{R^2}\right)\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right),\\ T_{\varphi \varphi }^{(0)}& =T_{\theta \theta }^{(0)}.\end{array}$$

$$\frac{\partial T_{rr}^{(0)}}{\partial r}=2\frac{T_{\theta \theta }^{(0)}-T_{rr}^{(0)}}{r}.$$

$$T_{rr}^{(0)}(r,t)=2\int \left({\beta }_1-{\beta }_{-1}\frac{r^2}{R^2}\right)\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{dr}{r}+\psi (t),$$

$$\begin{array}{rl}{T}_{rr}(r,t)& =-\rho \left(\frac{c^2\ddot{c}+2c{\dot{c}}^2}{r}-\frac{c^4{\dot{c}}^2}{2r^4}\right)+2\int \left({\beta }_1-{\beta }_{-1}\frac{r^2}{R^2}\right)\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{dr}{r}+\psi (t),\\ T_{\theta \theta }(r,t)& =T_{rr}(r,t)+\left({\beta }_1-{\beta }_{-1}\frac{r^2}{R^2}\right)\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right),\\ T_{\varphi \varphi }(r,t)& =T_{\theta \theta }(r,t).\end{array}$$

$$\begin{array}{rl}{T}_1(t)-T_2(t)& =\rho \left[\left(c^2\ddot{c}+2c{\dot{c}}^2\right)\left(\frac{1}{c}-\frac{1}{b}\right)-\frac{c^4{\dot{c}}^2}{2}\left(\frac{1}{c^4}-\frac{1}{b^4}\right)\right]+2{\int }_c^b\mu \left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{dr}{r}\\ & =\rho \left[\left(c\ddot{c}+2{\dot{c}}^2\right)\left(1-\frac{c}{b}\right)-\frac{{\dot{c}}^2}{2}\left(1-\frac{c^4}{b^4}\right)\right]+2{\int }_c^b\mu \left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{dr}{r}\\ & =\rho B^2\left[\left(\frac{c}{B}\frac{\ddot{c}}{B}+2\frac{{\dot{c}}^2}{B^2}\right)\left(1-\frac{c}{b}\right)-\frac{{\dot{c}}^2}{2B^2}\left(1-\frac{c^4}{b^4}\right)\right]+2{\int }_c^b\mu \left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{dr}{r}.\end{array}$$

$$u=\frac{r^{3}}{R^{3}}=\frac{r^{3}}{r^{3}-c^{3}},\qquad x=\frac{c}{B},\$$

$$\begin{array}{rl}\left(\frac{c}{B}\frac{\ddot{c}}{B}+2\frac{{\dot{c}}^2}{B^2}\right)\left(1-\frac{c}{b}\right)& -\frac{{\dot{c}}^2}{2B^2}\left(1-\frac{c^4}{b^4}\right)\\ & =\left(\ddot{x}x+2{\dot{x}}^2\right)\left[1-{\left(1+\frac{1}{x^3}\right)}^{-1/3}\right]-\frac{{\dot{x}}^2}{2}\left[1-{\left(1+\frac{1}{x^3}\right)}^{-4/3}\right]\\ & =\left(\ddot{x}x+\frac{3}{2}{\dot{x}}^2\right)\left[1-{\left(1+\frac{1}{x^3}\right)}^{-1/3}\right]-\frac{{\dot{x}}^2}{2}\frac{1}{x^3}{\left(1+\frac{1}{x^3}\right)}^{-4/3}\\ & =\frac{1}{2x^2}\frac{d}{dx}\left\{{\dot{x}}^2{x}^3\left[1-{\left(1+\frac{1}{x^3}\right)}^{-1/3}\right]\right\}\end{array}$$

$$\begin{array}{rl}{\int }_c^b\mu \left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{dr}{r}& ={\int }_c^b\mu \left[{\left(\frac{r^3}{r^3-c^3}\right)}^{2/3}-{\left(\frac{r^3-c^3}{r^3}\right)}^{4/3}\right]\frac{dr}{r}\\ & =\frac{1}{3}{\int }_{x^3+1}^{\infty }\mu \frac{1+u}{u^{7/3}}du.\end{array}$$

$$2x^2\frac{T_1(t)-T_2(t)}{\rho B^2}=\frac{d}{dx}\left\{{\dot{x}}^2{x}^3\left[1-{\left(1+\frac{1}{x^3}\right)}^{-1/3}\right]\right\}+\frac{4x^2}{3\rho B^2}{\int }_{x^3+1}^{\infty }\mu \frac{1+u}{u^{7/3}}du.$$

$$2x^2\frac{T_{rr}(b,t)}{\rho B^2}=\frac{d}{dx}\left\{{\dot{x}}^2{x}^3\left[1-{\left(1+\frac{1}{x^3}\right)}^{-1/3}\right]\right\}+\frac{4x^2}{3\rho B^2}{\int }_{x^3+1}^{\infty }\mu \frac{1+u}{u^{7/3}}du.$$

Full text

Key words: stochastic elasticity, finite strain analysis, cavitation, oscillations, quasi-equilibrated motion, probabilities.

“It is a problem of mechanics, highly complicated and irrelevant to probability theory except insofar as it forces us to think a little more carefully about how probability theory must be formulated if it is to be applicable to real situations.” - E.T. Jaynes (1996) [37]

1 Introduction

Cavitation in solids represents the void-formation within a solid under tensile loads, by analogy to the similar phenomenon observed in fluids. For rubber-like materials, this phenomenon was first reported by Gent and Lindley (1958) [24], whose experiments showed that rubber cylinders ruptured under relatively small tensile dead loads by opening an internal cavity. The nonlinear elastic analysis of Ball (1982) [7] provided the first theoretical explanation for the formation of a spherical cavity at the centre of a sphere of isotropic incompressible hyperelastic material in radially symmetric tension under prescribed surface displacements or dead loads. There, cavitation was treated as a bifurcation from the trivial state at a critical value of the surface traction or dis­placement, at which the trivial solution became unstable. This paved the way for numerous applied and theoretical studies devoted to this inherently nonlinear mechanical effect, which is not captured by the linear elasticity theory.

Many important results focusing on rubber-like materials are reviewed in [20, 23, 32]. Recent experimental studies regarding the onset, healing and growth of cavities in elastomers can be found in [34, 58, 59, 60]. Finite element simulations are presented in [40]. While most of this work was carried out in the setting of static deformations, in general, the elastodynamics of nonlinear hyperelastic bodies has been less investigated. The governing equations for dynamic large strain deformations leading to finite amplitude oscillations of spherical and cylindrical shells of homogeneous isotropic incompressible hyperelastic material, treated as particular cases of quasi-equilibrated motions defined by Truesdell (1962) [81], were reviewed in [82]. Under this type of dynamic deformation, the stress field is determined by the current configuration alone, and the body can be brought instantly to rest by applying a suitable pressure impulse on its boundary. Free and forced oscillations of spherical shells were first obtained in [29, 43, 84], using a similar approach as introduced in [41, 42] for the axially symmetric radial oscillations of infinitely long, isotropic incompressible cylindrical tubes. Oscillatory motion caused by the dynamic cavitation of a neo-Hookean sphere was considered in [13]. For a wide class of hyperelastic materials, the static and dynamic cavitation of homogeneous spheres were studied in [7]. For a hyperelastic sphere of Mooney-Rivlin material containing a cavity, the numerical solution to the nonlinear problem of large amplitude oscillations was presented in [6]. Theoretical and experimental studies of cylindrical and spherical shells of rubber-like material under external pressure were proposed in [85]. The finite amplitude radial oscillations of homogeneous isotropic incompressible hyperelastic spherical and cylindrical shells under a constant pressure difference between the inner and the outer surface were analysed theoretically in [11]. Radial oscillations of non-homogeneous thick-walled cylindrical and spherical shells of neo-Hookean material, with a material parameter varying continuously along the radial direction, were treated in [18]. In [83], the dynamic inflation of hyperelastic spherical membranes of Mooney-Rivlin material subjected to a uniform step pressure was considered, and the absence of damping was discussed. It was concluded that, as the amplitude and period of oscillations are strongly influenced by the rate of internal pressure, if the pressure was suddenly imposed and the inflation process was short, then sustained oscillations due to the dominant elastic effects could be observed. However, in many systems where the pressure is slowly increasing, strong damping would generally preclude oscillations [15]. The dynamic response of incompressible hyperelastic spherical shells subjected to periodic loading was examined in [62]. Radial oscillations of spherical shells of different hyperelastic materials were analysed in [10], where it was found that both the amplitude and period of oscillations decrease as the stiffness of the material increases. The influence of material constitutive law on the dynamic behaviour of spherical membranes was also considered in [65]. In [68], the nonlinear static and dynamic behaviour of a spherical membrane of neo-Hookean or Mooney-Rivlin material, subjected to a uniformly distributed radial pressure on the inner surface, was analysed, and the influence of the material constants was examined.

Despite these developments, there are important issues that remain unresolved. In particular, the quantification of uncertainties in material parameters and responses resulting from incomplete information remain largely unexplored [73]. However, for many materials, deterministic approaches, which are based on average data values, can greatly underestimate or overestimate their properties, and stochastic representations accounting also for data dispersion are needed to improve assessment and predictions [25, 33, 39, 55, 57, 80].

For nonlinear elastic materials, stochastic hyperelastic models have recently been described by strain-energy densities where the parameters are characterised by probability density functions (see [75, 76, 77, 78, 79] and also [50]). These are advanced phenomenological models, for which the parameters can be inferred from macroscopic load-deformation experiments, as shown, for example, in [19]. They rely on the finite elasticity theory as prior information, and on the notion of entropy (or uncertainty) [66, 74, 35, 36, 37] to enable the propagation of uncertainties from input data to output quantities of interest [72], and can also be incorporated into Bayesian approaches [8, 45] for model selection [19, 50, 55, 64].

To demonstrate the effect of probabilistic model parameters on predicted mechanical responses, for stochastic homogeneous incompressible hyperelastic bodies, we have analysed theoretically, so far, the static cavitation of a sphere under uniform tensile dead load [48], the inflation of pressurised spherical and cylindrical shells [46], the classic problem of the Rivlin cube [51], the radial oscillatory motion of cylindrical and spherical shells under a uniform pressure impulse [47], and the rotation and perversion of stochastic incompressible anisotropic hyperelastic cylindrical tubes [52]. In [48], in addition to the well-known case of stable cavitation post-bifurcation at the critical dead load, we show, for the first time, the existence of unstable (snap) cavitation for some (deterministic or stochastic) isotropic incompressible materials satisfying Baker-Ericksen (BE) inequalities [5]. In general, these problems, for which the elastic solutions are given explicitly, can offer significant insight into how the uncertainties in input parameters are propagated to output quantities. To investigate the effect of probabilistic parameters in the case of more complex geometries and loading conditions, numerical approaches have been proposed in [78, 79]. A similar stochastic methodology can be developed to study instabilities in other material systems (e.g., liquid crystal elastomers [16, 22, 86]), and may lead to more accurate assessment and prediction in many application areas [14, 30].

In this paper, we examine the static and dynamic cavitation and post-cavitation behaviour of an incompressible sphere of stochastic hyperelastic material, subject to either a uniform tensile surface dead load (which is constant in the reference configuration) or an impulse traction (which is kept constant in the current configuration), and each of them is applied uniformly in the radial direction. Here, cavitation represents the opening of a spherical cavity at the centre of the sphere as a bifurcated solution from the trivial solution, which becomes unstable at a critical value of the applied load. Moreover, the two different types of loads will lead to different static and dynamic solutions. First, in Section 2, we provide a summary of the stochastic elasticity prerequisites, where, in addition to the stochastic homogeneous material models, for which the elastic parameters are spatially-independent random variables at a macroscopic (non-molecular) level, we introduce a class of stochastic inhomogeneous models described by spatially-dependent strain-energy functions for which the parameters are non-Gaussian random fields. We then recall the notion of quasi-equilibrated motion, where the deformation field is circulation preserving and, at every time instant, the deformed state is possible also as a static equilibrium state under the same forces. Section 3 is devoted to the explicit analysis of cavitation and finite amplitude oscillations of a homogeneous sphere of stochastic neo-Hookean material. In Sections 4 and 5, we extend the analysis, respectively, to the cavitation and finite amplitude oscillations of a composite formed from two concentric homogeneous spheres of different stochastic neo-Hookean material, and to inhomogeneous spheres of locally neo-Hookean material with a radially varying parameter. In each case, the dynamic problem is considered first, then the static problem is treated by reducing the quasi-equilibrated motion to a static equilibrium state. In all cases, we show that the material at the centre determines the critical load at which a spherical cavity forms. However, there are important differences in the post-cavitation nonlinear elastic responses. While under dead-load traction, supercritical bifurcation, with stable cavitation, is obtained in a static sphere of stochastic neo-Hookean material, for the composite and radially inhomogeneous spheres, a subcritical bifurcation, with snap cavitation, is also possible. For the dynamic spheres, oscillatory motions are produced under suitable dead-load traction, such that a cavity forms and expands to a maximum radius, then collapses again to zero periodically, but not under impulse traction. Under a surface impulse, subcritical bifurcation is found in a static sphere of stochastic neo-Hookean material and also in an inhomogeneous sphere, whereas in composite spheres, supercritical bifurcation can occur as well. As the input material parameters are non-deterministic, the results can be characterised in terms of probability distributions. In Section 6, we draw concluding remarks and highlight potential extensions of the stochastic approach.

2 Prerequisites

In this section, we introduce a class of stochastic isotropic incompressible inhomogeneous hyperelastic models, where the constitutive parameter is spatially-dependent. This can be regarded as an extension of the stochastic homogeneous models in [19, 46, 47, 48, 50, 51], where both the mean value and variance of model parameters were constant in space, and also as a particular case of the inhomogeneous anisotropic ones discussed in [78, 79], which assumed heterogeneity in variance.

As both static and dynamic spheres of either homogeneous or inhomogeneous stochastic material will be considered in the next sections, we also recall the notion of quasi-equilibrated motion defined in [81, 82]. Quasi-equilibrated motions were first treated by us in [47], where the finite amplitude oscillations of cylindrical and spherical shells of stochastic isotropic incompressible homogeneous hyperelastic material under uniform pressure impulse were analysed. This type of motions can be reduced to static equilibrium states at any time instant under the given load.

2.1 Stochastic hyperelastic models

We are concerned here with the large strain analysis of spherical bodies of stochastic incompressible isotropic hyperelastic material, and rely on the following physically realistic assumptions [19, 46, 47, 48, 50, 51]:

• (A1)

Material objectivity, stating that constitutive equations must be invariant under changes of frame of reference. This requires that the scalar strain-energy function, $W=W(\textbf{F})$, depending only on the deformation gradient F, with respect to the reference configuration, is unaffected by a superimposed rigid-body transformation (which involves a change of position) after deformation, i.e., $W(\textbf{R}^{T}\textbf{F})=W(\textbf{F})$, where $\textbf{R}\in SO(3)$ is a proper orthogonal tensor (rotation). Material objectivity is guaranteed by defining strain-energy functions in terms of the scalar invariants.

• (A2)

Material isotropy, requiring that the strain-energy function is unaffected by a superimposed rigid-body transformation prior to deformation, i.e., $W(\textbf{F}\textbf{Q})=W(\textbf{F})$, where $\textbf{Q}\in SO(3)$. For isotropic materials, the strain-energy function is a symmetric function of the principal stretches $\{\lambda_{i}\}_{i=1,2,3}$ of F, i.e., $W(\textbf{F})=\mathcal{W}(\lambda_{1},\lambda_{2},\lambda_{3})$.

• (A3)

Baker-Ericksen (BE) inequalities, which state that the greater principal (Cauchy) stress occurs in the direction of the greater principal stretch [5, 44], and take the equivalent form

[TABLE]

When any two principal stretches are equal, the strict inequality “$>$” in (1) is replaced by “$\geq$”.

• (A4)

For any given finite deformation, at any point in the material, the shear modulus, $\mu$, and its inverse, $1/\mu$, are second-order random variables, i.e., they have finite mean value and finite variance [75, 76, 77, 78, 79].

Assumptions (A1)-(A3) are general requirements in isotropic finite elasticity [26, 49, 56, 82]. In particular, (A3) implies that the shear modulus is always positive, i.e., $\mu>0$ [49], while (A4) places random variables at the foundation on which to construct hyperelastic models [21, 37, 50, 54, 73].

To make our approach analytically tractable, we restrict our attention to a class of inhomogeneous models defined by the strain-energy density

[TABLE]

where the shear modulus, $\mu=\mu(R)>0$, is a random field variable depending on the radius $R$ of the undeformed sphere, and $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ are the principal stretch ratios. When $\mu$ is independent of $R$, the material model (2) reduces to the stochastic neo-Hookean model [75, 50].

For the shear modulus, $\mu=\mu(R)$, at any fixed $R$, we rely on the following available information

[TABLE]

Under the constraints (3), assumption (A4) is guaranteed [69, 70, 71, 73]. Then, by the maximum entropy principle, for any fixed $R$, $\mu$, follows a Gamma probability distribution with the shape and scale parameters, $\rho_{1}=\rho_{1}(R)>0$ and $\rho_{2}=\rho_{2}(R)>0$, respectively, satisfying

[TABLE]

where $\underline{\mu}$ is the mean value, $\text{Var}[\mu]=\|\mu\|^{2}$ is the variance, and $\|\mu\|$ is the standard deviation of $\mu$. The corresponding probability density function takes the form [1, 28, 38]

[TABLE]

where $\Gamma:\mathbb{R}^{*}_{+}\to\mathbb{R}$ is the complete Gamma function

[TABLE]

The word ‘hyperparameters’ is also used for the parameters $\rho_{1}$ and $\rho_{2}$ of the Gamma distribution, to distinguish them from the material parameter $\mu$ and other material constants [73, p. 8]. Examples where $\mu=\mu(R)$ takes specific forms are discussed in Section 5.

2.2 Quasi-equilibrated motion

For the large strain dynamic deformation of an elastic solid, Cauchy’s laws of motion (balance laws of linear and angular momentum) are governed by the following Eulerian field equations [82, p. 40],

[TABLE]

where $\rho$ is the material density, which is assumed constant, $\textbf{x}=\chi(\textbf{X},t)$ is the motion of the elastic solid, with velocity $\dot{\textbf{x}}=\partial\chi(\textbf{X},t)/\partial t$ and acceleration $\ddot{\textbf{x}}=\partial^{2}\chi(\textbf{X},t)/\partial t^{2}$, $\textbf{b}=\textbf{b}(\textbf{x},t)$ is the body force, $\textbf{T}=\textbf{T}(\textbf{x},t)$ is the Cauchy stress tensor, and the superscript $T$ defines the transpose.

To obtain possible dynamical solutions, one can solve Cauchy’s equation for particular motions, or generalise known static solutions to dynamical forms using the universal notion of quasi-equilibrated motion, which is defined as follows.

Definition 2.1

[82, p. 208]** A quasi-equilibrated motion, $\textbf{x}=\chi(\textbf{X},t)$, is the motion of an incompressible homogeneous elastic solid subject to a given body force, $\textbf{b}=\textbf{b}(\textbf{x},t)$, whereby, for each value of $t$, $\textbf{x}=\chi(\textbf{X},t)$ defines a static deformation that satisfies the equilibrium conditions under the body force $\textbf{b}=\textbf{b}(\textbf{x},t)$.

Theorem 2.2

[82, p. 208]** (see also the proof in [47]) A quasi-equilibrated motion, $\textbf{x}=\chi(\textbf{X},t)$, of an incompressible homogeneous elastic solid subject to a given body force, $\textbf{b}=\textbf{b}(\textbf{x},t)$, is dynamically possible, subject to the same body force, if and only if the motion is circulation preserving with a single-valued acceleration potential $\xi$, i.e.,

[TABLE]

For the condition (9) to be satisfied, it is necessary that

[TABLE]

Then, the Cauchy stress tensor takes the form

[TABLE]

where $\textbf{T}^{(0)}$ is the Cauchy stress for the equilibrium state at time $t$ and $\textbf{I}=\text{diag}(1,1,1)$ is the identity tensor. In this case, the stress field is determined by the present configuration alone. In particular, the shear stresses in the motion are the same as those of the equilibrium state at time $t$.

An immediate consequence of the above theorem is that a quasi-equilibrated motion is dynamically possible under a given body force in all elastic materials, if at every time instant the deformation is a possible equilibrium state under that body force in all those materials. Quasi-equilibrated motions of isotropic materials subject to surface tractions alone are obtained by taking the arbitrary constant in those deformations to be arbitrary functions of time. Under this type of motion, a body can be brought instantly to rest by applying a suitable pressure impulse on its boundary [82, p. 209]. Examples include the homogeneous motions that are possible in all homogeneous incompressible materials, and also those considered by us in the next sections.

2.3 Radial motion of stochastic hyperelastic spheres

Throughout our analysis, we assume that a sphere of stochastic hyperelastic material defined by (2) is subject to the following radially-symmetric dynamic deformation [47],

[TABLE]

where $(R,\Theta,\Phi)$ and $(r,\theta,\phi)$ are the spherical polar coordinates in the reference and current configuration, respectively, such that $0\leq R\leq B$, $B$ is the radius of the undeformed sphere, $c=c(t)\geq 0$ is the cavity radius to be calculated, and $b=b(t)=\sqrt[3]{B^{3}+c(t)^{3}}$ is the radius of the deformed sphere at time $t$ (see Figure 1).

By the governing equations (12),

[TABLE]

i.e., the condition (10) is valid for $\textbf{x}=(r,\theta,\phi)^{T}$. Therefore, (12) describes a quasi-equilibrated motion, such that

[TABLE]

where $\xi$ is the acceleration potential, satisfying (9). Integrating (14) gives

[TABLE]

For the deformation (12), the gradient tensor with respect to the polar coordinates $(R,\Theta,\Phi)$ is equal to

[TABLE]

and the principal stretch ratios are

[TABLE]

The corresponding Cauchy-Green tensor is

[TABLE]

with the principal invariants

[TABLE]

The associated Cauchy stress tensor then takes the form [27, pp. 87-91]

[TABLE]

where the scalar $p^{(0)}$, which is commonly referred to as the arbitrary hydrostatic pressure [26, pp. 286-287], [56, pp. 198-201], is the Lagrange multiplier for the internal constraint $I_{3}=1$ of incompressibility (i.e., all deformations are isochoric for incompressible materials) [82, pp. 71-72], and the coefficients

[TABLE]

are nonlinear material parameters, with $W=W(I_{1},I_{2},I_{3})=\mathcal{W}(\lambda_{1},\lambda_{2},\lambda_{2})$, and $I_{1}$, $I_{2}$, $I_{3}$ defined by (19). For this stress tensor, the principal components at time $t$ are as follows,

[TABLE]

As the stress components depend only on the radius $r$, the system of equilibrium equations reduces to

[TABLE]

Hence, by (22) and (23), the radial Cauchy stress for the equilibrium state at $t$ is equal to

[TABLE]

where $\psi=\psi(t)$ is an arbitrary function of time. Substitution of (15) and (24) into (11) then gives the following principal Cauchy stresses at time $t$,

[TABLE]

When $R^{2}/r^{2}\to 1$, all the stress components defined in (25) are equal.

3 Cavitation and radial motion of homogeneous spheres

We first analyse the static and dynamic radially-symmetric deformations of a stochastic homogeneous hyperelastic sphere, where the shear modulus $\mu$ of the stochastic model given by (2) is a space-invariant random variable. Hence, the model is of stochastic neo-Hookean type [75], with the shear modulus $\mu$ characterised by the Gamma distribution defined by (5). The radial motion of spheres of neo-Hookean material were treated deterministically in [13], while inhomogeneous cylindrical and spherical shells of neo-Hookean-like material with a radially varying material parameter were analysed in [18]. For homogeneous spheres of stochastic isotropic incompressible hyperelastic material, the static cavitation under uniform tensile dead load was investigated in [48], while radial oscillatory motions of homogeneous cylindrical and spherical shells of stochastic Mooney-Rivlin and neo-Hookean material, respectively, were treated in [47].

3.1 Oscillatory motion of a stochastic neo-Hookean sphere under dead-load traction

For a sphere of stochastic neo-Hookean material subject to the quasi-equilibrated motion (12), we denote the inner and outer radial pressures acting on the curvilinear surfaces, $r=c(t)$ and $r=b(t)$ at time $t$, as $T_{1}(t)$ and $T_{2}(t)$, respectively [82, pp. 217-219].

Evaluating $T_{1}(t)=-T_{rr}(c,t)$ and $T_{2}(t)=-T_{rr}(b,t)$, using (25), with $r=c$ and $r=b$, respectively, and subtracting the results, gives

[TABLE]

By setting the notation

[TABLE]

we can then rewrite

[TABLE]

and

[TABLE]

Thus (26) is equivalent to

[TABLE]

Next, assuming that the cavity surface is traction-free, $T_{1}(t)=0$ and (28) is equivalent to

[TABLE]

We now denote

[TABLE]

and set the uniform dead-load traction (see also eq. (2.7) of [13])

[TABLE]

where $p_{0}$ is constant, and $x=x(t)$ is the dimensionless cavity radius, as denoted in (27).

In this case, integrating (29) once gives

[TABLE]

with $H(x)$ defined by (30), and

[TABLE]

Then, after setting the initial conditions $x_{0}=x(0)=0$ and $\dot{x}_{0}=\dot{x}(0)=0$, equation (32) takes the form

[TABLE]

From (34), we obtain the velocity

[TABLE]

It is useful to note that this nonlinear elastic system is analogous to the motion of a point mass with energy

[TABLE]

where the energy is $E=C$, the potential is given by $V(x)=H(x)-\frac{2p_{0}}{\rho B^{2}}\left(x^{3}+1\right)^{1/3}$ and the position-dependent mass is $m(x)=x^{3}\left[1-\left(1+\frac{1}{x^{3}}\right)^{-1/3}\right]$. Due to the constraints on the function $H$, the system has simple dynamics, and the only solutions of interest are either static or periodic solutions.

For the given sphere to undergo an oscillatory motion, the following equation,

[TABLE]

must have exactly two finite distinct positive roots $x_{1},x_{2}$ such that $0\leq x_{1}<x_{2}<\infty$, as shown in Figure 2. Then, the minimum and maximum radii of the cavity in the oscillation are $x_{1}B$ and $x_{2}B$, respectively, and the period of oscillation is equal to

[TABLE]

For the stochastic sphere, the amplitude and the period of the oscillation are random variables characterised by probability distributions.

Evaluating the integral in (30) gives

[TABLE]

and assuming that the shear modulus, $\mu$, which is a random variable, is bounded from below, i.e.,

[TABLE]

for some constant $\eta>0$, it follows that $H(0)=0$ and $\lim_{x\to\infty}H(x)=\infty$.

When $p_{0}\neq 0$, substitution of (39) in (37) gives

[TABLE]

Equation (41) has one solution at $x_{1}=0$, while the second solution, $x_{2}$, is a root of

[TABLE]

As the right-hand side of equation (42) is an increasing function of $x$, this equation has a solution, $x_{2}>0$, if and only if (see also eq. (2.7) of [13])

[TABLE]

Then,

[TABLE]

In this case, at the centre of the given sphere, a spherical cavity forms and expands until its radius reaches the value $c=x_{2}B$, where $x_{2}$ is the root of (42), then contracts again to zero radius and repeats the cycle.

The critical dead load for the onset of cavitation is

[TABLE]

Regarding the random shear modulus, by (40) and (43), for the motion to be oscillatory, the following condition must hold,

[TABLE]

Then, as $\mu$ follows a Gamma distribution, defined by (5), the probability distribution of oscillatory motions occurring is

[TABLE]

and that of non-oscillatory motions is

[TABLE]

For example, when $\mu$ satisfies the Gamma distribution with $\rho_{1}=405$ and $\rho_{2}=0.01$, shown in Figure 3, the probability distributions given by (47)-(48) are shown in Figure 4 (with blue lines for $P_{1}$ and red lines for $P_{2}$). Specifically, the interval $(0,5\underline{\mu})$, where $\underline{\mu}=\rho_{1}\rho_{2}=4.05$ is the mean value of $\mu$, was divided into $100$ steps, then for each value of $p_{0}$, $100$ random values of $\mu$ were numerically generated from the specified Gamma distribution and compared with the inequalities defining the two intervals for values of $p_{0}$. For the deterministic elastic sphere, the critical value $p_{0}=5\underline{\mu}/2=10.125$ strictly divides the cases of oscillations occurring or not. For the stochastic problem, for the same critical value, there is, by definition, exactly 50% chance that the motion is oscillatory, and 50% chance that is not. To increase the probability of oscillatory motion ($P_{1}\approx 1$), one must apply a sufficiently small traction, $p_{0}$, below the expected critical point, whereas a non-oscillatory motion is certain to occur ($P_{2}\approx 1$) if $p_{0}$ is sufficiently large. However, the inherent variability in the probabilistic system means that there will also exist events where there is competition between the two cases.

We note that the Gamma distribution represented in Figure 3 is approximately a normal distribution. This is because, when $\rho_{1}$ much larger than $\rho_{2}$, the Gamma probability distribution is approximated by a normal distribution (see, e.g., [19, 46]). However, elastic moduli cannot be characterised by the normal distribution, since this distribution is defined on the entire real line whereas elastic moduli are typically positive. In practice, these moduli can meaningfully take on different values, corresponding to possible outcomes of experimental tests. Then, the principle of maximum entropy enables the construction of their probability distributions, given the available information. Approaches for the explicit derivation of probability distributions for the elastic parameters of stochastic homogeneous isotropic hyperelastic models calibrated to experimental data were developed in [50, 77].

In Figure 5, we illustrate the stochastic function $H(x)$, defined by (39), intersecting the curve $\frac{2p_{0}}{\rho B^{2}}\left[\left(x^{3}+1\right)^{1/3}-1\right]$, with $p_{0}=15$, to find the two distinct solutions of equation (37), and the associated velocity, given by (35), assuming that $\rho=1$, $B=1$, and $\mu$ is drawn from the Gamma distribution with $\rho_{1}=405$ and $\rho_{2}=0.01$ (see Figure 3). Each figure displays a probability histogram at each value of $x$. The histogram comprises of 1000 stochastic simulations and the colour bar defines the probability of finding a given value of $H(x)$, or of the associated velocity, respectively, at a given value of $x$. The dashed black line corresponds to the expected values based only on the mean value $\underline{\mu}=\rho_{1}\rho_{2}=4.05$, of $\mu$.

Throughout this paper, different simulations were produced by fixing the parameters given in each figure caption, and repeatedly drawing random samples from the underlying distribution. Our computer simulations were run in Matlab 2018a where we made specific use of inbuilt functions for random number generation.

3.2 Static deformation of a stochastic neo-Hookean sphere under dead-load traction

In view of our subsequent analysis, we review here the cavitation of a static sphere of stochastic neo-Hookean material, with the shear modulus $\mu$ following the Gamma distribution defined by (5). Incompressible spheres of different stochastic homogeneous hyperelastic material were treated in detail in [48] (see also [12] for the deterministic spheres). In this case, if the surface of the cavity is traction-free, then $T_{1}=0$ and (29) reduces to

[TABLE]

After evaluating the integral in (49), the required uniform dead-load traction at the outer surface, $R=B$, in the reference configuration, takes the form

[TABLE]

and increases as $x$ increases. The critical dead load for the onset of cavitation is then

[TABLE]

and is equal to that given by (45) for the dynamic sphere [7, 13].

To analyse the stability of this cavitation, we study the behaviour of the cavity opening, with radius $c$ as a function of $P$, in a neighbourhood of $P_{0}$. After differentiating the function given by (50), with respect to the dimensionless cavity radius $x=c/B$, we have

[TABLE]

i.e., the cavitation is stable (with 100% certainty), regardless of the material parameter $\mu$ [48]. For example, the post-cavitation stochastic behaviour of the static unit sphere when the shear modulus, $\mu$, follows a Gamma distribution with $\rho_{1}=405$ and $\rho_{2}=4.05/\rho_{1}=0.01$ is shown in Figure 6.

3.3 Non-oscillatory motion of a stochastic neo-Hookean sphere under impulse traction

For a sphere subject to the radially-symmetric dynamic deformation (12), an impulse (suddenly applied) traction, expressed in terms of the Cauchy stresses, is prescribed as follows [47],

[TABLE]

where $p_{0}$ is constant in time. Introducing the dimensionless cavity radius $x(t)=c(t)/B$ and setting the initial conditions $x_{0}=x(0)=0$ and $\dot{x}_{0}=\dot{x}(0)=0$, we obtain the following differential equation

[TABLE]

where $H(x)$ is given by (39). From (54), we obtain the velocity

[TABLE]

assuming that $\frac{p_{0}}{3}x^{3}-H(x)\geq 0$. For an oscillatory motion to occur, the equation

[TABLE]

must have exactly two distinct solutions, $x=x_{1}$ and $x=x_{2}$, such that $0\leq x_{1}<x_{2}<\infty$.

When $p_{0}\neq 0$, substitution of (39) in (56) gives

[TABLE]

which has a solution at $x_{1}=0$. The second solution, $x_{2}>0$, is a root of

[TABLE]

As the right-hand side of equation (58) is a decreasing function of $x$, this equation has a solution, $x_{2}>0$, if and only if

[TABLE]

and

[TABLE]

However, since this implies

[TABLE]

the sphere cannot oscillate when $p_{0}$ is constant in time, as assumed in (53) (see Figure 7).

By (53) and (58), the tensile traction, $T=T_{rr}(b,t)$, takes the form

[TABLE]

and decreases as $x$ increases. Thus, the critical tension for the onset of cavitation is

[TABLE]

as found also under dead loading.

3.4 Static deformation of a stochastic neo-Hookean sphere under impulse traction

For the static sphere subject to a uniform constant surface load in the current configuration, given in terms of the Cauchy stresses, we have

[TABLE]

which decreases as $x$ increases. Then, the critical tension for cavitation initiation is also given by (63).

The post-cavitation behaviour can be inferred from the sign of the derivative of $T$ in a neighbourhood of $T_{0}$. By differentiating (64) with respect to $x$, we obtain

[TABLE]

i.e., the cavitation is unstable, regardless of the material parameter $\mu$ [7]. The post-cavitation stochastic behaviour of the static homogeneous sphere is shown in Figure 8, for the unit sphere (with $B=1$) of stochastic neo-Hookean material, where the shear modulus, $\mu$, is drawn from a Gamma distribution with shape and scale parameters $\rho_{1}=405$ and $\rho_{2}=4.05/\rho_{1}=0.01$, respectively (see also Figure 2 of [7]).

4 Cavitation and radial motion of concentric homogeneous spheres

Next, we extend the approach developed in the previous section to investigate the behaviour under the quasi-equilibrated radial motion (12) of a composite formed from two concentric homogeneous spheres (see Figure 9). We restrict our attention to composite spheres with two stochastic neo-Hookean phases, similar to those containing two concentric spheres of different neo-Hookean material treated deterministically in [31] and [67]. In this case, we define the following strain-energy function,

[TABLE]

where $0<R<A$ and $A<R<B$ denote the radii of the inner and outer sphere in the reference configuration, and the corresponding shear moduli $\mu^{(1)}$ and $\mu^{(2)}$ are spatially-independent random variables characterised by the Gamma distributions $g(u;\rho^{(1)}_{1},\rho^{(1)}_{2})$ and $g(u;\rho^{(2)}_{1},\rho^{(2)}_{2})$, defined by (5), respectively.

4.1 Oscillatory motion of a sphere with two stochastic neo-Hookean phases under dead-load traction

For the composite sphere, we denote the radial pressures acting on the curvilinear surfaces $r=c(t)$ and $r=b(t)$ at time $t$ as $T_{1}(t)$ and $T_{2}(t)$, respectively, and impose the continuity condition for the stress components across their interface, $r=a(t)$. By analogy to (26), we obtain

[TABLE]

Using the notation (27), we can write (67) equivalently as follows,

[TABLE]

In this case, we denote

[TABLE]

and, after evaluating the integrals, we obtain the equivalent form

[TABLE]

Then, substitution of (70) in (37) gives

[TABLE]

This equation has one solution at $x_{1}=0$, while the second solution, $x_{2}>0$, is a root of

[TABLE]

Next, we expand $p_{0}$, given by (72), to the first order in $x^{3}$ to obtain

[TABLE]

The critical tension for cavity initiation is then

[TABLE]

and a comparison with (45) shows that it is the same as for the homogeneous sphere made of the same material as the inner sphere.

Considering the random shear parameters of the two concentric spheres, we distinguish the following two cases:

• (i)

If

[TABLE]

then the right-hand side of (73) is an increasing function of $x$, and

[TABLE]

By (76), for the motion of the composite sphere to be oscillatory, the shear modulus of the inner sphere must satisfy

[TABLE]

where the upper bound follows from (76), while the lower bound is assumed to be strictly positive, i.e., $\eta>0$. Thus, the probability distribution of oscillatory motions occurring is given by (47) and that of non-oscillatory motions by (48), with $g(u;\rho^{(1)}_{1},\rho^{(1)}_{2})$ instead of $g(u;\rho_{1},\rho_{2})$ (see Figure 4).

As the material parameters, $\mu^{(1)}$ and $\mu^{(2)}$, are positive, (75) is possible for any values $0<A<B$. In particular, if

[TABLE]

then (75) always holds, regardless of the material constants.

• (ii)

If

[TABLE]

then the right-hand side in (73) decreases as $x$ increases, hence

[TABLE]

As the material constants, $\mu^{(1)}$ and $\mu^{(2)}$, are positive, (79) is possible if and only if

[TABLE]

However, in this case,

[TABLE]

and the sphere cannot oscillate.

In Figure 11, we represent the stochastic function $H(x)$, defined by (70), intersecting the curve $\frac{2p_{0}}{\rho B^{2}}\left[\left(x^{3}+1\right)^{1/3}-1\right]$ , with $p_{0}=15$, to find the two distinct solutions to (37), and the associated velocity, given by (35), for a composite sphere with two concentric stochastic neo-Hookean phases, with radii $A=1/2$ and $B=1$, respectively, when the shear modulus of the inner phase, $\mu^{(1)}$, follows a Gamma distribution with $\rho^{(1)}_{1}=405$, $\rho^{(1)}_{2}=4.05/\rho^{(1)}_{1}=0.01$, while the shear modulus of the outer phase, $\mu^{(2)}$, is drawn from a Gamma distribution with $\rho^{(2)}_{1}=405$, $\rho^{(2)}_{2}=4.2/\rho^{(2)}_{1}$ (see Figure 10).

4.2 Static deformation of a sphere with two stochastic neo-Hookean phases under dead-load traction

For the static composite sphere, if the surface of the cavity is traction-free, then $T_{1}=0$ and (68) reduces to

[TABLE]

After evaluating the integrals in (83), the required dead-load traction at the outer surface, $R=B$, in the reference configuration, is equal to

[TABLE]

The critical dead load for the cavity formation is then

[TABLE]

and is the same as for the homogeneous sphere made entirely of the material of the inner sphere [31], as found also in the dynamic case.

To study the stability of this cavitation, we examine the behaviour of the cavity opening, with radius $c$ as a function of $P$, in a neighbourhood of $P_{0}$. Expanding the expression of $P$, given by (84), to the first order in $x^{3}$, we obtain [31]

[TABLE]

After differentiating the above function with respect to the dimensionless cavity radius $x=c/B$, we have

[TABLE]

Then, $\text{d}P/\text{d}x\to 0$ as $x\to 0_{+}$, and the bifurcation at the critical dead load, $P_{0}$, is supercritical (respectively, subcritical) if $\text{d}P/\text{d}x>0$ (respectively, $\text{d}P/\text{d}x<0$) for arbitrarily small $x$. From (87), the following two cases are distinguished (see also the discussion in [31]):

• (i)

If (75) is valid, then the solution presents itself as a supercritical bifurcation, and the cavitation is stable, in the sense that a new bifurcated solution exists locally for values of $P>P_{0}$, and the cavity radius monotonically increases with the applied load post-bifurcation. In particular, when $\mu^{(1)}=\mu^{(2)}$, the problem reduces to the case of a homogeneous sphere made of neo-Hookean material, for which stable cavitation is known to occur [48].

• (ii)

If (79) holds, then the bifurcation is subcritical, with the cavitation being unstable, in the sense that the required dead load starts to decrease post-bifurcation, causing a snap cavitation, i.e., a sudden jump in the cavity opening.

Thus, on the one hand, when (78) is satisfied, (75) is valid regardless of the values of $\mu^{(1)}$ and $\mu^{(2)}$, and the cavitation is guaranteed to be stable (with 100% certainty). On the other hand, if (81) holds, then, as either of the two inequalities (75) and (79) is possible, the probability distribution of stable cavitation in the composite sphere is equal to

[TABLE]

while that of unstable (snap) cavitation is

[TABLE]

For example, setting $A=1/2$ and $B=1$, if $\rho^{(1)}_{1}=405$, $\rho^{(1)}_{2}=4.05/\rho^{(1)}_{1}=0.01$ for the inner sphere, and $\rho^{(2)}_{1}=405$, $\rho^{(2)}_{2}=4.2/\rho^{(2)}_{1}$ for the outer sphere (see Figure 10), then the probability distributions given by equations (88)-(89) are illustrated numerically in Figure 12 (with blue lines for $P_{1}$ and red lines for $P_{2}$). For the numerical realisations in these plots, the interval $(0,8.1)=\left(0,2\underline{\mu}^{(1)}\right)$ was discretised into $100$ representative points, then for each value of $\mu^{(2)}$, $100$ random values of $\mu^{(1)}$ were numerically generated from the specified Gamma distribution and compared with the inequalities defining the two intervals for values of $\mu^{(2)}$. For the deterministic elastic case, which is based on the mean value of the shear modulus, $\underline{\mu}^{(1)}=\rho^{(1)}_{1}\rho^{(1)}_{2}=4.05$, the critical value of $27\underline{\mu}^{(1)}/28=3.9054$ strictly separates the cases where cavitation instability can, and cannot, occur. However, in the stochastic case, the two states compete. For example, at the same critical value, there is, by definition, exactly 50% chance that the cavitation is stable, and 50% chance that is not.

The post-cavitation stochastic behaviours shown in Figure 13 correspond to two different static composite sphere where the shear modulus of the inner phase, $\mu^{(1)}$, follows a Gamma distribution with $\rho^{(1)}_{1}=405$, $\rho^{(1)}_{2}=4.05/\rho^{(1)}_{1}=0.01$, while the shear modulus of the outer phase, $\mu^{(2)}$, is drawn from a Gamma distribution with $\rho^{(2)}_{1}=405$, $\rho^{(2)}_{2}=3.6/\rho^{(2)}_{1}$ (in Figure 13-left) or $\rho^{(2)}_{1}=405$, $\rho^{(2)}_{2}=4.2/\rho^{(2)}_{1}$ (in Figure 13-right). In each case, unstable or stable cavitation is expected, respectively, but there is also about 5% chance that the opposite behaviour is presented.

4.3 Non-oscillatory motion of a sphere with two stochastic neo-Hookean phases under impulse traction

In this case, setting the initial conditions $x_{0}=x(0)=0$ and $\dot{x}_{0}=\dot{x}(0)=0$, we obtain equation (54), where $H(x)$ is given by (70). Then, substitution of (70) in (56) gives

[TABLE]

Equivalently, by (53), in terms of the tensile traction $T=T_{rr}(b,t)$,

[TABLE]

This equation has one solution at $x_{1}=0$, while the second solution, $x_{2}>0$, is a root of

[TABLE]

After expanding $T$, given by (92), to the first order in $x^{3}$, we have

[TABLE]

The associated critical tension for cavity initiation is

[TABLE]

and a comparison with that given by (63) shows that it is the same as for the homogeneous sphere made of the same material as the inner sphere.

We now distinguish the following two cases:

• (i)

If

[TABLE]

then the right-hand side of (93) is an increasing function of $x$, hence

[TABLE]

By (96), for the motion of the composite sphere to be oscillatory, the shear modulus of the inner sphere must satisfy

[TABLE]

where the lower bound is assumed to be strictly positive, i.e., $\eta>0$. Equivalently, by (53), in terms of the constant $p_{0}$,

[TABLE]

However, in this case,

[TABLE]

for all $x>0$, and the sphere does not oscillate (see Figure 14-left).

• (ii)

If

[TABLE]

then the right-hand side of (93) decreases as $x$ increases, hence

[TABLE]

or equivalently, by (53),

[TABLE]

However, in this case,

[TABLE]

and the sphere cannot oscillate (see Figure 14-right).

Note that, since the material parameters, $\mu^{(1)}$ and $\mu^{(2)}$, are positive, both (95) and (100) are possible for any values $0<A<B$.

4.4 Static deformation of a sphere with two stochastic neo-Hookean phases under impulse traction

For the static composite sphere subject to a uniform constant surface load in the current configuration, given in terms of the Cauchy stresses, the tensile traction takes the form

[TABLE]

The critical traction for the onset of cavitation is thus also given by (94).

To investigate the stability of this cavitation, we examine the behaviour of the cavity opening, with scaled radius $x$ as a function of $T$, in a neighbourhood of $T_{0}$. Expanding the expression of $T$, given by (104), to the first order in $x^{3}$, we have

[TABLE]

and

[TABLE]

Then, $\text{d}T/\text{d}x\to 0$ as $x\to 0_{+}$, and the bifurcation at the critical impulse load, $T_{0}$, is supercritical (respectively, subcritical) if $\text{d}T/\text{d}x>0$ (respectively, $\text{d}T/\text{d}x<0$) for arbitrarily small $x$. From (106), the following two cases arise:

• (i)

If (95) holds, then the solution after bifurcation is supercritical, and the cavitation is stable, in the sense that a new bifurcated solution exists locally for values of $T>T_{0}$, and the cavity radius monotonically increases with the applied load post-bifurcation.

• (ii)

If (100) holds, then the bifurcated solution is subcritical, with the cavitation being unstable, in the sense that the required dead load starts to decrease post-bifurcation, causing a snap cavitation, i.e., a sudden jump in the cavity opening. In particular, when $\mu^{(1)}=\mu^{(2)}$, the problem reduces to the case of a homogeneous sphere of neo-Hookean material, for which unstable cavitation always occurs (see Section 3).

As either of the two inequalities (95) and (100) is possible, the probability distribution of stable cavitation in the composite sphere is equal to (see Figure 15)

[TABLE]

while that of unstable (snap) cavitation is

[TABLE]

The post-cavitation stochastic behaviours shown in Figure 16 correspond to two different static composite sphere where the shear modulus of the inner phase, $\mu^{(1)}$, follows a Gamma distribution with $\rho^{(1)}_{1}=405$, $\rho^{(1)}_{2}=4.05/\rho^{(1)}_{1}=0.01$, while the shear modulus of the outer phase, $\mu^{(2)}$, is drawn from a Gamma distribution with $\rho^{(2)}_{1}=405$, $\rho^{(2)}_{2}=0.005$ (in Figure 16-left) or $\rho^{(2)}_{1}=405$, $\rho^{(2)}_{2}=0.02$ (in Figure 16-right). Note that, in each case, unstable or stable cavitation is expected, respectively, but there is other values of the parameters lead to the opposite behaviour.

5 Cavitation and radial motion of inhomogeneous spheres

We further examine the cavitation of radially inhomogeneous, incompressible spheres of stochastic hyperelastic material. Cavitation of radially inhomogeneous hyperelastic spheres of compressible or incompressible material in static equilibrium was investigated in [67]. Here, we adopt a neo-Hookean-like model, where the constitutive parameter varies continuously along the radial direction. Our inhomogeneous model is similar to those proposed in [18], where the dynamic inflation of spherical shells was treated explicitly.

In particular, we define the class of stochastic inhomogeneous hyperelastic models (2), with the shear modulus taking the form (see also eq. (29) of [18])

[TABLE]

where $\mu(R)>0$, for all $0\leq R\leq B$. For simplicity, we further assume that, for any fixed $R$, the mean value $\underline{\mu}$ of $\mu=\mu(R)$ is independent of $R$, (namely the expected value of $C_{2}$ is 0), whereas its variance, $\text{Var}[\mu]$ changes with $R$. Then, $C_{1}=\mu(0)>0$ is a single-valued (deterministic) constant and $C_{2}$ is a random value, defined by a given probability distribution.

Alternative modelling formulations where both the mean value and variance of the shear modulus may vary with $R$ are discussed at the end of this section.

When the mean value of the shear modulus $\mu(R)$, described by (109), does not depend on $R$, it follows that, for any fixed $R$,

[TABLE]

where $\text{Var}[C_{2}]$ is the variance of $C_{2}$. Note that, at $R=0$, the value $\mu(0)=C_{1}$ represents the material property at the centre of the sphere, which is a unique point, whereas at any fixed $R>0$, the standard deviation $\|\mu\|=\sqrt{\text{Var}[\mu]}$ is proportional to the volume $4\pi R^{3}/3$ of the sphere with radius $R$.

By (4) and (110), the hyperparameters of the corresponding Gamma distribution, defined by (5), take the form

[TABLE]

For example, one can choose two constant values, $C_{0}>0$ and $C_{1}>0$, and set the hyperparameters for the Gamma distribution at any given $R$ as follows,

[TABLE]

By (109), $C_{2}=\left(\mu(R)-C_{1}\right)B^{3}/R^{3}$ is the shifted Gamma-distributed random variable with mean value $\underline{C}_{2}=0$ and variance $\text{Var}[C_{2}]=\rho_{1}\rho_{2}^{2}B^{6}/R^{6}=C_{0}C_{1}$. When $R$ decreases towards zero, by (112), $\rho_{1}$ increases, while $\rho_{2}$ decreases, hence, the Gamma distribution converges to a normal distribution [19, 46], which degenerates at the centre of the sphere, where $R=0$, as the shear modulus there takes the deterministic value $C_{1}$. In Figure 17, we show these distributions when $\rho_{1}=405\cdot B^{6}/R^{6}$ and $\rho_{2}=0.01\cdot R^{6}/B^{6}$. For this case, by (109), $C_{1}=\underline{\mu}=\rho_{1}\rho_{2}=4.05$ and $C_{2}=\mu(B)-C_{1}$.

5.1 Oscillatory motion of a stochastic radially inhomogeneous sphere under dead-load traction

The shear modulus of the form given by (109) can be expressed equivalently as follows,

[TABLE]

where $u=r^{3}/R^{3}$ and $x=c/B$, as denoted in (27).

Also, writing the invariants given by (19) in the equivalent forms

[TABLE]

and substituting these in (21) gives

[TABLE]

Hence,

[TABLE]

and

[TABLE]

After calculations analogous to those in Section 3, but with $\beta_{1}-\beta_{-1}\frac{r^{2}}{R^{2}}$ instead of $\mu$ in (26), (28), and (29), we denote

[TABLE]

Equivalently, by (117),

[TABLE]

Then, by setting the dead load (31) and substituting (119) in (37), we obtain

[TABLE]

Equation (120) has a solution at $x_{1}=0$, while the second solution, $x_{2}>0$, is a root of

[TABLE]

To obtain the right-hand side of the above equation, we have assumed that, in (120), $x$ is sufficiently small, such that, after expanding the respective functions to the second order in $(x^{3}+1)^{1/3}$,

[TABLE]

and

[TABLE]

Next, expanding the right-hand side of the equation (121) to the first order in $x^{3}$ gives

[TABLE]

The critical dead load for the onset of cavitation is then

[TABLE]

and is the same as for the homogeneous sphere made entirely from the material found at the centre of the inhomogeneous sphere.

Considering the parameters $C_{1}$ and $C_{2}$, we have the following two cases:

• (i)

When $C_{1}>2C_{2}$, the right-hand side in (124) is an increasing function of $x$, and therefore,

[TABLE]

• (ii)

When $C_{1}<2C_{2}$, the right-hand side in (124) decreases as $x$ increases, hence,

[TABLE]

However, since

[TABLE]

the sphere cannot oscillate.

In Figure 18, we illustrate the stochastic function $H(x)$, defined by (119), intersecting the (dashed red) curve $\frac{2p_{0}}{\rho B^{2}}\left[\left(x^{3}+1\right)^{1/3}-1\right]$, with $p_{0}=15$, to find the two distinct solutions to (37), and the associated velocity, given by (35), for a unit sphere (with $B=1$) of stochastic radially inhomogeneous neo-Hookean-like material when the random field shear modulus $\mu(R)$ is given by (109), with $\rho_{1}=405/R^{6}$ and $\rho_{2}=0.01\cdot R^{6}$ (see Figure 17). In particular, $\rho_{1}=405$ and $\rho_{2}=0.01$ at $R=B$, hence, by (109), $C_{1}=\underline{\mu}=\rho_{1}\rho_{2}=4.05$ and $C_{2}=\mu(B)-C_{1}$. For these functions, although the mean values, represented by black dashed lines, are the same as for those depicted in Figure 5, their stochastic behaviours are quite different compared to the case of a stochastic neo-Hookean sphere. Specifically, as we observe from

5.2 Static deformation of a stochastic radially inhomogeneous sphere under dead-load traction

For the static inhomogeneous sphere, when the cavity surface is traction-free, at the outer surface, by analogy to (49), and using (117), we have

[TABLE]

Equivalently, after evaluating the integrals,

[TABLE]

After multiplying by $\left(x^{3}+1\right)^{2/3}$ the above Cauchy stress, and denoting $\lambda_{b}=b/B=\left(x^{3}+1\right)^{1/3}$, the required tensile dead load at the outer surface, $R=B$, in the reference configuration, takes the form

[TABLE]

The critical dead-load traction for the initiation of cavitation is

[TABLE]

and is the same as for the homogeneous sphere made entirely from the material found at its centre [67], as found also for the dynamic sphere.

To examine the post-cavitation behaviour, we take the first derivative of $P$, given by (131), with respect to $\lambda_{b}$,

[TABLE]

Letting $\lambda_{b}\to 1$ in (133), we obtain

[TABLE]

In Figure 19, we represent the scaled dead load, $P/C_{2}$, with $P$ given by (131) as a function of $\lambda_{b}$, and its derivative with respect to $\lambda_{b}$, for different values of the ratio $C_{1}/C_{2}$.

The following two types of cavitated solution are now possible:

• (i)

If $C_{1}>2C_{2}$, or equivalently, if $\mu(B)=C_{1}+C_{2}<3C_{1}/2$, then cavitation is stable;

• (ii)

If $C_{1}<2C_{2}$, or equivalently, if $\mu(B)=C_{1}+C_{2}>3C_{1}/2$, then snap cavitation occurs.

Thus, the probability of stable cavitation (and also of oscillatory motion for the dynamic sphere) is equal to (see the example in Figure 20)

[TABLE]

and that of snap cavitation (and also of non-oscillatory motion for the dynamic sphere) is

[TABLE]

An example of post-cavitation stochastic behaviour of the static inhomogeneous unit sphere (with $B=1)$ is shown in Figure 21, for the case when the shear modulus $\mu(R)$ is given by (109), with $\rho_{1}=405/R^{6}$ and $\rho_{2}=0.01\cdot R^{6}$ (see Figure 17). Thus, $\rho_{1}=405$ and $\rho_{2}=0.01$ at $R=B$, and by (109), $C_{1}=\underline{\mu}=\rho_{1}\rho_{2}=4.05$ and $C_{2}=\mu(B)-C_{1}$.

5.3 Alternative modelling formulations of radially inhomogeneous spheres

In the foregoing analysis, by following the approach of [71, 78], we have assumed that, for the inhomogeneous shear modulus $\mu(R)$, the mean value $\underline{\mu}(R)=\underline{\mu}$ is independent of $R$. Then, the shear modulus can be described by (109), where $C_{1}$ is a deterministic constant and $C_{2}$ is a (shifted) Gamma-distributed random variable with zero mean value. In this case, at the centre of the sphere, the shear modulus is a deterministic constant, while throughout the rest of the sphere, this modulus has constant mean value and non-constant variance [63].

An alternative model, where the mean value of the shear modulus varies throughout the sphere (see also [71], can be constructed by taking $C_{2}$ as a Gamma-distributed random variable (with non-zero mean value, i.e., $\underline{C}_{2}\neq 0$), while $C_{1}$ remains a deterministic constant. Then, by (109), for any fixed $R$, $\mu(R)$ is a shifted Gamma-distributed random variable, with the mean value and variance, respectively, equal to

[TABLE]

A more general modelling approach, where the shear modulus is characterised by a probability distribution at the centre of the sphere as well, is to assume that $C_{1}$ and $C_{2}$ are described by independent Gamma distributions, with hyperparameters $\rho^{(1)}_{1}$, $\rho^{(1)}_{2}$, and $\rho^{(2)}_{1}$, $\rho^{(2)}_{2}$, respectively. Then, for any fixed $R$, the distribution of $\mu(R)$ is a linear combination of two Gamma distributions (see Theorem 1 of [53]), while at the centre of the sphere, the shear modulus $\mu(0)=C_{1}$ is Gamma-distributed.

In all of these cases, the analysis follows the exposition of Sections 5.1 and 5.2, but the numerical results will differ in each case.

For the two alternative formulations, the probability of stable cavitation is equal to (see Figure 22)

[TABLE]

and that of snap cavitation is

[TABLE]

If $C_{1}>2C_{2}$, then the probability distribution of oscillatory motions occurring is given by (47) and that of non-oscillatory motions by (48), with $C_{1}$ instead of $\mu$ and $g(u;\rho^{(1)}_{1},\rho^{(1)}_{2})$ instead of $g(u;\rho_{1},\rho_{2})$ (see Figure 4).

For illustration, we choose a unit sphere with shear modulus $\mu(R)$ given by (109), where $C_{2}$ follows a Gamma distribution with shape and scale parameters $\rho^{(2)}_{1}=405$ and $\rho^{(2)}_{2}=0.05$, respectively, while $C_{1}$ is either the deterministic constant $C_{1}=4.05$, or a Gamma-distributed random variable, with hyperparameters $\rho^{(1)}_{1}=405$ and $\rho^{(1)}_{2}=0.01$. For the dynamic sphere, the associated function $H(x)$, defined by (119), and the velocity, given by (35), are represented in Figures 23 and 24, respectively. The corresponding post-cavitation behaviour of the static sphere is shown in Figure 25.

5.4 Non-oscillatory motion of a stochastic radially inhomogeneous sphere under impulse traction

We now consider the cavitation of radially inhomogeneous, incompressible spheres of stochastic hyperelastic material with the shear modulus of the form (109), such that, for any fixed $R$, both the mean value, $\underline{\mu}$, and variance, $\text{Var}[\mu]$, of $\mu=\mu(R)$ may depend on $R$.

Setting the initial conditions $x_{0}=x(0)=0$ and $\dot{x}_{0}=\dot{x}(0)=0$, we obtain equation (54), where

[TABLE]

Then, by substitution of (140) in (56) and (53), we have

[TABLE]

Equation (141) has a solution at $x_{1}=0$, while the second solution, $x_{2}>0$, is a root of

[TABLE]

To obtain the right-hand side of the above equation we have assumed that, in (141), $x$ is sufficiently small such that after expanding the respective functions to the second order in $x^{3}$,

[TABLE]

and

[TABLE]

Next, expanding the right-hand side of the equation (142) to first order in $x^{3}$ gives

[TABLE]

The critical dead load for the onset of cavitation is then

[TABLE]

and is the same as for the homogeneous sphere made entirely from the material found at the centre of the inhomogeneous sphere.

Noting that $2C_{1}+C_{2}=\mu(0)+\mu(B)>0$, the right-hand side in (145) decreases as $x$ increases, hence,

[TABLE]

or equivalently, by (53),

[TABLE]

However,

[TABLE]

and the sphere cannot oscillate.

5.5 Static deformation of a stochastic radially inhomogeneous sphere under impulse traction

For the static inhomogeneous sphere, when the cavity surface is traction-free, at the outer surface, the tensile traction takes the form

[TABLE]

Denoting $\lambda_{b}=b/B=\left(x^{3}+1\right)^{1/3}$, the required tensile load at the outer surface, $R=B$, can be written equivalently as follows,

[TABLE]

Thus the critical load for the initiation of cavitation is

[TABLE]

and is the same as for the homogeneous sphere made entirely from the material found at its centre, as found also for the dynamic sphere.

To examine the post-cavitation behaviour, we take the first derivative of $T$, given by (151), with respect to $\lambda_{b}$,

[TABLE]

Letting $\lambda_{b}\to 1$ in (153), we obtain

[TABLE]

Hence, only snap cavitation can occur.

In Figures 26 and 27, we show examples of post-cavitation stochastic behaviour of the static inhomogeneous unit sphere (with $B=1)$.

6 Conclusion

We have presented here a theoretical study of static and dynamic cavitation, and possible oscillatory motions, of homogeneous spheres of stochastic neo-Hookean material, of composite spheres with two concentric stochastic neo-Hookean phases, and of inhomogeneous spheres of neo-Hookean-like material with the material parameter varying continuously in the radial direction. For these materials, the elastic parameters are either spatially-independent random variables, or spatially-dependent random fields, described by Gamma probability density functions at a continuum level.

In our approach, the radially-symmetric dynamic deformation of each sphere was treated as a quasi-equilibrated motion. Then, the static deformation was examined by reducing the quasi-equilibrated motion to a static equilibrium state. Universally, we find that the critical load at which a nontrivial cavitated solution appears is the same as for the homogeneous sphere composed entirely of the material found at its centre. However, there are important differences in the post-cavitation behaviours under the different constitutive formulations and loading conditions. These different conditions give us a variety of different behaviours, specifically:

• •

For a neo-Hookean sphere, cavitation is stable and the motion is oscillatory under dead-load traction, while under the impulse traction, cavitation is unstable and the radial motion is non-oscillatory;

• •

In composite spheres with two concentric neo-Hookean phases subject to either dead-load or impulse traction, both stable and unstable cavitation can be obtained, depending on the material parameters of the two phases. However, oscillatory motion only occurs when cavitation is stable under dead-load traction;

• •

For radially inhomogeneous spheres under dead-load traction, both stable and unstable cavitation occur, and oscillatory motion is found when cavitation is stable under dead-load traction. Under impulse traction, cavitation is unstable and the motion is non-oscillatory.

Usually, in many application areas, the Gaussian (normal) distribution is assumed to describe the random variation that occurs in the available data. However, many material measurements show a skewed distribution where mean values are low, variances are large, and values cannot be negative (see, e.g., [19]). These are important features of the data that need to be incorporated in models. For small variability, the Gaussian and Gamma distributions both fit well to experimental data collected from material tests [19, 46]. This guides us to choose the distribution that is adequate also for cases exhibiting large variability, i.e., the Gamma distribution. If the data cannot be smaller than a certain bound different from zero, then a shift parameter can be incorporated to define a three-parameter distribution. To compare directly our stochastic results with the deterministic ones, we sampled from distributions where the parameters were set to have mean values corresponding to the deterministic system. Thus, the mean value of the distribution is guaranteed to converge to the expected value.

We found that, in contrast to the deterministic problem where a single critical value strictly separates the stable and unstable cases, for the stochastic problem, there is a probabilistic interval where the stable and unstable states compete, in the sense that both have a quantifiable chance to be obtained with a given probability. This means that solutions exist in the stochastic case where deterministically they should not. Generally, for the stochastic solution, the variance changes non-uniformly around the mean value, suggesting that, as the variance increases, the average value becomes less significant from the physical point of view. A natural question that is yet to be consider is: what is the effect of individual parameters on the variance in the stochastic solution?

For the finite dynamic deformations, attention was given to the periodic (oscillatory) motion and the time-dependent stresses, while taking into account the stochastic model parameters. In this case, the amplitude and period of the oscillation of the stochastic bodies are characterised by probability distributions, and there is a parameter interval where both the oscillatory and non-oscillatory motions can occur with a given probability. For the oscillating spheres, a cavity would form and expand until its radius reaches a maximum value, then contract again to zero periodically. An immediate extension of this could be to consider similar spheres where there is pressure at the cavity surface as well.

One could also apply the stochastic approach to the dynamic cavitation of spheres where the ‘base’ neo-Hookean material is replaced, for example, by a material with the following strain-energy function [48],

[TABLE]

where the exponents $m$ and $n$ are fixed, and the coefficients $\mu_{1}$ and $\mu_{2}$ are given, such that the shear modulus in infinitesimal deformation satisfies $\mu=\mu_{1}+\mu_{2}>0$. As is well know, in this case, assuming that the BE inequalities (1) hold, cavitation occurs if and only if [7, 12, 32, 48]

[TABLE]

Therefore, cavitation is expected in a neo-Hookean sphere (with $m=1$ and $n=0$), but not in a Mooney-Rivlin sphere (with $m=1$ and $n=-1$). The cases where $m\in\{-1/2,1\}$ and $n=0$ are discussed also as particular examples in [7], and the spheres where $m\in\{1/2,3/4,1,5/4\}$ and $n=0$ are analysed explicitly in [12]. For these well-behaved spheres, under uniform dead-load traction, the cavitated solution is always stable, and it makes little difference whether the model coefficients are stochastic or not, although, in the stochastic case, the solution appears somewhat ’hazy’ (albeit with the variance changing non-uniformly) around the mean value [48].

Similarly, for a sphere made of a material described by (155) subject to uniform radial impulse traction, if the BE inequalities and (156) are satisfied, then cavitation always occurs and making the model parameters stochastic only generates some dispersion around the usual deterministic solution. However, in this sphere, cavitation can also happen through a subcritical bifurcation where unstable (snap) cavitation takes place for values less than the critical dead load that causes cavitation. Hence, a sudden jump to a cavitated solution with a finite internal radius can be observed. For example, when the strain-energy density takes the form [48]

[TABLE]

both stable and unstable cavitation can occur, depending on the values of the material coefficients $\mu_{1}$ and $\mu_{2}$ (see Figure 3 of [48]). Then, in a stochastic sphere, the stable and unstable cases compete [48].

We further note that the theorem on quasi-equilibrated motion cannot be applied when the material is compressible [82, pp. 208-209]. As our analysis relies on the notion of quasi-equilibrated motion, the same approach cannot be employed for compressible spheres. However, while stochastic versions of compressible hyperelastic materials can also be constructed [76], there are not many theoretical results available on oscillatory motions of finitely deformed compressible hyperelastic spheres [2]. Nevertheless, elastodynamic problems involving more complex constrained, or unconstrained materials, can still be formulated, where the motion is not quasi-equilibrated, and the corresponding solution can then be tackled computationally. More generally, radial asymmetry, various behaviour laws, with or without damping, and thermal effects may also be of great theoretical and physical interest.

Acknowledgement.

The support by the Engineering and Physical Sciences Research Council of Great Britain under research grants EP/R020205/1 to Alain Goriely and EP/S028870/1 to L. Angela Mihai is gratefully acknowledged.