Likely cavitation in stochastic elasticity

We revisit the classic problem of elastic cavitation within the framework of stochastic elasticity. For the deterministic elastic problem, involving homogeneous isotropic incompressible hyperelastic spheres under radially symmetric tension, there is a critical dead-load traction at which cavitation can occur for some materials. In addition to the well-known case of stable cavitation post-bifurcation at the critical dead load, we show the existence of unstable snap cavitation for some isotropic materials satisfying Baker-Ericksen inequalities. For the stochastic problem, we derive the probability distribution of the deformations after bifurcation. In this case, we find that, due to the probabilistic nature of the material parameters, there is always a competition between the stable and unstable states. Therefore, at a critical load, stable or unstable cavitation occurs with a given probability, and there is also a probability that the cavity may form under smaller or greater loads than the expected critical value. We refer to these phenomena as `likely cavitation'. Moreover, we provide examples of homogeneous isotropic incompressible materials exhibiting stable or unstable cavitation together with their stochastic equivalent.


Introduction
Experiments carried out by Gent and Lindley in 1958 [13], on rubber cylinders, revealed that some materials can rupture under relatively small tensile dead loads by opening an internal cavity. Following this work, the theoretical analysis of Ball (1982) [6] provided an explanation for the formation of a spherical cavity at the centre of a sphere of isotropic hyperelastic incompressible material in radially symmetric tension under prescribed surface displacements or dead loads. There, the word 'cavitation' was used to describe such void-formation within a solid by analogy to the similar phenomenon observed in fluids. As cavitation in solids is an inherently nonlinear mechanical effect, not captured by the linear elasticity theory, many different studies have been devoted to the modelling of this effect within the finite elasticity framework. For instance: spheres of particular homogeneous isotropic incompressible materials were discussed in [8]; homogeneous anisotropic spheres with transverse isotropy about the radial direction were examined in [2,31,41]; concentric homogeneous spheres of different hyperelastic material were analysed in [18,42,48]; non-spherical cavities were investigated in [21]; cavities with non-zero pressure were presented in [49]; cavitation in an elastic membrane was studied in [56]; the homogenisation problem of nonlinear elastic materials was treated in [26,27]; growth-induced cavitation in nonlinearly elastic solids was explored in [15,30,40]. Recent experimental results on the onset, healing and growth of cavities in elastomers were reported in [43,44]. For many other results on cavitation in solids, we refer to the review articles [11,12,19], focusing on rubberlike materials, and the references therein.
The present work focuses on the phenomenon of cavitation contained within the theoretical context of finite elastostatics. Finite elasticity theory covers the simplest case where internal forces only depend on the current deformation of the material and not on its history, and is based on average data values. Within this framework, hyperelastic materials are the class of material models described by a strain-energy function characterised by a set of deterministic model parameters. In addition, for solid elastic materials, uncertainties in the observational data generally arise from the inherent variation in material properties and testing protocols [7,10,20,38]. In view of these uncertainties, recently, stochastic representations of isotropic incompressible hyperelastic materials characterised by a stochastic strain-energy function, for which the model parameters are random variables following standard probability laws, were proposed in [52], while compressible versions of these models were constructed in [53]. Ogden-type stochastic strain-energy functions were calibrated to experimental data for rubber and soft tissue materials in [35,54], and anisotropic stochastic models with the model parameters as spatially-dependent random field variables were calibrated to vascular tissue data in [55]. These models employ the maximum entropy principle for a discrete probability distribution introduced by Jaynes (1957) [22][23][24] and based on the notion of entropy (or uncertainty) defined by Shannon (1948) [46,47]. Such models can be useful for stochastic finite element implementations [3,4,16,17].
For stochastic hyperelastic models, the immediate question is: what is the influence of the random model parameters on the predicted nonlinear elastic responses? This question was previously considered by us in [36], for the stochastic Rivlin cube, and in [37], for the symmetric inflation of internally pressurised stochastic spherical shells and tubes. These idealised problems illustrate some important effects on the likely elastic responses of stochastic hyperelastic materials under large strains.
Here, we address this question by employing a similar approach as in [36,37] to revisit, in the context of stochastic elasticity, the cavitation problem of incompressible spheres of stochastic homogeneous isotropic hyperelastic materials under uniform radial tensile dead loads. Moreover, for all homogeneous isotropic hyperelastic models considered so far in the literature, cavitation appears as a supercritical bifurcation, where typically, after bifurcation, the cavity radius monotonically increases as the applied load increases (see, e.g., [8]). However, as we demonstrate here, the usual restriction that a material satisfies the Baker-Ericksen (BE) inequalities [5] is not sufficient to exclude the possibility of a subcritical bifurcation. In this case, one expects a snap cavitation for which there is a jump in the radius of the cavity immediately after bifurcation. Indeed, we obtain the general conditions under which a cavitation can appear through a supercritical or subcritical bifurcation and construct explicitly, for the first time, examples of isotropic incompressible hyperelastic models that exhibit snap cavitation. The stochastic version of these models are then explored. In this case, we find that, due to the probabilistic nature of the model parameters, supercritical or subcritical bifurcation occurs with a given probability, and there is also a probability that the cavity may form under smaller or greater loads that the expected critical value. We refer to these phenomena as 'likely cavitation'.
We begin, in Section 2, with a detailed presentation of the stochastic elastic framework. Then, in Section 3, for the stochastic sphere, after we review the elastic solution to the cavitation problem under uniformly applied tensile dead load, we recast the problem in the stochastic setting, and find the probabilistic solution. Concluding remarks are provided in Section 4.

Stochastic isotropic hyperelastic models
We recall that a homogeneous hyperelastic model is described by a strain-energy function W (F) that depends on the deformation gradient tensor, F, with respect to a fixed reference configuration, and is characterised by a set of deterministic model parameters [14,39,57]. In contrast, a stochastic homogeneous hyperelastic model is defined by a stochastic strain-energy function, for which the model parameters are random variables that satisfy standard probability laws [35,[52][53][54]. In this case, each model parameter is described in terms of its mean value and its variance, which contains information about the range of values about the mean value. While it is rarely possible if ever to obtain complete information about a random quantity in an elastic sample of material, the partial information provided by the mean value and the variance is the most commonly used in many practical applications [9,20,29]. Here, we combine finite elasticity and information theory, and rely on the following general hypotheses [35][36][37]: (A1) Material objectivity: The principle of material objectivity (frame indifference) states that constitutive equations must be invariant under changes of frame of reference. It requires that the scalar strain-energy function, W = W (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 (R T F) = W (F), where R ∈ SO(3) is a proper orthogonal tensor (rotation). Material objectivity is guaranteed by considering strain-energy functions defined in terms of invariants.
(A3) Baker-Ericksen inequalities: In addition to the fundamental principles of objectivity and material symmetry, in order for the behaviour of a hyperelastic material to be physically realistic, there are some universally accepted constraints on the constitutive equations. Specifically, for a hyperelastic body, the Baker-Ericksen (BE) inequalities, which state that the greater principal (Cauchy) stress occurs in the direction of the greater principal stretch, are [5]: where {λ i } i=1,2,3 and {T i } i=1,2,3 denote the principal stretches and the principal Cauchy stresses, respectively, and "≥" replaces the strict inequality ">" if any two principal stretches are equal. The BE inequalities (1) take the equivalent form [5,28] where the strict inequality ">" is replaced by "≥" if any two principal stretches are equal.
(A4) Finite mean and variance for the random shear modulus: We assume that, for any given finite deformation, the random shear modulus, µ, and its inverse, 1/µ, are second-order random variables, i.e., they have finite mean value and finite variance [52][53][54].
While (A4) contains physically realistic expectations on the random shear modulus, which will be drawn from a probability distribution, assumptions (A1)-(A3) are well-known principles in isotropic finite elasticity [14,39,57]. Specifically, we focus our attention on homogeneous incompressible hyperelastic materials characterised by the following stochastic strain-energy function [35,52,54], where m and n are deterministic constants, and µ 1 and µ 2 are random variables following given probability distributions. In the deterministic elastic case, µ 1 , µ 2 , m and n are constants, and the model contains, as special cases, the neo-Hookean model, the Mooney-Rivlin model, and the one-and two-term Ogden models. In both the deterministic elastic and stochastic cases, the shear modulus for infinitesimal deformations of these models is defined as µ = µ 1 + µ 2 [34,35]. Note that we could easily extend our description to include m and n as stochastic variables as well. However, increasing the complexity in this way is not relevant for the present discussion. Including additional sources of randomness is an avenue of future research. As it is well known, the deformation of an homogeneous isotropic hyperelastic material under uniaxial tension is a simple extension in the direction of the tensile force if and only if the BE inequalities hold [28]. Under these conditions, the shear modulus is positive, but the individual coefficients may be either positive or negative, allowing for some interesting nonlinear elastic effects to be captured (see [32][33][34]36] and the references therein). In particular, in the present paper, the initiation of either stable or unstable snap cavitation in a homogeneous isotropic sphere will be presented.
For the stochastic materials described by (3), condition (A4) is guaranteed by setting the following mathematical expectations [35,36,[52][53][54]: Then, under the constraints (4), the random shear modulus, µ, with mean value µ and standard deviation µ = Var[µ], defined as the square root of the variance, Var[µ], follows a Gamma probability distribution [50,51], with hyperparameters ρ 1 > 0 and ρ 2 > 0 satisfying The corresponding probability density function takes the form [1,25] where Γ : R * + → R is the complete Gamma function For technical convenience, we set a finite constant value b > −∞, such that µ i > b, i = 1, 2 (e.g., b = 0 if µ 1 > 0 and µ 2 > 0, but b is not unique in general), and introduce the auxiliary random variable [35] such that 0 < R 1 < 1. Consequently, we can equivalently express the random model parameters µ 1 and µ 2 as follows, It is reasonable to assume [35,[52][53][54]] in which case, the random variable R 1 , with mean value R 1 and variance Var[R 1 ], follows a standard Beta distribution [1,25], with hyperparameters ξ 1 > 0 and ξ 2 > 0 satisfying The associated probability density function is where B : R * Thus, for the random coefficients given by (9), the corresponding mean values take the form, and the variances and covariance are, respectively, It should be noted that the random variables µ and R 1 are independent, depending on parameters (ρ 1 , ρ 2 ) and (ζ 1 , ζ 2 ), respectively, which are derived by fitting distributions to given data. However, µ 1 and µ 2 are dependent variables as they both require (µ, R 1 ) to be defined. Explicit derivations of the probability distributions for the random parameters when stochastic isotropic hyperelastic models are calibrated to experimental data are presented in [35,54]. Our aim here is to analyse the radially symmetric finite deformations of a sphere of stochastic hyperelastic material defined by (3), under tension, when subject to prescribed surface dead loads applied uniformly in the radial direction. One can view the stochastic sphere as an ensemble (or population) of spheres, where each sphere has the same initial radius and is made from a homogeneous isotropic incompressible hyperelastic material, with the elastic parameters not known with certainty, but drawn from known probability distributions. Then, for every hyperelastic sphere, the finite elasticity theory applies. For the stochastic hyperelastic body, the question is: what is the probability distribution of stable radially symmetric deformation under a given surface dead load?

Incompressible spheres
In this section, we consider a sphere of stochastic incompressible hyperelastic material described by (3), subject to a radially symmetric deformation, caused by the sole action of a given radial tensile dead load. As for the deterministic elastic sphere [6], we obtain conditions on the constitutive law, such that, setting the internal pressure equal to zero, where the radius tends to zero, the required external dead load is finite, and therefore cavitation occurs. We further analyse the stability of the cavitated solution, and distinguish between supercritical cavitation, where the cavity radius monotonically increases as the dead load increases, and subcritical (snap) cavitation, with a sudden jump to a finite internal radius immediately after initiation. To the best of our knowledge, in the deterministic elastic case, the onset of snap cavitation in a homogeneous isotropic sphere has not been discussed before. Therefore, we start our analysis in the deterministic elastic context before extending it to the stochastic case.
For the stochastic sphere, the radially symmetric deformation takes the form where (R, Θ, Φ) and (r, θ, φ) are the spherical polar coordinates in the reference and current configuration, respectively, such that 0 ≤ R ≤ B, and g(R) ≥ 0 is to be determined. The corresponding deformation gradient is equal to F = diag (λ 1 , λ 2 , λ 3 ), with where λ 1 and λ 2 = λ 3 are the radial and hoop stretches, respectively, and dg/dR denotes the derivative of g with respect to R. By (19), hence, where c ≥ 0 is a constant to be calculated. If c > 0, then g(R) → c > 0 as R → 0 + , and a spherical cavity of radius c forms at the centre of the sphere, from zero initial radius (see Figure 1), otherwise the sphere remains undeformed. Assuming that the deformation (18) is due to a prescribed radial tensile dead load, applied uniformly on the sphere surface in the reference configuration, in the absence of body forces, the radial equation of equilibrium is or equivalently, where P = (P ij ) i,j=1,2,3 is the first Piola-Kirchhoff stress tensor. For an incompressible material, Denoting where λ = r/R = g(R)/R = (1 + c 3 /R 3 ) 1/3 > 1, we obtain Then, setting the internal pressure (at R → 0 + ) equal to zero, by (23) and (25), the external tension (at R = B) is equal to and the applied dead load, in the reference configuration, is where λ c and λ b represent the stretches at the centre and outer surface, respectively. The value of the required dead load, P 0 , for the onset of cavitation (bifurcation from the reference state) is obtained by taking λ c → ∞ and λ b = 1 + c 3 /B 3 1/3 → 1 as c → 0 + in (27), i.e., The BE inequalities (2) imply dW dλ hence, P 0 > 0. Then, if the critical dead load given by (28) is finite, cavitation takes place, else, the sphere remains undeformed. Before considering the stochastic setting, we briefly revisit the deterministic elastic case.
In particular, cavitation is found 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 special cases when m ∈ {−1/2, 1} and n = 0, are given as examples in [6], and when m ∈ {1/2, 3/4, 1, 5/4} and n = 0, the explicit critical loads are provided in [8]. When these bounds and the BE inequalities are satisfied, the critical pressure P 0 is finite and the problem is to find the behavior of the cavity in a neighborhood of this critical value.
In each of those previously studied cases (see, e.g., Figure 2 of [8]), cavitation forms from zero radius and then presents itself as a supercritical bifurcation with stable cavitation (i.e. the new bifurcated solution exists locally for values of P > P 0 , and the radius of the cavity monotonically increases with the applied load post-bifurcation). Another theoretical possibility is that the bifurcation could be subcritical (i.e., the cavitated solution exists locally for values less than P 0 and is unstable). One would then expect an unstable snap cavitation with a sudden jump to a cavitated solution with a finite internal radius. This subcritical behaviour of the homogeneous isotropic elastic sphere has not been explicitly demonstrated in the literature before. Here, we show that, depending on the model parameters, the family of materials (24) can exhibit both behaviours. General conditions for a given material to exhibit either a subcritical or supercritical bifurcation are provided in Appendix A.
As an example, we illustrate the variety of behaviours when m = 1 and n = −1/2, such that (24) takes the form In this case, under the deformation (18), the BE inequalities (2) are reduced to The inequality (36) implies that, when λ → 1, the shear modulus must be positive, i.e., µ = µ 1 +µ 2 > 0, while if λ → ∞, then µ 1 +2µ 2 > 0. Noting that the function of λ on the left-hand side is monotonically   increasing when µ 2 is positive, and decreasing if µ 2 is negative, and taking µ 1 > 0, the two limits imply that the BE inequalities are satisfied for all values of λ if For sufficiently small c/B, the corresponding dead-load traction, defined by (27), is equal to Then, the critical tensile dead load given by (28) takes the form As P 0 is positive in (39), we have 0 < µ 1 /µ < 8/3, which is guaranteed by (37). The question is now to find the possible behaviour of the cavity opening c as a function of P in a neighbourhood of P 0 . On differentiating (38) with respect to c/B, we obtain Hence, by Proposition A.1 given in Appendix A (with n = 3), when where "inf" denotes infimum, the bifurcation is supercritical and the radius of the cavity monotonically increases as the tensile dead load increases. However, if there exists c 0 > 0, such that then the bifurcation is subcritical and the required applied load starts to decrease at c = c 0 , where there is a sudden jump in the opening of cavity. In particular, if (42) holds for c 0 = 0, i.e., then (37) is valid while the cavitation becomes unstable. Thus, dP/d(c/B) → 0 as c → 0 + , and by Proposition A.1, the bifurcation at the critical load, P 0 , is supercritical (respectively, subcritical) if dP/d(c/B) > 0 (respectively, dP/d(c/B) < 0) for arbitrarily small c/B. Examples of both these behaviours are illustrated in Figures 2 and 3.

Stochastic elastic spheres
We now turn our attention to the stochastic model described by (3), with m = 1 and n = −1/2, and the other parameters drawn from probability distributions. In this case, recalling that µ follows a Gamma distribution g(u; ρ 1 , ρ 2 ), defined by (6), the probability distribution of stable cavitation is equal to and that of unstable cavitation is   For example, taking ρ 1 = 405 and ρ 2 = 0.01 (see Figure 4), the mean value of the shear modulus is µ = ρ 1 ρ 2 = 4.05, and the probability distributions given by equations (44)- (45) are illustrated numerically in Figure 5 (with blue lines for P 1 and red lines for P 2 ). In this case, if µ 1 = 5 < 5.4 = 4µ/3 say, then stable cavitation is expected, but there is also about 10% chance that unstable snap cavitation occurs. Similarly, when 4µ/3 = 5.4 < µ 1 = 5.8 < 8.1 = 2µ, unstable cavitation is expected, but there is also about 10% chance that the cavitation is stable. Stable and unstable cavitation of a stochastic sphere are illustrated numerically in Figure 6. Specifically: (a) In Figure 6(a), b = 0 in (8), and the random variable R 1 = µ 1 /µ is drawn from a Beta distribution with ξ 1 = 287 and ξ 2 = 36. In this case, µ 1 = 3.6 < 5.4 = 4µ/3, and stable cavitation, with supercritical bifurcation after the spherical cavity opens, is expected.
(b) In Figure 6(b), b = −3 in (8), and the random variable R 1 = (µ 1 + 3)/(µ + 6) draws its values from a Beta distribution with ξ 1 = 325 and ξ 2 = 10. Thus, 4µ/3 = 5.4 < µ 1 = 6.75 < 8.1 = 2µ, and unstable cavitation, with subcritical bifurcation after the spherical cavity forms, is expected. For the numerical examples shown in Figure 6 also, the critical dead load is P 0 = 4µ − 3µ 1 /2, as given by (39), with µ and µ 1 following probability distributions. In each case, the expectation is that the onset of cavitation occurs at the mean value P 0 = 4µ − 3µ 1 /2, found at the intersection of the dashed black line with the horizontal axis. However, there is a chance that cavity can form under smaller or greater critical loads that the expected load value, as shown by the coloured interval about the mean value along the horizontal axis.
To summarise, for a stochastic elastic sphere under uniform tensile dead load, we obtain the probabilities of stable or unstable cavitation, given that the material parameters are generated from known probability density functions. In the deterministic elastic case, there is a single critical parameter value that strictly separates the cases where the initiation of either stable or unstable cavitation occurs. By contrast, in the stochastic case, there is a probabilistic interval, containing the deterministic critical value, where there is always a competition between the stable and unstable states in the sense that both have a quantifiable chance to be found. For the onset of cavitation, there is also a probabilistic interval where a cavity may form, with a given probability, under smaller or greater loads that the expected critical value.

Conclusion
This work is motivated by the fact that a crucial part in assessing the physical properties of many solid materials is to quantify the uncertainties in their mechanical responses, which cannot be ignored. In particular, the idealised problem of the formation of a spherical cavity at the centre of a solid sphere illustrates some important effects on the likely elastic responses of stochastic hyperelastic materials under large strains.
For homogeneous isotropic incompressible spheres of stochastic hyperelastic material, subject to radial tensile dead loads applied uniformly on the sphere surface, we examined the possible radially symmetric deformations and determined which of these deformations are stable. Homogeneous stochastic hyperelastic material models satisfying certain theoretical assumptions were recently introduced to capture the dispersion in experimental data in addition to the traditional mean-data values [35,54].
For the deterministic elastic problem, where the model parameters are single-valued constants, non-trivial deformations, whereby a spherical cavity forms at the centre of the sphere, are possible for some classes of materials when the applied tensile dead loads are sufficiently large [6]. In some materials, cavitation is stable, in the sense that the cavity radius monotonically increases as the applied dead load increases [8]. Here, we showed that a sudden jump in the cavity opening, causing unstable snap cavitation, at the critical dead load can also occur in a homogeneous isotropic incompressible sphere, provided that the material satisfies Baker-Ericksen inequalities. If such a material could be found, a sphere made of this material would suddenly increase its volume at a critical load and show some form of hysteresis as the load is removed.
In the stochastic case, the probabilistic nature of the solution reflects the probability in the constitutive law, and bifurcation and stability can be quantified in terms of probabilities. By contrast to the deterministic elastic problem, where deterministic critical parameter values strictly separate the cases where either the stable or unstable cavitation occurs, for the stochastic problem, we obtained probabilistic intervals where both states have a quantifiable chance to exist. For the onset of cavitation, there is a probabilistic interval where the cavity may form, with a given probability, under smaller or greater loads that the expected critical value.
As a direct application of our approach, one could consider the cavitation of an inhomogeneous sphere made of concentric homogeneous spheres of different stochastic material, similar to the concentric homogeneous spheres of deterministic elastic material treated in [18] and [48]. Such composite spheres would require comparing both ensemble and spatial averages.

A Stability analysis
We provide a corrected version of Proposition 5.2 of [6] and its proof. In particular, we show that, in the deterministic elastic case, both subcritical and supercritical behaviours close to the bifurcation are possible, depending on the material.
These cases are illustrated, for the particular example of material presented in this paper, in Figures 2  and 3.
Note that the difference between this result and Proposition 5.2 of [6] comes from the (correct) minus sign between the two terms on the right-hand side of (48) (whereas a plus sign is found in the corresponding unlabelled expression appearing between Equations (5.25) and (5.26) of [6]).