Likely oscillatory motions of stochastic hyperelastic solids

Abstract

Stochastic homogeneous hyperelastic solids are characterized by strain-energy densities where the parameters are random variables defined by probability density functions. These models allow for the propagation of uncertainties from input data to output quantities of interest. To investigate the effect of probabilistic parameters on predicted mechanical responses, we study radial oscillations of cylindrical and spherical shells of stochastic incompressible isotropic hyperelastic material, formulated as quasi-equilibrated motions where the system is in equilibrium at every time instant. Additionally, we study finite shear oscillations of a cuboid, which are not quasi-equilibrated. We find that, for hyperelastic bodies of stochastic neo-Hookean or Mooney–Rivlin material, the amplitude and period of the oscillations follow probability distributions that can be characterized. Further, for cylindrical tubes and spherical shells, when an impulse surface traction is applied, there is a parameter interval where the oscillatory and non-oscillatory motions compete, in the sense that both have a chance to occur with a given probability. We refer to the dynamic evolution of these elastic systems, which exhibit inherent uncertainties due to the material properties, as ‘likely oscillatory motions’.

stochastic hyperelastic models, dynamic finite strain deformation, quasi-equilibrated motion, finite amplitude oscillations, incompressibility, applied probability

“Denominetur motus talis, qualis omni momento temporis |$t$| praebet configurationem capacem aequilibrii corporis iisdem viribus massalibus sollicitati, ‘motus quasi aequilibratus’. Generatim motus quasi aequilibratus non congruet legibus dynamicis et proinde motus verus corporis fieri non potest, manentibus iisdem viribus masalibus.” - Truesdell (1962)

1. Introduction

Motivated by numerous long-standing and modern engineering problems, oscillatory motions of cylindrical and spherical shells made of linear elastic material (Krauss, 1967; Love, 1888, 1944; Reissner, 1941) have generated a wide range of experimental, theoretical and computational studies (Alijani & Amabili, 2014; Amabili, 2008; Amabili & Païdoussis, 2003; Breslavsky & Amabili, 2018; Dong et al., 2018). In contrast, time-dependent finite oscillations of cylindrical tubes and spherical shells of nonlinear hyperelastic material, relevant to the modelling of physical responses in many biological and synthetic systems (Ahamed et al., 2018; Aranda-Iglesias et al., 2017; Destrade et al., 2011; Haas & Goldstein, 2015, 2019; Haslach & Humphrey, 2004; Kumar & DasGupta, 2013), have been less investigated, and much of the work in finite nonlinear elasticity has focused on the static stability of pressurized shells (Adkins & Rivlin, 1952; Biscari & Omati, 2010; Bucchi & Hearn, 2013a,b; Carroll, 1987; Fu et al., 2016; Goncalves et al., 2008; Goriely et al., 2006; Green & Shield, 1950; Mangan & Destrade, 2015; Müller & Struchtrup, 2002; Rivlin, 1949; Shield, 1972; Zamani & Pence, 2017), or on wave-type solutions in infinite media (Il’ichev & Fu, 2014; Pearce & Fu, 2010).

The governing equations for large amplitude oscillations of cylindrical tubes and spherical shells of homogeneous isotropic incompressible nonlinear hyperelastic material, formulated as special cases of quasi-equilibrated motions (Truesdell, 1962), were reviewed in Truesdell & Noll (2004). These are the class of motions for which the deformation field is circulation preserving, and at every time instant, the current configuration is a possible static configuration under the given forces. The free and forced axially symmetric radial oscillations of infinitely long, isotropic incompressible circular cylindrical tubes, with arbitrary wall thickness, were described for the first time in Knowles (1960, 1962). In Heng & Solecki (1963), Knowles & Jakub (1965) and Wang (1965), free and forced oscillations of spherical shells were derived analogously. For the combined radial–axial large amplitude oscillations of hyperelastic cylindrical tubes, in Shahinpoor (1973), the surface tractions necessary to maintain the periodic motions were discussed, and the results were applied to a tube sealed at both ends and filled with an incompressible fluid. The dynamic deformation of cylindrical tubes of Mooney–Rivlin material in finite amplitude radial oscillation was obtained in Shahinpoor & Nowinski (1971), Shahinpoor (1973) and Shahinpoor & Balakrishnan (1978). Oscillatory motion caused by the dynamic cavitation of a neo-Hookean sphere was considered in Chou-Wang & Horgan (1989). For a wide class of hyperelastic materials, both the static and dynamic cavitation of homogeneous spheres were analysed in Ball (1982). For a hyperelastic sphere of Mooney–Rivlin material, with a cavity, the solution to the nonlinear problem of large amplitude oscillations was computed numerically in Balakrishnan & Shahinpoor (1978). Theoretical and experimental studies of cylindrical and spherical shells of rubberlike material under external pressure were presented in Wang & Ertepinar (1972). In Calderer (1983), 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 studied theoretically. The finite longitudinal, or ‘telescopic’, oscillations of infinitely long cylindrical tubes were investigated in Nowinski & Schultz (1964). In Nowinski (1966), the oscillatory motions of cylindrical and prismatic bodies of incompressible hyperelastic material under dynamic finite shear deformation were analysed. Other dynamic shear deformations were considered in Wang (1969), where it was emphasized that such shear motions were not quasi-equilibrated. In Huilgol (1967), the dynamic problem of axially symmetric oscillations of cylindrical tubes of transversely isotropic incompressible material, with radial transverse isotropy, was treated. The dynamic deformation of a longitudinally anisotropic thin-walled cylindrical tube under radial oscillations was obtained in Shahinpoor (1974). In Ertepinar & Akay (1976), radial oscillations of non-homogeneous thick-walled cylindrical and spherical shells of neo-Hookean material, with a material constant varying continuously along the radial direction, were studied. In Akyüz & Ertepinar (1998), for pressurized homogeneous isotropic compressible hyperelastic tubes of arbitrary wall thickness under uniform radial dead-load traction, the stability of the finitely deformed state and small radial vibrations about this state were treated, using the theory of small deformations superposed on large elastic deformations, while the governing equations were solved numerically. In Verron et al. (1999), the dynamic inflation of hyperelastic spherical membranes of Mooney–Rivlin material subjected to a uniform step pressure was studied, and the absence of damping in these models 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, for many systems under slowly increasing pressure, strong damping would generally preclude oscillations (De Pascalis et al., 2018). More recently, the dynamic response of incompressible hyperelastic cylindrical and spherical shells subjected to periodic loading was discussed in Ren (2008, 2009). Radial oscillations of cylindrical tubes and spherical shells of neo-Hookean (Treloar, 1944), Mooney–Rivlin (Mooney, 1940; Rivlin, 1948) and Gent (1996) hyperelastic materials were analysed in Beatty (2007, 2011), where it was inferred that, in general, both the amplitude and period of oscillations decrease when the stiffness of the material increases. The influence of material constitutive law on the dynamic behaviour of cylindrical and spherical shells was also examined in Aranda-Iglesias et al. (2018, 2015), Rodriíguez-Martiínez et al. (2015), and Yuan et al. (2008), where the results for Yeoh (1993) and Mooney–Rivlin material models were compared. In Breslavsky et al. (2016), the static and dynamic behaviour of circular cylindrical shells of homogeneous isotropic incompressible hyperelastic material modelling arterial walls were considered. In (Soares et al., 2019), 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 its inner surface, was studied, and a parametric analysis of the influence of the material constants was presented.

For the assessment and prediction of the mechanical responses of engineered and natural materials, additional challenges arise from the uncertainties in their elastic properties inferred from sparse and approximate observational data (Ghanem et al., 2017; Hughes & Hase, 2010; Kaminski & Lauke, 2018; Oden, 2018; Ostoja-Starzewski, 2007; Sullivan, 2015). For these 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 significantly improve assessment and predictions. In response to this challenge, stochastic elasticity is a fast-developing field that combines nonlinear elasticity and stochastic theories in order to significantly improve model predictions by accounting for uncertainties in the mechanical responses of materials. Within this framework, stochastic hyperelastic materials are advanced phenomenological models described by a strain-energy density where the parameters are characterized by probability density functions, as constructed in Staber & Guilleminot (2015, 2016, 2017, 2018), Staber et al. (2019) and Mihai et al. (2018c). These models rely on the notion of entropy (or uncertainty) (Shannon, 1948; Soni & Goodman, 2017) and on the maximum entropy principle for a discrete probability distribution (Jaynes, 1957a,b, 2003), and allow for the propagation of uncertainties from input data to output quantities of interest (Soize, 2013). They are also suitable for incorporation into Bayesian methodologies (Bayes, 1763; McGrayne, 2012) for models selection or updates (Mihai et al., 2018c; Oden, 2018; Robert, 2007).

To study the effect of probabilistic model parameters on predicted mechanical responses, in Mihai et al. (2018a,b, 2019a,b), for different bodies with simple geometries at finite strain deformations, it was shown explicitly that, in contrast to the deterministic elastic 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 always compete, in the sense that both have a quantifiable chance to be found. In addition, revisiting these problems from a novel perspective offered fresh opportunities for gaining new insights into the fundamental elastic solutions, and correcting some inconsistencies found in the previous works. Specific case studies, so far, include the cavitation of a sphere under uniform tensile dead load (Mihai et al., 2018a), the inflation of pressurized spherical and cylindrical shells (Mihai et al., 2018b), the classical problems of the Rivlin cube (Mihai et al., 2019a) and the rotation and perversion of anisotropic hyperelastic cylindrical tubes (Mihai et al., 2019b).

In this paper, we extend the stochastic framework developed in Mihai et al. (2018a,b, 2019a,b) to study radial oscillations of cylindrical and spherical shells of stochastic incompressible isotropic hyperelastic material formulated as quasi-equilibrated motions. For these motions, the system is in equilibrium at every time instant. We consider also finite shear oscillations of a cuboid, which are not quasi-equilibrated. We find that, for hyperelastic bodies of stochastic neo-Hookean or Mooney–Rivlin material, the amplitude and period of the oscillations follow probability distributions that can be fully characterized. Further, for cylindrical tubes and spherical shells, when an impulse surface traction is applied, there is a parameter interval where the oscillatory and non-oscillatory motions compete in the sense that both have a chance to occur with a given probability. We refer to the dynamic evolution of these elastic systems, which exhibit inherent uncertainties due to the material properties, as ‘likely oscillatory motions’. Section 2 provides a summary of the stochastic elasticity prerequisites. Section 3 is devoted to the oscillatory motions of a stochastic hyperelastic cuboid under dynamic generalized shear. This is followed, in Sections 4 and 5, by the radial oscillatory motions of stochastic cylindrical and spherical shells with bounded wall thickness, respectively. The limiting cases of thin- and infinitely thick-walled structures are also discussed. Some less straightforward calculations, inherent for these problems, are deferred to Appendix  A. Concluding remarks are drawn in Section 6.

2. Prerequisites

In this section, we recall the notion of (universal) quasi-equilibrated motion in finite elasticity, introduced in Truesdell (1962) and reviewed in Truesdell & Noll (2004) and summarize the stochastic finite elasticity framework developed in Mihai et al. (2018c) and applied to various static stability problems in Mihai et al. (2018a,b, 2019a).

2.1 Quasi-equilibrated motion

For the large strain time-dependent behaviour of an elastic solid, Cauchy’s laws of motion (balance laws of linear and angular momentum) are governed by the following Eulerian field equations (Truesdell & Noll, 2004, p. 40):

$$\begin{eqnarray}&&\rho\ddot{\textbf{x}}=\mathrm{div}\ \textbf{T}+\rho\textbf{b},\end{eqnarray}$$

(1)

$$\begin{eqnarray}&&\textbf{T}=\textbf{T}^{T},\end{eqnarray}$$

(2)

where |$\textbf{x}=\chi (\textbf{X},t)$| is the motion of the elastic solid, |$\rho$| is the material density, which is assumed constant, |$\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 generalize known static solutions to dynamical forms, using the so-called quasi-equilibrated motion, which is defined as follows.

 

Definition 2.1

(Truesdell & Noll, 2004, 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

(Truesdell & Noll, 2004, p. 208) 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.,

$$\begin{equation}\ddot{\textbf{x}}=-\mathrm{grad}\ \xi.\end{equation}$$

(3)

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

$$\begin{equation}\mathrm{curl}\ \ddot{\textbf{x}}=\textbf{0}.\end{equation}$$

(4)

Then, the Cauchy stress tensor takes the form

$$\begin{equation}\boldsymbol{T}=-\rho\xi\boldsymbol{I}+\boldsymbol{T}^{(0)},\end{equation}$$

(5)

where |$\boldsymbol{T}^{(0)}$| is the Cauchy stress for the equilibrium state at time |$t$| and |$\boldsymbol{I}=\textrm{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$|⁠.

 

Proof.

The Cauchy stress |$\textbf{T}^{(0)}$| for the equilibrium state under the body force |$\textbf{b}=\textbf{b}(\textbf{x},t)$| at time |$t$| satisfies

$$\begin{equation}-\mathrm{div}\ \boldsymbol{T}^{(0)}=\rho\textbf{b}.\end{equation}$$

(6)

First, we assume that the motion |$\textbf{x}=\chi (\boldsymbol{X},t)$| is quasi-equilibrated under the body force |$\textbf{b}=\textbf{b}(\textbf{x},t)$| and deduce that there is a single-valued function |$\xi$|⁠, such that (3) holds. Indeed, if the motion is quasi-equilibrated, then Definition 2.1 implies that, at any fixed time-instant |$t$|⁠, the Cauchy stress takes the form (5), where |$\xi =\xi (t)$| is a single-valued function of |$t$|⁠. Substituting (5) in (1) gives

$$\begin{equation}\rho\ddot{\textbf{x}}=-\rho\ \mathrm{grad}\ \xi+\mathrm{div}\ \boldsymbol{T}^{(0)}+\rho\textbf{b}.\end{equation}$$

(7)

Then, (3) follows from (6) and (7).

Conversely, if (3) holds, with |$\xi$| a single-valued function, then, substitution of (3) and (6) in (1) gives

$$\begin{equation}-\rho\ \mathrm{grad}\ \xi=\mathrm{div}\ \left(\boldsymbol{T}-\boldsymbol{T}^{(0)}\right),\end{equation}$$

(8)

at any time-instant |$t$|⁠. From (8), it follows that the Cauchy stress |$\boldsymbol{T}$| takes the form (5). Hence, the motion is quasi-equilibrated according to Definition 2.1.

Theorem 2.2 may only be applicable to specific quasi-equilibrated motions of specific materials. Nevertheless, by the above theorem, for a quasi-equilibrated motion to be dynamically possible under a given body force in all elastic materials, it is necessary that, 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. Some examples are the homogeneous motions that are possible in all homogeneous incompressible materials, and also those considered by us in Sections 4 and 5 (see also Truesdell & Noll, 2004, p. 209).

2.2 Stochastic isotropic incompressible hyperelastic models

A stochastic homogeneous hyperelastic model is defined by a stochastic strain-energy function, for which the model parameters are random variables, drawn from probability distributions (Mihai et al., 2018c; Staber & Guilleminot, 2015, 2016, 2017). In this case, each model parameter is usually described in terms of its mean value and its variance, which contains information about the range of values about the mean value (Brewick & Teferra, 2018; Caylak et al., 2018; Hughes & Hase, 2010; McCoy, 1973; Nörenberg & Mahnken, 2015). Here, we combine finite elasticity (Goriely, 2017; Ogden, 1997; Truesdell & Noll, 2004) and probability theory (Grimmett & Stirzaker, 2001; Jaynes, 2003) and rely on the following general assumptions (Mihai et al., 2018a,b,c, 2019a):

• (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 |$\textbf{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 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 |$\textbf{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, are (Baker & Ericksen, 1954)

  $$\begin{equation}\left(T_{i}-T_{j}\right)\left(\lambda_{i}-\lambda_{j}\right)>0\quad \mbox{if}\quad \lambda_{i}\neq\lambda_{j},\quad i,j=1,2,3,\end{equation}$$

  (9)

  where |$\{\lambda _{i}\}_{i=1,2,3}$| and |$\{T_{i}\}_{i=1,2,3}$| denote the principal stretches and the principal Cauchy stresses, respectively. The BE inequalities (9) take the equivalent form

  $$\begin{equation}\left(\lambda_{1}\frac{\partial\mathcal{W}}{\partial\lambda_{1}}-\lambda_{2}\frac{\partial\mathcal{W}}{\partial\lambda_{2}}\right)\left(\lambda_{1}-\lambda_{2}\right)>0\quad \mbox{if}\quad \lambda_{i}\neq\lambda_{j},\quad i,j=1,2,3.\end{equation}$$

  (10)

  In (9) and (10), the strict inequality ‘>’ is replaced by ‘|$\geq$|’ if any two principal stretches are equal.

• (A4) Finite mean and variance of the random shear modulus, i.e., at any given deformation, the random shear modulus, |$\mu$|⁠, and its inverse, |$1/\mu$|⁠, are second-order random variables (Staber & Guilleminot, 2015, 2016, 2017).

Assumptions (A1)–(A3) are well-known principles in isotropic finite elasticity (Goriely, 2017; Ogden, 1997; Truesdell & Noll, 2004). In particular, regarding (A3), we recall that, for a homogeneous hyperelastic body under uniaxial tension, the deformation is a simple extension in the direction of the tensile force if and only if the BE inequalities (9) hold (Marzano, 1983). Another important deformation is that of simple shear superposed on axial stretch, defined by

$$\begin{equation}x_1=\alpha X_1+k\frac{X_2}{\alpha^2},\qquad x_2=\frac{X_2}{\alpha^2},\qquad x_3=\alpha X_3,\end{equation}$$

(11)

where |$(X_1,X_2,X_3)$| and |$(x_1,x_2,x_3)$| are the Cartesian coordinates for the reference (Lagrangian) and the current (Eulerian) configuration, respectively, and |$k>0$| and |$\alpha>0$| are positive constants representing the shear parameter and the axial stretch (⁠|$0<\alpha <1$| for axial tension and |$\alpha>1$| for axial compression), respectively. For this deformation, the principal stretches, |$\{\lambda _{i}\}_{i=1,2,3}$|⁠, satisfy

$$\begin{equation}\begin{split} \lambda_{1}^2&=\frac{\alpha^6+k^2+1+\sqrt{\left(\alpha^6+k^2+1\right)^2-4\alpha^6}}{2\alpha^4}, \\ \lambda_{2}^2&=\frac{\alpha^6+k^2+1-\sqrt{\left(\alpha^6+k^2+1\right)^2-4\alpha^6}}{2\alpha^4}, \\ \lambda_{3}^2&=\alpha^2.\end{split}\end{equation}$$

(12)

Then assuming that the material is incompressible, the associated principal Cauchy stresses take the form

$$\begin{equation}T_{i}=\lambda_{i}\frac{\partial\mathcal{W}}{\partial{\lambda_{i}}}-p,\qquad i=1,2,3,\end{equation}$$

(13)

where |$p$| is the Lagrange multiplier for the incompressibility constraint. In this case, if the BE inequalities (9) hold, then the nonlinear shear modulus defined by Mihai & Goriely (2017) and Mihai et al. (2018c)

$$\begin{equation}\widetilde{\mu}=\frac{T_{1}-T_{2}}{\lambda_{1}^2-\lambda_{2}^2}\end{equation}$$

(14)

is positive, i.e., |$\widetilde{\mu }>0$|⁠, for all |$k>0$| and |$\alpha>0$|⁠. In the linear elastic limit, where |$k\to 0$| and |$\alpha \to 1$|⁠, the nonlinear shear modulus given by (14) converges to the classical shear modulus under infinitesimal deformation, i.e.,

$$\begin{equation}\lim_{a\to1,k\to0}\widetilde{\mu}=\mu.\end{equation}$$

(15)

Assumption (A4) then contains physically realistic expectations on the (positive) random shear modulus |$\mu>0$|⁠, which will be characterized by a suitable probability density function.

In the next sections, we analyse the dynamic generalized shear deformation of a cuboid and the radially symmetric motion of cylindrical tube and spherical shells of stochastic isotropic incompressible hyperelastic material. One can regard a stochastic hyperelastic body as an ensemble of bodies with the same geometry where each individual body is made from a homogeneous isotropic incompressible hyperelastic material, with the elastic parameters drawn from probability distributions. Then for the individual hyperelastic bodies, the finite elasticity theory applies.

Throughout this paper, we confine our attention to a class of stochastic homogeneous incompressible hyperelastic materials described by the Mooney–Rivlin-like constitutive law (Mihai et al., 2018c; Staber & Guilleminot, 2015),

$$\begin{equation}\mathcal{W}(\lambda_{1},\lambda_{2},\lambda_{3})=\frac{\mu_{1}}{2}\left(\lambda_{1}^2+\lambda_{2}^2+\lambda_{3}^2-3\right) +\frac{\mu_{2}}{2}\left(\lambda_{1}^{-2}+\lambda_{2}^{-2}+\lambda_{3}^{-2}-3\right),\end{equation}$$

(16)

where |$\mu _{1}$| and |$\mu _{2}$| are random variables. The non-deterministic model (16) reduces to a stochastic neo-Hookean model if |$\mu _{2}=0$|⁠. If the random parameters |$\mu _{1}$| and |$\mu _{2}$| are replaced by their respective mean values, |$\underline{\mu }_{1}$| and |$\underline{\mu }_{2}$|⁠, then the resulting mean hyperelastic model coincides with the usual deterministic one.

For the stochastic material model described by (16), the shear modulus in infinitesimal deformation is defined as |$\mu =\mu _{1}+\mu _{2}$|⁠. For this modulus, we set the following mathematical constraints, which ensure that the assumption (A4), made in Section 2 is satisfied (Mihai et al., 2018c),

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

(17)

i.e., the mean value |$\underline{\mu }$| of the shear modulus |$\mu$| is fixed and greater than zero, and the mean value of |$\log \ \mu$| is fixed and finite. It follows that |$\mu$| and |$1/\mu$| are second-order random variables, i.e., they have finite mean and finite variance (Soize, 2000, 2001). Under the constraints (??), |$\mu$| follows a Gamma probability distribution with hyperparameters |$\rho _{1}>0$| and |$\rho _{2}>0$|⁠, such that

$$\begin{equation}\underline{\mu}=\rho_{1}\rho_{2},\qquad \textrm{Var}[\mu]=\rho_{1}\rho_{2}^2,\end{equation}$$

(18)

where |$\underline{\mu }$| is the mean value, |$\textrm{Var}[\mu ]$| is the variance and |$\|\mu \|=\sqrt{\textrm{Var}[\mu ]}$| is the standard deviation of |$\mu$|⁠. The corresponding probability density function takes the form (Abramowitz & Stegun, 1964; Johnson et al., 1994)

$$\begin{equation}g(\mu;\rho_{1},\rho_{2})=\frac{\mu^{\rho_{1}-1}e^{-\mu/\rho_{2}}}{\rho_{2}^{\rho_{1}}\varGamma(\rho_{1})},\qquad\mbox{for}\ \mu>0\ \mbox{and}\ \rho_{1}, \rho_{2}>0,\end{equation}$$

(19)

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

$$\begin{equation}\varGamma(z)=\int_{0}^{+\infty}t^{z-1}e^{-t}\,\textrm{d}t.\end{equation}$$

(20)

When |$\mu _{1}>0$| and |$\mu _{2}>0$|⁠, we can define the auxiliary random variable (Mihai et al., 2018c)

$$\begin{equation}R_{1}=\frac{\mu_{1}}{\mu},\end{equation}$$

(21)

such that |$0<R_{1}<1$|⁠. Then

$$\begin{equation}\mu_{1}=\mu R_{1},\qquad \mu_{2}=\mu-\mu_{1}=\mu(1-R_{1}).\end{equation}$$

(22)

Setting the realistic constraints (Mihai et al., 2018c; Staber & Guilleminot, 2015, 2016, 2017),

$$\begin{eqnarray}\begin{cases} E\left[\log\ R_{1}\right]=\nu_{1},& \mbox{such that {$|\nu_{1}|<+\infty$}},\\ E\left[\log(1-R_{1})\right]=\nu_{2},& \mbox{such that {$|\nu_{2}|<+\infty$}},\end{cases}\end{eqnarray}$$

(23)

we obtain that |$R_{1}$| follows a standard Beta distribution (Abramowitz & Stegun, 1964; Johnson et al., 1994), with hyperparameters |$\xi _{1}>0$| and |$\xi _{2}>0$| satisfying

$$\begin{equation}\underline{R}_{1}=\frac{\xi_{1}}{\xi_{1}+\xi_{2}},\qquad \textrm{Var}[R_{1}]=\frac{\xi_{1}\xi_{2}}{\left(\xi_{1}+\xi_{2}\right)^2\left(\xi_{1}+\xi_{2}+1\right)},\end{equation}$$

(24)

where |$\underline{R}_{1}$| is the mean value, |$\textrm{Var}[R_{1}]$| is the variance and |$\|R_{1}\|=\sqrt{\textrm{Var}[R_{1}]}$| is the standard deviation of |$R_{1}$|⁠. The corresponding probability density function takes the form

$$\begin{equation}\beta(r;\xi_{1},\xi_{2})=\frac{r^{\xi_{1}-1}(1-r)^{\xi_{2}-1}}{B(\xi_{1},\xi_{2})},\qquad \qquad\mbox{for}\ r\in(0,1)\ \mbox{and}\ \xi_{1}, \xi_{2}>0,\end{equation}$$

(25)

where |$B:\mathbb{R}^{*}_{+}\times \mathbb{R}^{*}_{+}\to \mathbb{R}$| is the Beta function

$$\begin{equation}B(x,y)=\int_{0}^{1}t^{x-1}(1-t)^{y-1}\,\mathrm{d}t.\end{equation}$$

(26)

Thus, for the random coefficients given by (22), the corresponding mean values take the form,

$$\begin{equation}\underline{\mu}_{1}=\underline{\mu}\underline{R}_{1}, \qquad \underline{\mu}_{2}=\underline{\mu}-\underline{\mu}_{1}=\underline{\mu}(1-\underline{R}_{1}),\end{equation}$$

(27)

and the variances and covariance are, respectively,

$$\begin{eqnarray}&&\textrm{Var}\left[\mu_{1}\right] =(\underline{\mu})^2\textrm{Var}[R_{1}]+(\underline{R}_{1})^2\textrm{Var}[\mu]+\textrm{Var}[\mu]\textrm{Var}[R_{1}],\end{eqnarray}$$

(28)

$$\begin{eqnarray}&&\textrm{Var}\left[\mu_{2}\right] =(\underline{\mu})^2\textrm{Var}[R_{1}]+(1-\underline{R}_{1})^2\textrm{Var}[\mu]+\textrm{Var}[\mu]\textrm{Var}[R_{1}],\end{eqnarray}$$

(29)

$$\begin{eqnarray}&&\textrm{Cov}[\mu_{1},\mu_{2}]=\frac{1}{2}\left(\textrm{Var}[\mu]-\textrm{Var}[\mu_{1}]-\textrm{Var}[\mu_{2}]\right).\end{eqnarray}$$

(30)

Note that the random variables |$\mu$| and |$R_1$| are independent, depending on parameters |$(\rho _1,\rho _2)$| and |$(\zeta _1,\zeta _2)$|⁠, respectively, whereas |$\mu _1$| and |$\mu _2$| are codependent variables as they both require |$(\mu ,R_1)$| to be defined.

For the numerical illustration of our subsequent results, throughout this paper, we assume that the random shear modulus |$\mu$| follows the Gamma distribution represented in Fig. 1, where the shape and scale parameters are |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠, respectively (Mihai et al., 2018b). Different simulations were created 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.

Fig. 1.

Example of Gamma distribution with hyperparameters |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠.

Note also that the Gamma distribution represented in Fig. 1, for which |$\rho _1$| is large compared to |$\rho _2$|⁠, appears to be approximately a normal distribution. However, despite known convergence results and the qualitative agreement between the two density functions for large values of the mean (see Mihai et al., 2018b for a detailed discussion), in general, the normal distribution cannot be used to model material parameters. This is due to the fact that the normal 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 the experiments. Then the maximum entropy principle allows for the explicit construction of their probability laws, given the available information. Explicit derivations of probability distributions for the elastic parameters of stochastic homogeneous isotropic hyperelastic models, calibrated to available experimental data, are presented in Staber & Guilleminot (2017) and Mihai et al. (2018c).

3. Generalized shear motion of a stochastic hyperelastic cuboid

First, we consider a stochastic hyperelastic cuboid subject to dynamic generalized shear.

3.1 Dynamic generalized shear

The generalized shear motion of an elastic body is described by Destrade et al. (2011)

$$\begin{equation}x=\frac{X}{\sqrt{\alpha}},\qquad y=\frac{Y}{\sqrt{\alpha}},\qquad z=\alpha Z+u\left(X,Y,t\right),\end{equation}$$

(31)

where |$(X,Y,Z)$| and |$(x,y,z)$| are the Cartesian coordinates for the reference (Lagrangian, material) and current (Eulerian, spatial) configuration, respectively, |$\alpha>0$| is a given constant, and |$u=z-Z$|⁠, representing the displacement in the third direction, is a time-dependent function to be determined. Here we assume that the edges of the cuboid are aligned with the directions of the Cartesian axes in the undeformed state (see Fig. 2).

Fig. 2.

Schematic of generalized shear of a cuboid, showing the reference state (left) and the deformed state (right), respectively.

By the governing equations (31), the condition (4) is valid for |$\textbf{x}=(x,y,z)^{T}$| if and only if

$$\begin{equation}\textbf{0}=\mathrm{curl}\ \ddot{\textbf{x}}= \left[\begin{array}{c} \partial\ddot{z}/\partial y-\partial\ddot{y}/\partial z\\ \partial\ddot{x}/\partial z-\partial\ddot{z}/\partial x\\ \partial\ddot{y}/\partial x-\partial\ddot{x}/\partial y\end{array} \right]= \left[\begin{array}{c} \partial\ddot{u}/\partial Y\\ -\partial\ddot{u}/\partial X\\ 0\end{array} \right].\end{equation}$$

(32)

This condition imposes very strict constraints on the motion. Yet, we will see that even though the generalized shear motion (31) is not quasi-equilibrated, exact solutions are still available, although these solutions are not universal (Nowinski, 1966; Wang, 1969).

For the deformation (31), the gradient tensor is equal to

$$\begin{equation*}\textbf{F} =\left[\begin{array}{ccc} 1/\sqrt{\alpha} & 0 & 0\\ 0 & 1/\sqrt{\alpha} & 0\\ u_{X} & u_{Y} & \alpha\end{array} \right],\end{equation*}$$

where |$u_{X}$| and |$u_{Y}$| denote the partial first derivatives of |$u$| with respect to |$X$| and |$Y$|⁠, respectively. The corresponding left Cauchy–Green tensor is

$$\begin{equation}\boldsymbol{B}=\boldsymbol{F}\boldsymbol{F}^{T}=\left[\begin{array}{ccc} 1/\alpha & 0 & u_{X}/\sqrt{\alpha} \\ 0 & 1/\alpha & u_{Y}/\sqrt{\alpha} \\ u_{X}/\sqrt{\alpha} & u_{Y}/\sqrt{\alpha} & u_{X}^2+ u_{Y}^2+\alpha^2\end{array} \right]\end{equation}$$

(33)

and has the principal invariants

$$\begin{equation}\begin{split} I_{1}=&\mathrm{tr}\ (\boldsymbol{B})=u_{X}^2+u_{Y}^2+\frac{2}{\alpha}+\alpha^2,\\ I_{2}=&\frac{1}{2}\left[\left(\mathrm{tr}\,\boldsymbol{B}\right)^{2}-\mathrm{tr}\left(\boldsymbol{B}^{2}\right)\right]=\frac{u_{X}^2}{\alpha}+\frac{u_{Y}^2}{\alpha}+\frac{1}{\alpha^2}+2\alpha,\\ I_{3}=&\det\boldsymbol{B}=1.\end{split}\end{equation}$$

(34)

The associated Cauchy stress tensor takes the form (Green & Adkins, 1970, pp. 87–91)

$$\begin{equation}\boldsymbol{T}=-p\boldsymbol{I}+\beta_{1}\boldsymbol{B}+\beta_{-1}\boldsymbol{B}^{-1},\end{equation}$$

(35)

where |$p$| is the Lagrange multiplier for the incompressibility constraint (⁠|$I_{3}=1$|⁠), and

$$\begin{equation}\beta_{1}={2}\frac{\partial W}{\partial I_{1}},\qquad \beta_{-1}=-2\frac{\partial W}{\partial I_{2}}\end{equation}$$

(36)

are the nonlinear material parameters, with |$I_{1}$|⁠, |$I_{2}$| given by (34).

3.2 Shear oscillations of a cuboid of stochastic neo-Hookean material

We now specialize to the case of a cuboid of stochastic neo-Hookean material, with |$\mu _{1}=\mu>0$| and |$\mu _{2}=0$| in (16), where the non-zero components of the Cauchy stress tensor given by (35) are as follows:

$$\begin{equation}\begin{split} T_{xx}&=T_{yy}=-p+\frac{\mu}{\alpha},\\ T_{zz}&=-p+\mu\left(u_{X}^2+u_{Y}^2+\alpha^2\right),\\ T_{xz}&=\frac{\mu}{\sqrt{\alpha}} u_{X},\\ T_{yz}&=\frac{\mu}{\sqrt{\alpha}} u_{Y}.\end{split}\end{equation}$$

(37)

Then by the equation of motion (1),

$$\begin{equation}\begin{split} &\frac{\partial p}{\partial x}=0,\\ &\frac{\partial p}{\partial y}=0,\\ &\frac{\partial p}{\partial z}=-\rho\ddot{u}+\mu\left(u_{XX}+u_{YY}\right),\end{split}\end{equation}$$

(38)

where |$u_{XX}$| and |$u_{YY}$| represent the second derivatives of |$u$| with respect to |$X$| and |$Y$|⁠, respectively. Hence, |$p$| is independent of |$x$| and |$y$|⁠.

We consider the undeformed cuboid to be long in the |$Z$|-direction and impose an initial displacement |$u_{0}(X,Y)=u(X,Y,0)$| and velocity |$\dot{u}_{0}(X,Y)=\dot{u}(X,Y,0)$|⁠. For the boundary condition, we distinguish the following two cases: (i) If we impose null normal Cauchy stresses, |$T_{xx}=T_{yy}=0$|⁠, on the faces perpendicular to the |$X$|- and |$Y$|-directions, at all time, then |$p=\mu /\alpha$| is constant and |$T_{zz}=\mu \left (u_{X}^2+u_{Y}^2+\alpha ^2-1/\alpha \right )$|⁠.

(ii) If |$T_{xx}=T_{yy}\neq 0$|⁠, as |$T_{zz}$| cannot be made point-wise zero, we denote the normal force acting on the cross-sections of area |$A$| in the |$z$|-direction at time |$t$| by

$$\begin{equation}N_{z}(t)=\int_{A} T_{zz}\,\mathrm{d}A,\end{equation}$$

(39)

and consider this force to be zero, i.e., |$N_{z}(t)=0$| at all time. Then |$p$| is independent of |$z$|⁠, and, by (38), it is also independent of |$x$| and |$y$|⁠; hence, |$p=p(t)$|⁠.

In both the above cases, (i) and (ii), respectively, by (38),

$$\begin{equation}\ddot{u}=\frac{\mu}{\rho}\left(u_{XX}+u_{YY}\right).\end{equation}$$

(40)

It remains to solve, by standard procedures, the linear wave equation (40), describing the propagation of waves, subject to the given initial and boundary conditions. To solve this equation, we let the shear stresses |$T_{xz}$| and |$T_{yz}$|⁠, defined by (37), vanish at the sides, i.e.,

$$\begin{equation}\begin{split} T_{xz}(0,Y,Z,t)=T_{xz}(1,Y,Z,t)=0\qquad&\iff\qquad u_{X}(0,Y,t)=u_{X}(1,Y,t)=0,\\ T_{yz}(X,0,Z,t)=T_{yz}(X,1,Z,t)=0\qquad&\iff\qquad u_{Y}(X,0,t)=u_{Y}(X,1,t)=0.\end{split}\end{equation}$$

(41)

In this case, the general solution takes the form

$$\begin{equation}u(X,Y,t)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\left[A_{mn}\cos\left(\omega_{mn}t\right)+B_{mn}\sin\left(\omega_{mn}t\right)\right]\cos\left(\pi mX\right)\cos\left(\pi nY\right),\end{equation}$$

(42)

where

$$\begin{equation}\omega_{mn}=\pi\sqrt{\left(m^2+n^2\right)\frac{\mu}{\rho}},\end{equation}$$

(43)

and

$$\begin{eqnarray}&&A_{mn}=4\int_{0}^{1}\left[\int_{0}^{1}u_{0}(X,Y)\cos\left(\pi mX\right)\,\mathrm{d}X\right]\cos\left(\pi nY\right)\,\mathrm{d}Y,\end{eqnarray}$$

(44)

$$\begin{eqnarray}&&B_{mn}=\frac{4}{\omega_{mn}}\int_{0}^{1}\left[\int_{0}^{1}\dot{u}_{0}(X,Y)\cos\left(\pi mX\right)\,\mathrm{d}X\right]\cos\left(\pi nY\right)\,\mathrm{d}Y.\end{eqnarray}$$

(45)

These oscillations under the generalized shear motion (31) cannot be completely ‘free’, due to the non-zero tractions corresponding to the cases (i) and (ii), respectively. Note that the condition (32) is not satisfied.

Fig. 3.

Stochastic displacement |$u(X,Y,t)$| of the edges |$(X,Y,Z)\in \{(0,0,Z), (1,1,Z)\}$| of the cuboid in dynamic generalized shear, when |$m=n=1$|⁠, |$A_{11}=1$|⁠, |$B_{11}=0$|⁠, |$\rho =1$| and |$\mu$| is drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠. The top figure illustrates the displacement over time of two cuboids, with randomly chosen values of |$\mu$|⁠, derived from the specified Gamma distribution. The middle figure illustrates a probability histogram at each time instant. Specifically, the integral of the probabilities over the displacements at any given time instant is equal to |$1$|⁠. The histogram comprises of |$1000$| stochastic simulations, and the colour bar defines the probability of finding a given displacement at a given time. 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$|⁠. The bottom two figures illustrate specific histogram distributions at two given times (noted above each figure). These are the distributions that would be seen if the middle figure was cut along the green and magenta arrows, respectively.

As |$\mu$| is a random variable, it follows that the speed of wave propagation, |$\sqrt{\mu /\rho }$|⁠, is stochastic. Hence, both the period and the amplitude of the oscillations are stochastic. As an example, we consider the initial data |$u_{0}(X,Y)=\cos (\pi X)\cos (\pi Y)$| and |$\dot{u}_{0}(X,Y)=0$| leading to |$A_{11}=1$| and |$B_{11}=0$|⁠. In Fig. 3, we illustrate the stochastic dynamic displacement on the edges |$(X,Y,Z)\in \{(0,0,Z), (1,1,,Z)\}$| when |$m=n=1$|⁠, |$A_{11}=1$|⁠, |$B_{11}=0$|⁠, |$\rho =1$| and |$\mu$| is drawn from the Gamma distribution with hyperparameters |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠, as represented in Fig. 1. The top plot of Fig. 3 represents two single simulations, with two different values of |$\mu$| drawn from the distribution, illustrating the variety of outcomes that can be obtained. The middle plot of Fig. 3 then represents histograms of the ensemble data. Namely, since not all material parameters are equally likely, not all outcomes are equally likely. Specifically, the values of |$u(0,0,t)$| are most likely going to be near the mean value (dashed line) with the probability of observing alternative values of |$u$| decreasing as we tend away from the mean. We note from Fig. 3 that, as we might expect, extremal probabilities always occur at the extreme displacement of the oscillations, i.e., when the cuboid is slowest. However, in between these probability maxima the variance grows over time. Thus, although the displacements are initially close (seen explicitly at the top of Fig. 3 and by the tight distribution around the mean at the bottom left of Fig. 3), eventually, the phase difference dominates causing the displacements to diverge (top of Figure 3) and an increase in the variance of the distribution (bottom right of Fig. 3).

4. Quasi-equilibrated radial–axial motion of a stochastic hyperelastic cylindrical tube

In this section, we analyse the stability and finite amplitude oscillations of a stochastic hyperelastic cylindrical tube subject to the combined radial and axial quasi-equilibrated dynamic deformation.

4.1 Dynamic radial–axial deformation of a cylindrical tube

For a circular cylindrical tube, the combined radial and axial motion is described by (see Fig. 4)

$$\begin{equation}r^2=a^2+\frac{R^2-A^2}{\alpha},\qquad \theta=\varTheta,\qquad z=\alpha Z,\end{equation}$$

(46)

where |$(R,\varTheta ,Z)$| and |$(r,\theta ,z)$| are the cylindrical polar coordinates in the reference and current configuration, respectively, such that |$A\leq R\leq B$|⁠, |$A$| and |$B$| are the inner and outer radii in the undeformed state, respectively, |$a=a(t)$| and |$b=b(t)=\sqrt{a^2+\left (B^2-A^2\right )/\alpha }$| are the inner and outer radii at time |$t,$| respectively, and |$\alpha>0$| is a given constant (when |$\alpha <0$|⁠, the tube is everted, so that the inner surface becomes the outer surface). When |$\alpha =1$|⁠, the time-dependent deformation (46) simplifies to that studied also in Knowles (1960, 1962) and Beatty (2007). The case when |$\alpha$| is time dependent was considered in Shahinpoor (1973).

Fig. 4.

Schematic of inflation of a cylindrical tube, showing the reference state, with inner radius |$A$| and outer radius |$B$| (left), and the deformed state, with inner radius |$a$| and outer radius |$b$| (right), respectively.

The radial–axial motion (46) of the cylindrical tube is fully determined by the inner radius |$a$| at time |$t$|⁠, which in turn is obtained from the initial conditions. Thus, the acceleration |$\ddot{r}$| can be computed in terms of the acceleration |$\ddot{a}$| on the inner surface. By the governing equations (46), the condition (4) is valid for |$\textbf{x}=(r,\theta ,z)^{T}$|⁠, since

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

(47)

and the acceleration potential, |$\xi$|⁠, satisfies (3). Hence, this is a quasi-equilibrated motion, such that

$$\begin{equation}-\frac{\partial\xi}{\partial r}=\ddot{r}=\frac{\dot{a}^2}{r}+\frac{a\ddot{a}}{r}-\frac{a^2\dot{a}^2}{r^3},\end{equation}$$

(48)

and, by integrating (48), the acceleration potential, |$\xi$|⁠, is given by (Truesdell & Noll, 2004, p. 215)

$$\begin{equation}-\xi=\dot{a}^2\log r+a\ddot{a}\log r+\frac{a^2\dot{a}^2}{2r^2}=\dot{r}^2\log r+r\ddot{r}\log r+\frac{1}{2}\dot{r}^2.\end{equation}$$

(49)

The deformation gradient of (46), with respect to the polar coordinates |$(R,\varTheta ,Z)$|⁠, is equal to

$$\begin{equation}\boldsymbol{F}=\mathrm{diag}\left(\frac{R}{\alpha r}, \frac{r}{R}, \alpha\right),\end{equation}$$

(50)

the Cauchy–Green deformation tensor is

$$\begin{equation}\boldsymbol{B}=\boldsymbol{F}^2=\mathrm{diag}\left(\frac{R^2}{\alpha^2 r^2}, \frac{r^2}{R^2}, \alpha^2\right),\end{equation}$$

(51)

and the principal invariants take the form

$$\begin{equation}\begin{split} I_{1}=&\mathrm{tr}\ (\boldsymbol{B})=\frac{R^2}{\alpha^2 r^2}+\frac{r^2}{R^2}+\alpha^2,\\ I_{2}=&\frac{1}{2}\left[\left(\mathrm{tr}\,\boldsymbol{B}\right)^{2}-\mathrm{tr}\left(\boldsymbol{B}^{2}\right)\right]=\frac{\alpha^2r^2}{R^2}+\frac{R^2}{r^2}+\frac{1}{\alpha^2},\\ I_{3}=&\det\boldsymbol{B}=1.\end{split}\end{equation}$$

(52)

Thus, the principal components of the equilibrium Cauchy stress tensor at time |$t$| are

$$\begin{equation}\begin{split} T^{(0)}_{rr}&=-p^{(0)}+\beta_{1}\frac{R^2}{\alpha^2r^2}+\beta_{-1}\frac{\alpha^2r^2}{R^2},\\ T^{(0)}_{\theta\theta}&=T^{(0)}_{rr}+\left(\beta_{1}-\beta_{-1}\alpha^2\right)\left(\frac{r^2}{R^2}-\frac{R^2}{\alpha^2r^2}\right),\\ T^{(0)}_{zz}&=T^{(0)}_{rr}+\left(\beta_{1}-\beta_{-1}\frac{r^2}{R^2}\right)\left(\alpha^2-\frac{R^2}{\alpha^2r^2}\right),\end{split}\end{equation}$$

(53)

where |$p^{(0)}$| is the Lagrangian multiplier for the incompressibility constraint (⁠|$I_{3}=1$|⁠), and

$$\begin{equation}\beta_{1}=2\frac{\partial W}{\partial I_{1}},\qquad \beta_{-1}=-2\frac{\partial W}{\partial I_{2}}\end{equation}$$

(54)

are the nonlinear material parameters, with |$I_{1}$| and |$I_{2}$| given by (52).

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

$$\begin{equation}\frac{\partial T^{(0)}_{rr}}{\partial r}=\frac{T^{(0)}_{\theta\theta}-T^{(0)}_{rr}}{r}.\end{equation}$$

(55)

Hence, by (53) and (55), the radial Cauchy stress for the equilibrium state at time |$t$| is equal to

$$\begin{equation}T^{(0)}_{rr}(r,t)=\int\left(\beta_{1}-\beta_{-1}\alpha^2\right)\left(\frac{r^2}{R^2}-\frac{R^2}{\alpha^2r^2}\right)\frac{\mathrm{d}r}{r}+\psi(t),\end{equation}$$

(56)

where |$\psi =\psi (t)$| is an arbitrary function of time. Substitution of (49) and (56) into (5) then gives the principal Cauchy stress components at time |$t$| as follows:

$$\begin{equation}\begin{split} T_{rr}(r,t)&=\rho\left(a\ddot{a}\log r+\dot{a}^2\log r+\frac{a^2\dot{a}^2}{2r^2}\right)+\int\left(\beta_{1}-\beta_{-1}\alpha^2\right)\left(\frac{r^2}{R^2}-\frac{R^2}{\alpha^2r^2}\right)\frac{\mathrm{d}r}{r}+\psi(t),\\ T_{\theta\theta}(r,t)&=T_{rr}(r,t)+\left(\beta_{1}-\beta_{-1}\alpha^2\right)\left(\frac{r^2}{R^2}-\frac{R^2}{\alpha^2r^2}\right),\\ T_{zz}(r,t)&=T_{rr}(r,t)+\left(\beta_{1}-\beta_{-1}\frac{r^2}{R^2}\right)\left(\alpha^2-\frac{R^2}{\alpha^2r^2}\right).\end{split}\end{equation}$$

(57)

In (57), the function |$\beta _{1}-\alpha ^2\beta _{-1}$| can be interpreted as the following nonlinear shear modulus (Mihai & Goriely, 2017)

$$\begin{equation}\widetilde{\mu}=\beta_{1}-\beta_{-1}\alpha^2,\end{equation}$$

(58)

corresponding to the combined deformation of simple shear superposed on axial stretch, described by (11), with shear parameter |$k=\sqrt{\alpha ^2R^2/r^2+\alpha ^4r^2/R^2-\alpha ^6-1}$| and stretch parameter |$\alpha$|⁠. As shown in Mihai & Goriely (2017), this modulus is positive if the BE inequalities (9) hold. In this case, the integrand is negative for |$0<r^2/R^2<1/\alpha$| and positive for |$r^2/R^2>1/\alpha$|⁠. Using the first equation in (46), it is straightforward to show that |$0<r^2/R^2<1/\alpha$| (respectively, |$r^2/R^2>1/\alpha$|⁠) is equivalent to |$0<a^2/A^2<1/\alpha$| (respectively, |$a^2/A^2>1/\alpha$|⁠). When |$\alpha =1$|⁠, the modulus defined by (58) coincides with the generalized shear modulus defined in Truesdell & Noll (2004, p. 174), and also in Beatty (2007).

In the limiting case when |$\alpha \to 1$| and |$k\to 0$|⁠, the nonlinear shear modulus defined by (58) converges to the classical shear modulus from the infinitesimal theory (Mihai & Goriely, 2017),

$$\begin{equation}\mu=\lim_{\alpha\to1}\lim_{k\to0}\widetilde{\mu}.\end{equation}$$

(59)

In this case, as |$R^2/r^2\to 1$|⁠, the three stress components defined by (57) are equal.

Next for the cylindrical tube deforming by (46), we set the inner and outer radial pressures acting on the curvilinear surfaces |$r=a(t)$| and |$r=b(t)$| at time |$t$| (measured per unit area in the present configuration), as |$T_{1}(t)$| and |$T_{2}(t),$| respectively (Truesdell & Noll, 2004, pp. 214–217). Evaluating |$T_{1}(t)=-T_{rr}(a,t)$| and |$T_{2}(t)=-T_{rr}(b,t)$|⁠, using (57), with |$r=a$| and |$r=b,$| respectively, then subtracting the results, then gives

$$\begin{equation}\begin{split} T_{1}(t)-T_{2}(t)&=\frac{\rho}{2}\left[\left(a\ddot{a}+\dot{a}^2\right)\log\frac{b^2}{a^2}+\dot{a}^2\left(\frac{a^2}{b^2}-1\right)\right] +\int_{a}^{b}\widetilde{\mu}\left(\frac{r^2}{R^2}-\frac{R^2}{\alpha^2r^2}\right)\frac{\mathrm{d}r}{r}\\ &=\frac{\rho A^2}{2}\left[\left(\frac{a}{A}\frac{\ddot{a}}A+\frac{\dot{a}^2}{A^2}\right)\log\frac{b^2}{a^2}+\frac{\dot{a}^2}{A^2}\left(\frac{a^2}{b^2}-1\right)\right] +\int_{a}^{b}\widetilde{\mu}\left(\frac{r^2}{R^2}-\frac{R^2}{\alpha^2r^2}\right)\frac{\mathrm{d}r}{r}.\end{split}\end{equation}$$

(60)

Setting the notation

$$\begin{equation}u=\frac{r^2}{R^2}=\frac{r^2}{\alpha\left(r^2-a^2\right)+A^2},\qquad x=\frac{a}{A},\qquad \gamma=\frac{B^2}{A^2}-1,\end{equation}$$

(61)

we can rewrite

$$\begin{equation*}\begin{split} \left(\frac{a}{A}\frac{\ddot{a}}A+\frac{\dot{a}^2}{A^2}\right)\log\frac{b^2}{a^2}+\frac{\dot{a}^2}{A^2}\left(\frac{a^2}{b^2}-1\right) &=\left(\ddot{x}x+\dot{x}^2\right)\log\left(1+\frac{\gamma}{\alpha x^2}\right)-\dot{x}^2\frac{\frac{\gamma}{\alpha x^2}}{1+\frac{\gamma}{\alpha x^2}}\\ &=\frac{1}{2x}\frac{d}{dx}\left[\dot{x}^2x^2\log\left(1+\frac{\gamma}{\alpha x^2}\right)\right]\end{split}\end{equation*}$$

and

$$\begin{equation*}\begin{split} \int_{a}^{b}\widetilde{\mu}\left(\frac{r^2}{R^2}-\frac{R^2}{\alpha^2r^2}\right)\frac{\mathrm{d}r}{r} &=\int_{a}^{b}\widetilde{\mu}\left[\frac{r^2}{\alpha\left(r^2-a^2\right)+A^2}-\frac{\alpha\left(r^2-a^2\right)+A^2}{\alpha^2r^2}\right]\frac{\mathrm{d}r}{r}\\ &=\frac{1}{2}\int_{\frac{x^2+\frac{\gamma}{\alpha}}{1+\gamma}}^{x^2}\widetilde{\mu}\frac{1+\alpha u}{\alpha^2u^2}\,\mathrm{d}u.\end{split}\end{equation*}$$

Then we can express the equation (60) equivalently as follows:

$$\begin{equation}2x\frac{T_{1}(t)-T_{2}(t)}{\rho A^2}=\frac{1}{2}\frac{d}{dx}\left[\dot{x}^2x^2\log\left(1+\frac{\gamma}{\alpha x^2}\right)\right] +\frac{x}{\rho A^2}\int_{\frac{x^2+\frac{\gamma}{\alpha}}{1+\gamma}}^{x^2}\widetilde{\mu}\frac{1+\alpha u}{\alpha^2u^2}\,\mathrm{d}u.\end{equation}$$

(62)

Note that, when the BE inequalities (9) hold, |$\widetilde{\mu }>0$|⁠, and the integral in (60), or equivalently in (62), is negative if |$0<u<1/\alpha$| (i.e., if |$0<x<1/\sqrt{\alpha }$|⁠) and positive if |$u>1/\alpha$| (i.e., if |$x>1/\sqrt{\alpha }$|⁠).

In the static case, where |$\dot{a}=0$| and |$\ddot{a}=0$|⁠, (60) becomes

$$\begin{equation}T_{1}(t)-T_{2}(t)=\int_{a}^{b}\widetilde{\mu}\left(\frac{r^2}{R^2}-\frac{R^2}{\alpha^2r^2}\right)\frac{\mathrm{d}r}{r},\end{equation}$$

(63)

and (62) reduces to

$$\begin{equation}2\frac{T_{1}(t)-T_{2}(t)}{\rho A^2}=\frac{1}{\rho A^2}\int_{\frac{x^2+\frac{\gamma}{\alpha}}{1+\gamma}}^{x^2}\widetilde{\mu}\frac{1+\alpha u}{\alpha^2u^2}\,\mathrm{d}u.\end{equation}$$

(64)

For the cylindrical tube in finite dynamic deformation, we set

$$\begin{equation}G(x,\gamma)=\frac{1}{\rho A^2}\int_{1/\sqrt{\alpha}}^{x}\left(\zeta\int_{\frac{\zeta^2+\frac{\gamma}{\alpha}}{1+\gamma}}^{\zeta^2}\widetilde{\mu}\frac{1+\alpha u}{\alpha^2u^2}\,\mathrm{d}u\right)\,\mathrm{d}\zeta,\end{equation}$$

(65)

and find that |$G(x,\gamma )$| is monotonically decreasing when |$0<x<1/\sqrt{\alpha }$| and increasing when |$x>1/\sqrt{\alpha }$|⁠. This function will be useful in establishing whether the radial motion is oscillatory or not.

We also set the pressure impulse (suddenly applied pressure difference)

$$\begin{equation}2\alpha\frac{T_{1}(t)-T_{2}(t)}{\rho A^2}=\left\{ \begin{array}{cc} 0 & \mbox{if}\ t\leq0,\\ p_{0} & \mbox{if}\ t>0,\end{array} \right.\end{equation}$$

(66)

where |$p_{0}$| is constant in time. Then integrating (62) once gives

$$\begin{equation}\frac{1}{2}\dot{x}^2x^2\log\left(1+\frac{\gamma}{\alpha x^2}\right)+G(x,\gamma)=\frac{p_{0}}{2\alpha}\left(x^2-\frac{1}{\alpha}\right)+C,\end{equation}$$

(67)

with |$G(x,\gamma )$| defined by (65) and

$$\begin{equation}C=\frac{1}{2}\dot{x}_{0}^2x_{0}^2\log\left(1+\frac{\gamma}{\alpha x_{0}^2}\right)+G(x_{0},\gamma)-\frac{p_{0}}{2\alpha}\left(x_{0}^2-\frac{1}{\alpha}\right),\end{equation}$$

(68)

where |$x(0)=x_{0}$| and |$\dot{x}(0)=\dot{x}_{0}$| are the initial conditions. By (67),

$$\begin{equation}\dot{x}=\pm\sqrt{\frac{\frac{p_{0}}{\alpha}\left(x^2-\frac{1}{\alpha}\right)+2C-2G(x,\gamma)}{x^2\log\left(1+\frac{\gamma}{\alpha x^2}\right)}}.\end{equation}$$

(69)

Physically, this system is analogous to the motion of a point mass with energy

$$\begin{equation}E=\frac{1}{2} m(x)\dot x^{2}+V(x).\end{equation}$$

(70)

The energy is |$E=C$|⁠, the potential is given by |$V(x)=G(x,\gamma )-\frac{p_{0}}{2\alpha }\left (x^2-\frac{1}{\alpha }\right )$| and the position-dependent mass is |$m(x)={x^2\log \left (1+\frac{\gamma }{\alpha x^2}\right )}$|⁠. Due to the constraints on the function |$G$|⁠, this system has simple dynamics. Depending on the constant |$\mu$|⁠, the system may have a static state or periodic motion. Indeed, the radial motion is periodic if and only if the following equation,

$$\begin{equation}G(x,\gamma)=\frac{p_{0}}{2\alpha}\left(x^2-\frac{1}{\alpha}\right)+C,\end{equation}$$

(71)

has exactly two distinct solutions, representing the amplitudes of the oscillation, |$x=x_{1}$| and |$x=x_{2}$|⁠, such that |$0<x_{1}<x_{2}<\infty$|⁠. Then by (61), the minimum and maximum radii of the inner surface in the oscillation are equal to |$x_{1}A$| and |$x_{2}A,$| respectively, and by (69), the period of oscillation is equal to

$$\begin{equation}T=2\left|\int_{x_{1}}^{x_{2}}\frac{\mathrm{d}x}{\dot{x}}\right|=2\left|\int_{x_{1}}^{x_{2}}\sqrt{\frac{x^2\log\left(1+\frac{\gamma}{\alpha x^2}\right)}{\frac{p_{0}}{\alpha}\left(x^2-\frac{1}{\alpha}\right)+2C-2G(x,\gamma)}}\,\mathrm{d}x\right|.\end{equation}$$

(72)

Note that both the amplitudes and period of the oscillation are random variables described in terms of probability distributions.

4.2 Radial oscillations of a cylindrical tube of stochastic Mooney–Rivlin material

For cylindrical tubes of stochastic Mooney–Rivlin material defined by (16), with |$\mu =\mu _{1}+\mu _{2}>0$|⁠, evaluating the integral in (65) gives (see Appendix  A for detailed calculations)

$$\begin{equation}G(x,\gamma)=\frac{\widetilde{\mu}}{2\alpha\rho A^2}\left(x^2-\frac{1}{\alpha}\right)\log\frac{1+\gamma}{1+\frac{\gamma}{\alpha x^2}},\end{equation}$$

(73)

where |$\widetilde{\mu }=\mu _{1}+\mu _{2}\alpha ^2$|⁠. In this case, assuming that the nonlinear shear modulus |$\mu$| has a uniform lower bound, i.e.,

$$\begin{equation}\mu>\eta,\end{equation}$$

(74)

for some constant |$\eta>0$|⁠, it follows that

$$\begin{equation}\lim_{x\to0}G(x,\gamma)= \lim_{x\to\infty}G(x,\gamma)=\infty.\end{equation}$$

(75)

(i) If |$p_{0}=0$| and |$C>0$|⁠, then equation (71) has exactly two solutions, |$x=x_{1}$| and |$x=x_{2}$|⁠, satisfying |$0<x_{1}<1/\sqrt{\alpha }<x_{2}<\infty$|⁠, for any positive constant |$C$|⁠. It should be noted that, by (57), if |$T_{rr}(r,t)=0$| at |$r=a$| and |$r=b$|⁠, so that |$T_{1}(t)=T_{2}(t)=0$|⁠, then |$T_{\theta \theta }(r,t)\neq 0$| and |$T_{zz}(r,t)\neq 0$| at |$r=a$| and |$r=b$|⁠, unless |$\alpha \to 1$| and |$r^2/R^2\to 1$|⁠. Thus, in general, these oscillations cannot be ‘free’ (Shahinpoor, 1973).

In Fig. 5, for example, we represent the stochastic function |$G(x,\gamma )$|⁠, defined by (73), intersecting the line |$C=10$|⁠, to solve equation (71) when |$p_{0}=0$|⁠, and the associated velocity, given by (69), assuming that |$\alpha =1$|⁠, |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$|⁠, and |$\mu$| follows the Gamma distribution with hyperparameters |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠.

Fig. 5.

The function |$G(x,\gamma )$|⁠, defined by (73), intersecting the (dashed red) line |$C=10$| when |$p_{0}=0$| (left), and the associated velocity, given by (69) (right), for a cylindrical tube of stochastic Mooney–Rivlin material when |$\alpha =1$|⁠, |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$|⁠, and |$\widetilde{\mu }=\mu =\mu _1+\mu _{2}$| is drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠. The dashed black lines correspond to the expected values based only on the mean value, |$\underline{\mu }=\rho _{1}\rho _{2}=4.05$|⁠, of |$\mu$|⁠. Each distribution was calculated from the average of |$1000$| stochastic simulations.

For a thin-walled tube (Knowles, 1960; Shahinpoor & Nowinski, 1971), where |$\alpha =1$| and |$\gamma \to 0$|⁠, equation (67) takes the form

$$\begin{equation}\dot{x}^2+\frac{\mu}{\rho A^2}\left(x^2+\frac{1}{x^2}\right)=\dot{x}_{0}^2+\frac{\mu}{\rho A^2}\left(x_{0}^2+\frac{1}{x_{0}^2}\right),\end{equation}$$

(76)

and has the explicit solution (Shahinpoor & Nowinski, 1971)

$$\begin{equation}x=\sqrt{\left[x_{0}\cos\left(\frac{t}{A}\sqrt{\frac{\mu}{\rho}}\right)+\dot{x}_{0}A\sqrt{\frac{\rho}{\mu}}\sin\left(\frac{t}{A}\sqrt{\frac{\mu}{\rho}}\right)\right]^2+\frac{1}{x_{0}^2}\sin^2\left(\frac{t}{A}\sqrt{\frac{\mu}{\rho}}\right)}.\end{equation}$$

(77)

In this case, assuming that the shear modulus, |$\mu$|⁠, has a uniform lower bound, equation (71) becomes (Knowles, 1960)

$$\begin{equation}x^2+\frac{1}{x^2}=\frac{\rho A^2}{\mu}\dot{x}_{0}^2+x_{0}^2+\frac{1}{x_{0}^2}.\end{equation}$$

(78)

This equation can be solved directly to find the amplitudes

$$\begin{equation}x_{1,2}=\sqrt{\frac{\frac{\rho A^2}{\mu}\dot{x}_{0}^2+x_{0}^2+\frac{1}{x_{0}^2}\pm\sqrt{\left(\frac{\rho A^2}{\mu}\dot{x}_{0}^2+x_{0}^2+\frac{1}{x_{0}^2}\right)^2-4}}{2}}.\end{equation}$$

(79)

Noting that |$x_{2}=1/x_{1}$|⁠, the period of the oscillations can be calculated as

$$\begin{equation}T=2\sqrt{\frac{\rho A^2}{\mu}}\left|\int_{x_{1}}^{1/x_{1}}\frac{\mathrm{d}x}{\sqrt{\frac{\rho A^2}{\mu}\dot{x}_{0}^2+x_{0}^2+\frac{1}{x_{0}^2}-x^2-\frac{1}{x^2}}}\right|=\pi A\sqrt{\frac{\rho}{\mu}}.\end{equation}$$

(80)

In Fig. 6, we illustrate the stochastic solution given by (77), with the initial conditions |$x_{0}=1$| and |$\dot{x}_{0}=4.5$|⁠, assuming that |$\rho =1$|⁠, |$A=1$| and |$\mu$| satisfies the Gamma distribution with hyperparameters |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠.

Fig. 6.

Stochastic solution given by (77), with the initial conditions |$x_{0}=1$| and |$\dot{x}_{0}=4.5$|⁠, for a thin-walled tube, where |$\rho =1$|⁠, |$A=1$| and |$\mu$| is drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠. 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$|⁠. The distribution was calculated from the average of |$1000$| stochastic simulations.

(ii) When |$p_{0}\neq 0$| and |$C\geq 0$|⁠, substitution of (73) in (71) gives

$$\begin{equation}p_{0}=\frac{\widetilde{\mu}}{\rho A^2}\log\frac{1+\gamma}{1+\frac{\gamma}{\alpha x^2}}-\frac{2\alpha C}{x^2-\frac{1}{\alpha}}.\end{equation}$$

(81)

As the right-hand side of the above equation is a monotonically increasing function of |$x$|⁠, there exists a unique positive |$x$| satisfying (81) if and only if the following condition holds:

$$\begin{equation*}\lim_{x\to0}\left(\frac{\widetilde{\mu}}{\rho A^2}\log\frac{1+\gamma}{1+\frac{\gamma}{\alpha x^2}}-\frac{2\alpha C}{x^2-\frac{1}{\alpha}}\right)<p_{0}<\lim_{x\to\infty}\left(\frac{\widetilde{\mu}}{\rho A^2}\log\frac{1+\gamma}{1+\frac{\gamma}{\alpha x^2}}-\frac{2\alpha C}{x^2-\frac{1}{\alpha}}\right),\end{equation*}$$

that is,

$$\begin{equation}-\infty<p_{0}<\frac{\widetilde{\mu}}{\rho A^2}\log\left(1+\gamma\right).\end{equation}$$

(82)

Then by (61), (66) and (82), the necessary and sufficient condition that oscillatory motions occur is that the nonlinear shear modulus, |$\widetilde{\mu }$|⁠, is uniformly bounded from below as follows:

$$\begin{equation}\widetilde{\mu}>\frac{p_{0}\rho A^2}{\log\left(1+\gamma\right)}=\alpha\frac{T_{1}(t)-T_{2}(t)}{\log B-\log A}.\end{equation}$$

(83)

By (58),

$$\begin{equation*}\widetilde{\mu}=\mu_{1}+\mu_{2}\alpha^2=\mu_{1}+\left(\mu-\mu_{1}\right)\alpha^2=\mu\alpha^2+\mu_{1}\left(1-\alpha^2\right).\end{equation*}$$

Hence, (83) is equivalent to

$$\begin{equation}\mu>\frac{p_{0}\rho A^2}{\alpha^2\log\left(1+\gamma\right)}+\mu_{1}\frac{1-\alpha^2}{\alpha^2}.\end{equation}$$

(84)

Then, the probability distribution of oscillatory motions occurring is

$$\begin{equation}P_{1}(p_0)=1-\int_{0}^{\frac{p_{0}\rho A^2}{\alpha^2\log\left(1+\gamma\right)}+\mu_{1}\frac{1-\alpha^2}{\alpha^2}}g(u;\rho_{1},\rho_{2})\,\mathrm{d}u,\end{equation}$$

(85)

where |$g(u;\rho _{1},\rho _{2})$| is the Gamma probability density function defined by (19), and that of non-oscillatory motions is

$$\begin{equation}P_{2}(p_0)=1-P_{1}(p_0)=\int_{0}^{\frac{p_{0}\rho A^2}{\alpha^2\log\left(1+\gamma\right)}+\mu_{1}\frac{1-\alpha^2}{\alpha^2}}g(u;\rho_{1},\rho_{2})\,\mathrm{d}u.\end{equation}$$

(86)

For example, when |$\alpha =1$|⁠, |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$| and |$\widetilde{\mu }=\mu =\mu _1+\mu _{2}$| satisfies the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠, the probability distributions given by (85) and (86) are shown in Fig. 7 (blue lines for |$P_{1}$| and red lines for |$P_{2}$|⁠). Specifically, |$(0,\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 tube, the critical value |$p_{0}=\underline{\mu }\log 2\approx 2.8072$| 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 of 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 impulse, |$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.

Fig. 7.

Probability distributions of whether oscillatory motions can occur or not for a cylindrical tube of stochastic Mooney–Rivin material, with |$\alpha =1$|⁠, |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$| and the shear modulus, |$\mu$|⁠, following the Gamma distribution with |$\rho _{1}=405$|⁠, |$\rho _{2}=0.01$|⁠. Dark-coloured lines represent analytically derived solutions, given by equations (85) and (86), whereas the lighter versions represent stochastically generated data. The vertical line at the critical value, |$p_{0}=2.8072$|⁠, separates the expected regions based only on the mean value of the shear modulus, |$\underline{\mu }=\rho _{1}\rho _{2}=4.05$|⁠. The probabilities were calculated from the average of 100 stochastic simulations.

In Fig. 8, we illustrate the stochastic function |$G(x,\gamma )$|⁠, defined by (73), intersecting the curve |$p_{0}\left (x^2-1/\alpha \right )/(2\alpha )+C$|⁠, with |$p_{0}=1$| and |$C=7$|⁠, to find the solutions of equation (71), and the associated velocity, given by (69), assuming that |$\alpha =1$|⁠, |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$|⁠, and |$\mu$| satisfies the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠.

Fig. 8.

The function |$G(x,\gamma )$|⁠, defined by (73), intersecting the (dashed red) curve |$p_{0}\left (x^2-1/\alpha \right )/(2\alpha )+C$|⁠, with |$p_{0}=1$| and |$C=7$|⁠, (left), and the associated velocity, given by (69) (right), for a cylindrical tube of stochastic Mooney–Rivlin material when |$\alpha =1$|⁠, |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$| and |$\widetilde{\mu }=\mu =\mu _1+\mu _{2}$| is drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠. The dashed black lines correspond to the expected values based only on the mean value, |$\underline{\mu }=\rho _{1}\rho _{2}=4.05$|⁠, of |$\mu$|⁠. Each distribution was calculated from the average of |$1000$| stochastic simulations.

When |$C=0$|⁠, equation (81) can be solved explicitly to find the amplitude

$$\begin{equation}x_{1}=\sqrt{\frac{\gamma/\alpha}{\left(1+\gamma\right)\exp\left[-\left(p_{0}\rho A^2\right)/(\widetilde{\mu})\right]-1}} =\sqrt{\frac{\left(B^2-A^2\right)/\alpha}{B^2\exp\left[-2\alpha\left(P_{1}-P_{2}\right)/\widetilde{\mu}\right]-A^2}}.\end{equation}$$

(87)

Note that, in the static case, by (64) and (66), at |$x=x_{1}$|⁠, the required pressure takes the form

$$\begin{equation}p_{0}^{(s)}=\frac{\widetilde{\mu}}{\alpha x^2\rho A^2}\frac{\gamma-\frac{\gamma}{\alpha x^2}}{1+\frac{\gamma}{\alpha x^2}}+\frac{\widetilde{\mu}}{\rho A^2}\log\frac{1+\gamma}{1+\frac{\gamma}{\alpha x^2}}.\end{equation}$$

(88)

Thus, the difference between the applied pressure in the static and dynamic case, given by (88) and (81), with |$C=0$|⁠, respectively, is

$$\begin{equation}p_{0}^{(s)}-p_{0}=\frac{\widetilde{\mu}}{\alpha x^2\rho A^2}\frac{\gamma-\frac{\gamma}{\alpha x^2}}{1+\frac{\gamma}{\alpha x^2}}.\end{equation}$$

(89)

Hence, |$p_{0}^{(s)}<p_{0}$| if |$0<x_{1}<\sqrt{\alpha }$|⁠, and |$p_{0}^{(s)}>p_{0}$| if |$x_{1}>\sqrt{\alpha }$|⁠.

If the tube wall is thin (Knowles, 1962; Shahinpoor & Nowinski, 1971), then |$0<\gamma \ll 1$| and |$\alpha =1$|⁠, and (81) becomes

$$\begin{equation}\frac{p_{0}}{\gamma}=\frac{\mu}{\rho A^2}\left(1-\frac{1}{x^2}\right)-\frac{2\alpha C}{x^2-\frac{1}{\alpha}}.\end{equation}$$

(90)

Then the necessary and sufficient condition that oscillatory motions occur is that

$$\begin{equation}-\infty=\lim_{x\to0}\left[\frac{\mu}{\rho A^2}\left(1-\frac{1}{x^2}\right)-\frac{2\alpha C}{x^2-\frac{1}{\alpha}}\right]<\frac{p_{0}}{\gamma}<\lim_{x\to\infty}\left[\frac{\mu}{\rho A^2}\left(1-\frac{1}{x^2}\right)-\frac{2\alpha C}{x^2-\frac{1}{\alpha}}\right]=\frac{\mu}{\rho A^2}.\end{equation}$$

(91)

Thus, for the motion to be oscillatory, the shear modulus must be bounded from below as follows:

$$\begin{equation}\mu>\frac{p_{0}}{\gamma}\rho A^2=\frac{2}{\gamma}\left(T_{1}(t)-T_{2}(t)\right).\end{equation}$$

(92)

Then the probability distribution of oscillatory motions occurring is

$$\begin{equation}P_{1}(p_0/\gamma)=1-\int_{0}^{\frac{p_{0}}{\gamma}\rho A^2}g(u;\rho_{1},\rho_{2})\,\mathrm{d}u,\end{equation}$$

(93)

and that of non-oscillatory motions is

$$\begin{equation}P_{2}(p_0/\gamma)=1-P_{1}(p_0/\gamma)=\int_{0}^{\frac{p_{0}}{\gamma}\rho A^2}g(u;\rho_{1},\rho_{2})\,\mathrm{d}u.\end{equation}$$

(94)

For |$\rho =1$|⁠, |$A=1$| and |$\widetilde{\mu }=\mu =\mu _1+\mu _{2}$| drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠, the probability distributions given by (93) and (94) are shown in Fig. 9 (blue lines for |$P_{1}$| and red lines for |$P_{2}$|⁠). For the deterministic thin-walled tube, the critical value |$p_{0}/\gamma =\underline{\mu }=4.05$| strictly separates the cases of oscillations occurring or not. However, in the stochastic case, the two cases compete.

Fig. 9.

Probability distributions of whether oscillatory motions can occur or not for a thin-walled cylindrical tube of stochastic Mooney–Rivin material, with |$\rho =1$|⁠, |$A=1$| and the shear modulus, |$\mu$|⁠, following the Gamma distribution with |$\rho _{1}=405$|⁠, |$\rho _{2}=0.01$|⁠. Dark-coloured lines represent analytically derived solutions, given by equations (85) and (86), whereas the lighter versions represent stochastically generated data. The vertical line at the critical value, |$p_{0}/\gamma =4.05$|⁠, separates the expected regions based only on the mean value of the shear modulus, |$\underline{\mu }=\rho _{1}\rho _{2}=4.05$|⁠. The probabilities were calculated from the average of 100 stochastic simulations.

If |$C=0$|⁠, then setting |$x_{0}=1$| and |$\dot{x}_{0}=0$|⁠, the equation of motion has the explicit solution (Shahinpoor & Nowinski, 1971)

$$\begin{equation}x=\sqrt{\frac{\frac{\mu}{\rho A^2}-\frac{p_{0}}{2\gamma}}{\frac{\mu}{\rho A^2}-\frac{p_{0}}{\gamma}}-\frac{\frac{p_{0}}{2\gamma}}{\frac{\mu}{\rho A^2}-\frac{p_{0}}{\gamma}}\cos\left(2t\sqrt{\frac{\mu}{\rho A^2}-\frac{p_{0}}{\gamma}}\right)}.\end{equation}$$

(95)

In Fig. 10, we illustrate the stochastic solution given by (95), with |$p_{0}/\gamma =1$|⁠, assuming that |$\rho =1$|⁠, |$A=1$|⁠, and |$\mu$| satisfies the Gamma distribution with hyperparameters |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠.

Fig. 10.

Stochastic solution given by (95), with |$p_{0}/\gamma =1$|⁠, for a thin-walled tube, where |$\rho =1$|⁠, |$A=1$| and |$\mu$| is drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠. 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$|⁠. The distribution was calculated from the average of |$1000$| stochastic simulations.

If the tube wall is infinitely thick (Shahinpoor, 1973), then |$\gamma \to \infty$|⁠, and assuming that the nonlinear shear modulus, |$\widetilde{\mu }$|⁠, has a uniform lower bound, (82) becomes

$$\begin{equation}-\infty=\lim_{x\to0}\left[\frac{\widetilde{\mu}}{\rho A^2}\log\left(\alpha x^2\right)-\frac{2\alpha C}{x^2-\frac{1}{\alpha}}\right]<p_{0}<\lim_{x\to\infty}\left[\frac{\widetilde{\mu}}{\rho A^2}\log\left(\alpha x^2\right)-\frac{2\alpha C}{x^2-\frac{1}{\alpha}}\right]=\infty.\end{equation}$$

(96)

Hence, the motion is always oscillatory for any value of the applied impulse.

5. Quasi-equilibrated radial motion of a stochastic hyperelastic spherical shell

Next we examine the stability and finite amplitude oscillations of a stochastic hyperelastic spherical shell under quasi-equilibrated dynamic radial deformation.

5.1 Dynamic radial deformation of a spherical shell

For a spherical shell, the radial motion is described by (Balakrishnan & Shahinpoor, 1978; Beatty, 2011; Heng & Solecki, 1963; Knowles & Jakub, 1965) (see Fig. 11)

$$\begin{equation}r^3=a^3+R^3-A^3,\qquad \theta=\varTheta,\qquad \phi=\varPhi,\end{equation}$$

(97)

where |$(R,\varTheta ,\varPhi )$| and |$(r,\theta ,\phi )$| are the spherical polar coordinates in the reference and current configuration, respectively, such that |$A\leq R\leq B$|⁠, |$A$| and |$B$| are the inner and outer radii in the undeformed state, and |$a=a(t)$| and |$b=b(t)=\sqrt [3]{a^3+B^3-A^3}$| are the inner and outer radii at time |$t,$| respectively.

Fig. 11.

Schematic of inflation of a spherical shell, showing the reference state, with inner radius |$A$| and outer radius |$B$| (left), and the deformed state, with inner radius |$a$| and outer radius |$b$| (right), respectively.

As for the cylindrical tube, the radial motion (97) of the spherical shell is determined entirely by the inner radius |$a$| at time |$t$|⁠. By the governing equations (97), the condition (4) is valid for |$\textbf{x}=(r,\theta ,\phi )^{T}$|⁠, since

$$\begin{equation}\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],\end{equation}$$

(98)

and the acceleration potential, |$\xi$|⁠, satisfies (3). Hence, this is a quasi-equilibrated motion, such that

$$\begin{equation}-\frac{\partial\xi}{\partial r}=\ddot{r}=\frac{2a\dot{a}^2+a^2\ddot{a}}{r^2}-\frac{2a^4\dot{a}^2}{r^5},\end{equation}$$

(99)

and integrating (99) gives (Truesdell & Noll, 2004, p. 217)

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

(100)

For the deformation (97), the gradient tensor with respect to the polar coordinates |$(R,\varTheta ,\varPhi )$| takes the form

$$\begin{equation}\boldsymbol{F}=\mathrm{diag}\left(\frac{R^2}{r^2}, \frac{r}{R}, \frac{r}{R}\right),\end{equation}$$

(101)

the Cauchy–Green tensor is equal to

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

(102)

and the corresponding principal invariants are

$$\begin{equation}\begin{split} I_{1}=&\mathrm{tr}\ (\boldsymbol{B})=\frac{R^4}{r^4}+2\frac{r^2}{R^2},\\ I_{2}=&\frac{1}{2}\left[\left(\mathrm{tr}\,\boldsymbol{B}\right)^{2}-\mathrm{tr}\left(\boldsymbol{B}^{2}\right)\right]=\frac{r^4}{R^4}+2\frac{R^2}{r^2},\\ I_{3}=&\det\boldsymbol{B}=1.\end{split}\end{equation}$$

(103)

Then the principal components of the equilibrium Cauchy stress at time |$t$| are

$$\begin{equation}\begin{split} T^{(0)}_{rr}&=-p^{(0)}+\beta_{1}\frac{R^4}{r^4}+\beta_{-1}\frac{r^4}{R^4},\\ T^{(0)}_{\theta\theta}&=T^{(0)}_{rr}+\left(\beta_{1}-\beta_{-1}\frac{r^2}{R^2}\right)\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right),\\ T^{(0)}_{\phi\phi}&=T^{(0)}_{\theta\theta},\end{split}\end{equation}$$

(104)

where |$p^{(0)}$| is the Lagrangian multiplier for the incompressibility constraint (⁠|$I_{3}=1$|⁠), and

$$\begin{equation}\beta_{1}=2\frac{\partial W}{\partial I_{1}},\qquad \beta_{-1}=-2\frac{\partial W}{\partial I_{2}},\end{equation}$$

(105)

with |$I_{1}$| and |$I_{2}$| given by (103).

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

$$\begin{equation}\frac{\partial T^{(0)}_{rr}}{\partial r}=2\frac{T^{(0)}_{\theta\theta}-T^{(0)}_{rr}}{r}.\end{equation}$$

(106)

Hence, by (104) and (106), the radial Cauchy stress for the equilibrium state at |$t$| is equal to

$$\begin{equation}T^{(0)}_{rr}(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{\mathrm{d}r}{r}+\psi(t),\end{equation}$$

(107)

where |$\psi =\psi (t)$| is an arbitrary function of time. Substitution of (100) and (107) into (5) gives the following principal Cauchy stresses at time |$t$|⁠:

$$\begin{equation}\begin{split} T_{rr}(r,t)&=-\rho\left(\frac{a^2\ddot{a}+2a\dot{a}^2}{r}-\frac{a^4\dot{a}^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{\mathrm{d}r}{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_{\phi\phi}(r,t)&=T_{\theta\theta}(r,t).\end{split}\end{equation}$$

(108)

In (108), the function |$\beta _{1}-\beta _{-1}\left (r^2/R^2\right )$| can be regarded as the following nonlinear shear modulus (Beatty, 2011; Mihai & Goriely, 2017):

$$\begin{equation}\widetilde{\mu}=\beta_{1}-\beta_{-1}\frac{r^2}{R^2},\end{equation}$$

(109)

corresponding to the combined deformation of infinitesimal shear superposed on finite axial stretch, defined by (11), with the shear parameter satisfying |$k\to 0$| and the stretch parameter |$\alpha =r/R$|⁠. This modulus is positive if the BE inequalities (9) hold (Mihai & Goriely, 2017). In this case, the integrand in (108) is negative for |$0<r^2/R^2<1$| (i.e., when |$0<a^2/A^2<1$|⁠) and positive for |$r^2/R^2>1$| (i.e., when |$a^2/A^2>1$|⁠).

When |$R^{2}/r^{2}\to 1$|⁠, the nonlinear elastic modulus given by (109) converges to the shear modulus from linear elasticity,

$$\begin{equation}\mu=\lim_{R^{2}/r^{2}\to 1}\widetilde{\mu}.\end{equation}$$

(110)

In this case, the stress components given by (108) are equal.

For the spherical shell deforming by (97), we set the inner and outer radial pressures acting on the curvilinear surfaces, |$r=a(t)$| and |$r=b(t)$| at time |$t$|⁠, as |$T_{1}(t)$| and |$T_{2}(t),$| respectively (Truesdell & Noll, 2004, pp. 217–219). Then evaluating |$T_{1}(t)=-T_{rr}(a,t)$| and |$T_{2}(t)=-T_{rr}(b,t)$|⁠, using (108), with |$r=a$| and |$r=b,$| respectively, and subtracting the results, gives

$$\begin{equation}\begin{split} T_{1}(t)-T_{2}(t)&=\rho\left[\left(a^2\ddot{a}+2a\dot{a}^2\right)\left(\frac{1}{a}-\frac{1}{b}\right)-\frac{a^4\dot{a}^2}{2}\left(\frac{1}{a^4}-\frac{1}{b^4}\right)\right]+2\int_{a}^{b}\widetilde{\mu}\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{\mathrm{d}r}{r}\\ &=\rho\left[\left(a\ddot{a}+2\dot{a}^2\right)\left(1-\frac{a}{b}\right)-\frac{\dot{a}^2}{2}\left(1-\frac{a^4}{b^4}\right)\right]+2\int_{a}^{b}\widetilde{\mu}\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{\mathrm{d}r}{r}\\ &=\rho A^2\left[\left(\frac{a}{A}\frac{\ddot{a}}A+2\frac{\dot{a}^2}{A^2}\right)\left(1-\frac{a}{b}\right)-\frac{\dot{a}^2}{2A^2}\left(1-\frac{a^4}{b^4}\right)\right]+2\int_{a}^{b}\widetilde{\mu}\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{\mathrm{d}r}{r}.\end{split}\end{equation}$$

(111)

Setting the notation

$$\begin{equation}u=\frac{r^3}{R^3}=\frac{r^3}{r^3-a^3+A^3},\qquad x=\frac{a}{A},\qquad \gamma=\frac{B^3}{A^3}-1,\end{equation}$$

(112)

we can rewrite

$$\begin{equation*}\begin{split} \left(\frac{a}{A}\frac{\ddot{a}}A+2\frac{\dot{a}^2}{A^2}\right)\left(1-\frac{a}{b}\right)&-\frac{\dot{a}^2}{2A^2}\left(1-\frac{a^4}{b^4}\right)\\ &=\left(\ddot{x}x+2\dot{x}^2\right)\left[1-\left(1+\frac{\gamma}{x^3}\right)^{-1/3}\right]-\frac{\dot{x}^2}{2}\left[1-\left(1+\frac{\gamma}{x^3}\right)^{-4/3}\right]\\ &=\left(\ddot{x}x+\frac{3}{2}\dot{x}^2\right)\left[1-\left(1+\frac{\gamma}{x^3}\right)^{-1/3}\right]-\frac{\dot{x}^2}{2}\frac{\gamma}{x^3}\left(1+\frac{\gamma}{x^3}\right)^{-4/3}\\ &=\frac{1}{2x^2}\frac{d}{dx}\left\{\dot{x}^2x^3\left[1-\left(1+\frac{\gamma}{x^3}\right)^{-1/3}\right]\right\}\end{split}\end{equation*}$$

and

$$\begin{equation*}\begin{split} \int_{a}^{b}\widetilde{\mu}\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{\mathrm{d}r}{r} &=\int_{a}^{b}\widetilde{\mu}\left[{\left(\frac{r^3}{r^3-a^3+A^3}\right)^{2/3}}-\left(\frac{r^3-a^3+A^3}{r^3}\right)^{4/3}\right]\frac{\mathrm{d}r}{r}\\ &=\frac{1}{3}\int_{\frac{x^3+\gamma}{1+\gamma}}^{x^3}\widetilde{\mu}\frac{1+u}{u^{7/3}}\mathrm{d}u.\end{split}\end{equation*}$$

Hence, (111) can be written equivalently as follows:

$$\begin{equation}2x^2\frac{T_{1}(t)-T_{2}(t)}{\rho A^2} =\frac{d}{dx}\left\{\dot{x}^2x^3\left[1-\left(1+\frac{\gamma}{x^3}\right)^{-1/3}\right]\right\} +\frac{4x^2}{3\rho A^2}\int_{\frac{x^3+\gamma}{1+\gamma}}^{x^3}\widetilde{\mu}\frac{1+u}{u^{7/3}}\mathrm{d}u.\end{equation}$$

(113)

Note that, when the BE inequalities (9) hold, |$\widetilde{\mu }>0$|⁠, and the integral in (113) is negative if |$0<x<1$| and positive if |$x>1$|⁠.

In the static case, (111) reduces to

$$\begin{equation}T_{1}(t)-T_{2}(t)=\int_{a}^{b}\widetilde{\mu}\left(\frac{r^2}{R^2}-\frac{R^4}{r^4}\right)\frac{\mathrm{d}r}{r},\end{equation}$$

(114)

and (113) becomes

$$\begin{equation}2\frac{T_{1}(t)-T_{2}(t)}{\rho A^2}=\frac{4}{3\rho A^2}\int_{\frac{x^3+\gamma}{1+\gamma}}^{x^3}\widetilde{\mu}\frac{1+u}{u^{7/3}}\mathrm{d}u.\end{equation}$$

(115)

For the dynamic spherical shell, we set

$$\begin{equation}H(x,\gamma)=\frac{4}{3\rho A^2}\int_{1}^{x}\left(\zeta^2\int_{\frac{\zeta^3+\gamma}{1+\gamma}}^{\zeta^3}\widetilde{\mu}\frac{1+u}{u^{7/3}}\mathrm{d}u\right)\mathrm{d}\zeta,\end{equation}$$

(116)

and obtain that |$H(x,\gamma )$| is monotonically decreasing when |$0<x<1$| and increasing when |$x>1$|⁠.

We also set a pressure impulse that is constant in time,

$$\begin{equation}2\frac{T_{1}(t)-T_{2}(t)}{\rho A^2}=\left\{ \begin{array}{cc} 0 & \mbox{if}\ t\leq0,\\ p_{0} & \mbox{if}\ t>0.\end{array} \right.\end{equation}$$

(117)

Then integrating (113) once gives

$$\begin{equation}\dot{x}^2x^3\left[1-\left(1+\frac{\gamma}{x^3}\right)^{-1/3}\right]+H(x,\gamma)=\frac{p_{0}}{3}\left(x^3-1\right)+C,\end{equation}$$

(118)

with |$H(x,\gamma )$| defined by (116), and

$$\begin{equation}C=\dot{x}_{0}^2x_{0}^3\left[1-\left(1+\frac{\gamma}{x^3}\right)^{-1/3}\right]+H(x_{0},\gamma)-\frac{p_{0}}{3}\left(x_{0}^3-1\right),\end{equation}$$

(119)

where |$x(0)=x_{0}$| and |$\dot{x}(0)=\dot{x}_{0}$| are the initial conditions. From (118), we obtain

$$\begin{equation}\dot{x}=\pm\sqrt{\frac{\frac{p_{0}}{3}\left(x^3-1\right)+C-H(x,\gamma)}{x^3\left[1-\left(1+\frac{\gamma}{x^3}\right)^{-1/3}\right]}}.\end{equation}$$

(120)

The analogy with the motion of a point mass in a potential still holds with appropriate modification. Hence, oscillatory motion of the spherical shell occurs if and only if the following equation,

$$\begin{equation}H(x,\gamma)=\frac{p_{0}}{3}\left(x^3-1\right)+C,\end{equation}$$

(121)

has exactly two distinct solutions, representing the amplitudes of the oscillation, |$x=x_{1}$| and |$x=x_{2}$|⁠, such that |$0<x_{1}<x_{2}<\infty$|⁠. In this case, the minimum and maximum radii of the inner surface in the oscillation are given by |$x_{1}A$| and |$x_{2}A,$| respectively, and the period of oscillation is equal to

$$\begin{equation}T=2\left|\int_{x_{1}}^{x_{2}}\frac{\mathrm{d}x}{\dot{x}}\right|=2\left|\int_{x_{1}}^{x_{2}}\sqrt{\frac{x^3\left[1-\left(1+\frac{\gamma}{x^3}\right)^{-1/3}\right]}{\frac{p_{0}}{3}\left(x^3-1\right)+C-H(x,\gamma)}}\mathrm{d}x\right|.\end{equation}$$

(122)

Note that the amplitude and the period of the oscillation are random variables characterized by probability distributions.

5.2 Radial oscillations of a spherical shell of stochastic neo-Hookean material

For a spherical shell of stochastic neo-Hookean material, with |$\mu _{1}=\mu>0$| and |$\mu _{2}=0$| in (16), evaluating the integral in (116) gives (see Appendix  A for a detailed derivation)

$$\begin{equation}H(x,\gamma)=\frac{\mu}{\rho A^2}\left(x^3-1\right)\left[\frac{2x^3-1}{x^3+x^2+x} -\frac{2\frac{x^3+\gamma}{1+\gamma}-1}{\frac{x^3+\gamma}{1+\gamma}+\left(\frac{x^3+\gamma}{1+\gamma}\right)^{2/3}+\left(\frac{x^3+\gamma}{1+\gamma}\right)^{1/3}}\right].\end{equation}$$

(123)

Assuming that the nonlinear shear modulus, |$\mu$|⁠, is uniformly bounded from below, i.e.,

$$\begin{equation}\mu>\eta,\end{equation}$$

(124)

for some constant |$\eta>0$|⁠, it follows that

$$\begin{equation}\lim_{x\to0}H(x,\gamma)=\lim_{x\to\infty}H(x,\gamma)=\infty.\end{equation}$$

(125)

(i)When |$p_{0}=0$| and |$C>0$|⁠, equation (121) has exactly two solutions, |$x=x_{1}$| and |$x=x_{2}$|⁠, satisfying |$0<x_{1}<1<x_{2}<\infty$|⁠, for any positive constant |$C$|⁠. In this case, it should be noted that, by (108), if |$T_{rr}(r,t)=0$| at |$r=a$| and |$r=b$|⁠, so that |$T_{1}(t)=T_{2}(t)=0$|⁠, then, |$T_{\theta \theta }(r,t)=T_{\phi \phi }(r,t)\neq 0$| at |$r=a$| and |$r=b$|⁠, unless |$r^3/R^3\to 1$|⁠. Thus, the oscillations cannot be considered as ‘free’ in general.

In Fig. 12, we show the stochastic function |$H(x,\gamma )$|⁠, defined by (123), intersecting the line |$C=10$|⁠, to solve equation (121) when |$p_{0}=0$|⁠, and the associated velocity, given by (120), assuming that |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$| and |$\mu$| follows the Gamma distribution with hyperparameters |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠.

Fig. 12.

The function |$H(x,\gamma )$|⁠, defined by (123), intersecting the (dashed red) line |$C=10$|⁠, when |$p_{0}=0$| (left), and the associated velocity, given by (120) (right), for the spherical shell of stochastic neo-Hookean material, where |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$| and |$\mu$| is drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠. The dashed black lines correspond to the expected values based only on the mean value, |$\underline{\mu }=\rho _{1}\rho _{2}=4.05$|⁠, of |$\mu$|⁠. Each distribution was calculated from the average of |$1000$| stochastic simulations.

(ii)When |$p_{0}\neq 0$| and |$C\geq 0$|⁠, substitution of (123) in (121) gives

$$\begin{equation}p_{0}=\frac{3\mu}{\rho A^2}\left[\frac{2x^3-1}{x^3+x^2+x} -\frac{2\frac{x^3+\gamma}{1+\gamma}-1}{\frac{x^3+\gamma}{1+\gamma}+\left(\frac{x^3+\gamma}{1+\gamma}\right)^{2/3}+\left(\frac{x^3+\gamma}{1+\gamma}\right)^{1/3}}\right]-\frac{3C}{x^3-1}.\end{equation}$$

(126)

As the right-hand side of (126) is function of |$x$| that monotonically increases from |$-\infty$| as |$x\to 0$| to |$\infty$| as |$x\to \infty$|⁠, the motion is oscillatory for all values of the given pressure difference.

In the static case, by (115) and (117), the applied pressure takes the form

$$\begin{equation}p_{0}^{(s)}=\frac{\mu}{\rho A^2}\left[\left(\frac{1+\gamma}{x^3+\gamma}\right)^{4/3}+4\left(\frac{1+\gamma}{x^3+\gamma}\right)^{1/3}-\frac{1}{x^4}-\frac{4}{x}\right].\end{equation}$$

(127)

In Fig. 13, we represent the stochastic function |$H(x,\gamma )$|⁠, defined by (123), intersecting the curve |$p_{0}\left (x^3-1\right )/3+C$|⁠, with |$p_{0}=1$| and |$C=3$|⁠, to obtain the solutions of equation (121), and the associated velocity, given by (120), assuming that |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$| and |$\mu$| follows the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠.

Fig. 13.

The function |$H(x,\gamma )$|⁠, defined by (123), intersecting the (dashed red) curve |$p_{0}\left (x^3-1\right )/3+C$|⁠, with |$p_{0}=1$| and |$C=3$| (left), and the associated velocity, given by (120) (right), for the spherical shell of stochastic neo-Hookean material, where |$\rho =1$|⁠, |$A=1$|⁠, |$\gamma =1$| and |$\mu$| is drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠. The dashed black lines correspond to the expected values based only on the mean value, |$\underline{\mu }=\rho _{1}\rho _{2}=4.05$|⁠, of |$\mu$|⁠. Each distribution was calculated from the average of |$1000$| stochastic simulations.

If the spherical shell has an infinitely thick wall (Balakrishnan & Shahinpoor, 1978; Knowles & Jakub, 1965), then |$\gamma \to \infty$|⁠, and the necessary and sufficient condition for the motion to be oscillatory becomes

$$\begin{equation*}\lim_{x\to0}\left[\frac{3\mu}{\rho A^2}\left(\frac{2x^3-1}{x^3+x^2+x}-\frac{1}{3}\right)-\frac{3C}{x^3-1}\right]<p_{0}<\lim_{x\to\infty}\left[\frac{3\mu}{\rho A^2}\left(\frac{2x^3-1}{x^3+x^2+x}-\frac{1}{3}\right)-\frac{3C}{x^3-1}\right],\end{equation*}$$

that is

$$\begin{equation}-\infty<p_{0}<\frac{5\mu}{\rho A^2}.\end{equation}$$

(128)

Thus, for the oscillations to occur, the shear modulus must satisfy (Knowles & Jakub, 1965)

$$\begin{equation}\mu>p_{0}\frac{\rho A^2}{5}=\frac{2}{5}\left(T_{1}(t)-T_{2}(t)\right).\end{equation}$$

(129)

Then the probability distribution of oscillatory motions occurring is

$$\begin{equation}P_{1}(p_0)=1-\int_{0}^{p_{0}\frac{\rho A^2}{5}}g(u;\rho_{1},\rho_{2})\,\mathrm{d}u,\end{equation}$$

(130)

and that of non-oscillatory motions is

$$\begin{equation}P_{2}(p_0)=1-P_{1}(p_0)=\int_{0}^{p_{0}\frac{\rho A^2}{5}}g(u;\rho_{1},\rho_{2})\,\mathrm{d}u.\end{equation}$$

(131)

For |$\rho =1$|⁠, |$A=1$| and |$\widetilde{\mu }=\mu =\mu _1+\mu _{2}$| drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠, the probability distributions given by (130) and (131) are shown in Fig. 14 (blue lines for |$P_{1}$| and red lines for |$P_{2}$|⁠). For the deterministic thin-walled tube, the critical value |$p_{0}=5\underline{\mu }=20.25$| strictly separates the cases of oscillations occurring or not. However, in the stochastic case, there is competition between the two cases.

Fig. 14.

Probability distributions of whether oscillatory motions can occur or not for an infinitely thick-walled spherical shell of stochastic neo-Hookean material, with |$\rho =1$|⁠, |$A=1$| and the shear modulus, |$\mu$|⁠, following the Gamma distribution with |$\rho _{1}=405$|⁠, |$\rho _{2}=0.01$|⁠. Dark-coloured lines represent analytically derived solutions, given by equations (130) and (131), whereas the lighter versions represent stochastically generated data. The vertical line at the critical value, |$p_{0}=20.25$|⁠, separates the expected regions based only on the mean value of the shear modulus, |$\underline{\mu }=\rho _{1}\rho _{2}=4.05$|⁠. The probabilities were calculated from the average of 100 stochastic simulations.

If the spherical shell wall is thin (Beatty, 2011; Verron et al., 1999; Wang, 1965), then |$0<\gamma \ll 1$|⁠, and setting |$C=0$| for example, the necessary and sufficient condition for the oscillatory motions to occur becomes

$$\begin{equation}-\infty=\lim_{x\to0}\frac{\mu}{\rho A^2}\frac{\left(x+1\right)\left(2x^4-x^2-1\right)}{x^3\left(x^3+x^2+x\right)}<\frac{p_{0}}{\gamma}<\sup\frac{\mu}{\rho A^2}\frac{\left(x+1\right)\left(2x^4-x^2-1\right)}{x^3\left(x^3+x^2+x\right)}\approx 0.7414\frac{\mu}{\rho A^2}.\end{equation}$$

(132)

where sup denotes supremum.

Hence, for the motion to be oscillatory, the shear modulus must be uniformly bounded from below as follows,

$$\begin{equation}\mu>\frac{p_{0}}{\gamma}\frac{\rho A^2}{0.7414}\approx\frac{2.7}{\gamma}\left(T_{1}(t)-T_{2}(t)\right).\end{equation}$$

(133)

Then, the probability distribution of oscillatory motions occurring is

$$\begin{equation}P_{1}(p_0/\gamma)=1-\int_{0}^{\frac{p_{0}}{\gamma}\frac{\rho A^2}{0.7414}}g(u;\rho_{1},\rho_{2})\,\mathrm{d}u,\end{equation}$$

(134)

and that of non-oscillatory motions is

$$\begin{equation}P_{2}(p_0/\gamma)=1-P_{1}(p_0/\gamma)=\int_{0}^{\frac{p_{0}}{\gamma}\frac{\rho A^2}{0.7414}}g(u;\rho_{1},\rho_{2})\,\mathrm{d}u.\end{equation}$$

(135)

For |$\rho =1$|⁠, |$A=1$| and |$\widetilde{\mu }=\mu =\mu _1+\mu _{2}$| drawn from the Gamma distribution with |$\rho _{1}=405$| and |$\rho _{2}=0.01$|⁠, the probability distributions given by (134) and (135) are shown in Fig. 15 (blue lines for |$P_{1}$| and red lines for |$P_{2}$|⁠). For the deterministic thin-walled tube, the critical value |$p_{0}/\gamma =0.7414\underline{\mu }=3.0027$| strictly separates the cases of oscillations occurring or not. However, in the stochastic case, the two cases compete.

Fig. 15.

Probability distributions of whether oscillatory motions can occur or not for a thin-walled spherical shell of stochastic neo-Hookean material, with |$\rho =1$|⁠, |$A=1$| and the shear modulus, |$\mu$|⁠, following the Gamma distribution with |$\rho _{1}=405$|⁠, |$\rho _{2}=0.01$|⁠. Dark-coloured lines represent analytically derived solutions, given by equations (134) and (135), whereas the lighter versions represent stochastically generated data. The vertical line at the critical value, |$p_{0}/\gamma =3.0027$|⁠, separates the expected regions based only on the mean value of the shear modulus, |$\underline{\mu }=\rho _{1}\rho _{2}=4.05$|⁠. The probabilities were calculated from the average of 100 stochastic simulations.

6. Conclusion

We provided here a synthesis on the analysis of finite amplitude oscillations resulting from dynamic finite deformations of given isotropic incompressible nonlinear hyperelastic solids and extended this to non-deterministic oscillatory motions of stochastic isotropic incompressible hyperelastic solids with similar geometries. Specifically, we treated in a unified manner the generalized shear motion of a cuboid and the radial motion of inflated cylindrical and spherical shells of stochastic neo-Hookean or Mooney–Rivlin material. For these finite dynamic problems, attention was given to the periodic motion and the time-dependent stresses, while taking into account the stochastic model parameters, which are random variables described by given probability laws. We found that, in this case, the amplitude and period of the oscillation of the stochastic bodies are also characterized by probability distributions, and, for cylindrical tubes and spherical shells, when an impulse surface traction is applied, there is a parameter interval where both the oscillatory and non-oscillatory motions can occur with a given probability. This is in contrast to the deterministic problem where a single critical parameter value strictly separates the cases where oscillations can or cannot occur.

The finite dynamic analysis presented here can be extended (albeit numerically) to other stochastic homogeneous hyperelastic materials (for example, using the stochastic strain-energy functions derived from experimental data in Mihai et al., 2018c), or to inhomogeneous incompressible bodies similar to those considered deterministically in Ertepinar & Akay (1976). For incompressible bodies with inhomogeneous material parameters, the constitutive parameters of the stochastic hyperelastic models can be treated as random fields, as described in Staber & Guilleminot (2018)and Staber et al. (2019). Clearly, the combination of knowledge from elasticity, statistics and probability theories offers a richer set of tools compared to the elastic framework alone and would logically open the way to further considerations of this type. However, as the role of stochastic effects on instabilities in finite strain elastodynamics is still in its infancy, it is important to consider the homogeneous case in the first instance.

If the material is compressible (unconstrained), then the theorem on quasi-equilibrated dynamics given by Truesdell (1962), and recalled by us in Section 2, is not applicable (Truesdell & Noll, 2004, p. 209). As we relied on the notion of quasi-equilibrated motion to derive our analytical results for incompressible cylindrical tubes and spherical shells, the same approach cannot be used for the compressible case. Nevertheless, as seen from the generalized shear motion of a cuboid, presented in Section 3, more general elastodynamic problems can still be formulated where the motion is not quasi-equilibrated. However, while stochastic versions of compressible hyperelastic materials can also be obtained, as shown in Staber & Guilleminot (2016), few theoretical results are available on the oscillatory motion of finitely deformed compressible hyperelastic solids (see, e.g., Akyüz & Ertepinar, 1998).

The analysis presented here is timely not only because ‘Today, it is well understood that as soon as the probability theory can be used, then the probabilistic approach of uncertainties is certainly the most powerful, efficient and effective tool for modelling and for solving direct and inverse problems’ (Soize, 2013), but also because time-dependent finite elastic deformations, although relevant to the modelling of various physical systems, have seldom been considered in more recent studies, which focused primarily on static elastic deformations or on dynamic viscoelasticity problems. Clearly, further numerical and experimental investigations of oscillatory finite deformations could help to bridge the gap between these popular areas and add some valuable insight into specific applications as well.

Funding

Engineering and Physical Sciences Research Council of Great Britain (EP/R020205/1 to A.G. and EP/S028870/1 to L.A.M.).

A. Additional detailed calculations

In this appendix, for the stochastic cylindrical and spherical shells discussed in Sections 4 and 5, respectively, we provide detailed derivations of the general functions |$G(x,\gamma )$|⁠, defined by (73), and |$H(x,\gamma )$|⁠, defined by (123), and calculate the limits of these functions in the particular cases of thin-walled and infinitely thick-walled shells.

• (I) For a Mooney-type model, the function |$G(x,\gamma )$| is defined by (65), where |$\widetilde{\mu }=\mu _{1}+\mu _{2}\alpha ^2$|t. In this case, we obtain

  $$\begin{equation*}\begin{split} G(x,\gamma)&=\frac{1}{\rho A^2}\int_{1/\sqrt{\alpha}}^{x}\left(\zeta\int_{\frac{\zeta^2+\frac{\gamma}{\alpha}}{1+\gamma}}^{\zeta^2}\widetilde{\mu}\frac{1+\alpha u}{\alpha^2u^2}\,\mathrm{d}u\right)\,\mathrm{d}\zeta\\ &=\frac{\widetilde{\mu}}{\rho A^2}\int_{1/\sqrt{\alpha}}^{x}\left(\zeta\int_{\frac{\zeta^2+\frac{\gamma}{\alpha}}{1+\gamma}}^{\zeta^2}\frac{1+\alpha u}{\alpha^2u^2}\,\mathrm{d}u\right)\,\mathrm{d}\zeta\\ &=\frac{\widetilde{\mu}}{\rho A^2}\int_{1/\sqrt{\alpha}}^{x}\left\{\frac{1}{\alpha^2}\left[(1+\gamma)\frac{\zeta}{\zeta^2+\frac{\gamma}{\alpha}}-\frac{1}{\zeta}\right]+\frac{1}{\alpha}\left[\zeta\log\zeta^2-\zeta\log\frac{\zeta^2+\frac{\gamma}{\alpha}}{1+\gamma}\right]\right\}\,\mathrm{d}\zeta\\ &=\frac{\widetilde{\mu}}{2\alpha\rho A^2}\left(\frac{1+\gamma}{\alpha}\log\frac{x^2+\frac{\gamma}{\alpha}}{\frac{1}{\alpha}+\frac{\gamma}{\alpha}} -\frac{1}{\alpha}\log x^2+\frac{1}{\alpha}\log\frac{1}{\alpha}\right)\\ &+\frac{\widetilde{\mu}}{2\alpha\rho A^2}\left(x^2\log x^2-x^2-\frac{1}{\alpha}\log\frac{1}{\alpha}+\frac{1}{\alpha}\right)\\ &-\frac{\widetilde{\mu}}{2\alpha\rho A^2}\left(x^2\log\frac{x^2+\frac{\gamma}{\alpha}}{1+\gamma}-x^2+\frac{\gamma}{\alpha}\log\frac{x^2+\frac{\gamma}{\alpha}}{\frac{1}{\alpha}+\frac{\gamma}{\alpha}}-\frac{1}{\alpha}\log\frac{1}{\alpha}+\frac{1}{\alpha}\right)\\ &=\frac{\widetilde{\mu}}{2\alpha\rho A^2}\left(x^2-\frac{1}{\alpha}\right)\log\frac{1+\gamma}{1+\frac{\gamma}{\alpha x^2}}.\end{split}\end{equation*}$$

  For the thin-walled tube (Knowles, 1962; Shahinpoor & Nowinski, 1971), |$\alpha =1$| and |$0<\gamma \ll 1$|⁠, and approximating |$\log (1+\gamma )$| by |$\gamma$| and |$\log \left [1+\gamma /\left (\alpha x^2\right )\right ]$| by |$\gamma /\left (\alpha x^2\right )$|⁠, we find

  $$\begin{equation*}G(x,\gamma)=\gamma\frac{\widetilde{\mu}}{2\rho A^2}\left(x^2-1\right)\left(1-\frac{1}{x^2}\right).\end{equation*}$$

  For the cylindrical cavity (Shahinpoor, 1973), |$\gamma \to \infty ;$| hence,

  $$\begin{equation*}G(x,\gamma)=\frac{\widetilde{\mu}}{2\alpha\rho A^2}\left(x^2-\frac{1}{\alpha}\right)\log\left(\alpha x^2\right).\end{equation*}$$
• (II) For a neo-Hookean-type model, the function |$H(x,\gamma )$| is defined by (116), where |$\widetilde{\mu }=\mu$|⁠. Following Knowles & Jakub (1965), we set the corresponding strain-energy density in the form

  $$\begin{equation*}W_{0}(u)=\frac{\mu}{2}\left(u^{-4/3}+2u^{2/3}-3\right),\end{equation*}$$

  and denote by |$W_{0}^{\prime}(u)$| its first derivative with respect to |$u$|⁠. Then by standard calculations (involving integration by parts and change of variables), we obtain

  $$\begin{equation*}\begin{split} H(x,\gamma)&=\frac{4}{3\rho A^2}\int_{1}^{x}\left(\zeta^2\int_{\frac{\zeta^3+\gamma}{1+\gamma}}^{\zeta^3}\widetilde{\mu}\frac{1+u}{u^{7/3}}\,\mathrm{d}u\right)\,\mathrm{d}\zeta\\[3pt] &=\frac{2}{\rho A^2}\int_{1}^{x}\left(\zeta^2\int_{\frac{\zeta^3+\gamma}{1+\gamma}}^{\zeta^3}\frac{W_{0}^{\prime}(u)}{u-1}\,\mathrm{d}u\right)\,\mathrm{d}\zeta\\[3pt] &=\frac{2}{\rho A^2}\int_{1}^{x}\left\{\zeta^2\left[\frac{W_{0}(\zeta^3)}{\zeta^3-1}-\frac{W_{0}\left(\frac{\zeta^3+\gamma}{1+\gamma}\right)}{\frac{\zeta^3+\gamma}{1+\gamma}-1}+\int_{\frac{\zeta^3+\gamma}{1+\gamma}}^{\zeta^3}\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u\right]\right\}\,\mathrm{d}\zeta\\[3pt] &=\frac{2}{\rho A^2}\left\{\int_{1}^{x}\zeta^2\left[\frac{W_{0}(\zeta^3)}{\zeta^3-1}-\frac{W_{0}\left(\frac{\zeta^3+\gamma}{1+\gamma}\right)}{\frac{\zeta^3+\gamma}{1+\gamma}-1}\right]\,\mathrm{d}\zeta +\int_{1}^{x}\left[\zeta^2\int_{\frac{\zeta^3+\gamma}{1+\gamma}}^{\zeta^3}\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u\right]\,\mathrm{d}\zeta\right\}\\[3pt] &=\frac{2}{3\rho A^2}\left[\int_{1}^{x^3}\frac{W_{0}(u)}{u-1}\,\mathrm{d}u+x^3\int_{1}^{x^3}\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u-\int_{1}^{x^3}\frac{uW_{0}(u)}{(u-1)^2}\,\mathrm{d}u\right]\\[3pt] &+\frac{2}{3\rho A^2}\left[\int_{\frac{x^3+\gamma}{1+\gamma}}^{1}(1+\gamma)\frac{W_{0}(u)}{u-1}\,\mathrm{d}u+x^3\int_{\frac{x^3+\gamma}{1+\gamma}}^{1}\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u-\int_{\frac{x^3+\gamma}{1+\gamma}}^{1} \right. \\[3pt] & \left.\frac{\left[u(1+\gamma)-\gamma\right]W_{0}(u)}{(u-1)^2}\,\mathrm{d}u\right]\\[3pt] &=\frac{2}{3\rho A^2}\left[x^3\int_{1}^{x^3}\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u-\int_{1}^{x^3}\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u\right]\\[3pt] &+\frac{2}{3\rho A^2}\left[x^3\int_{\frac{x^3+\gamma}{1+\gamma}}^{1}\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u-\int_{\frac{x^3+\gamma}{1+\gamma}}^{1}(1+\gamma)\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u+\int_{\frac{x^3+\gamma}{1+\gamma}}^{1}\gamma\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u\right]\\[3pt] &=\frac{2}{3\rho A^2}\left(x^3-1\right)\int_{\frac{x^3+\gamma}{1+\gamma}}^{x^3}\frac{W_{0}(u)}{(u-1)^2}\,\mathrm{d}u\\[3pt] &=\frac{\mu}{3\rho A^2}\left(x^3-1\right)\int_{\frac{x^3+\gamma}{1+\gamma}}^{x^3}\frac{2u^{4/3}+4u+3u^{2/3}+2u^{1/3}+1}{u^{2/3}(u+u^{2/3}+u^{1/3})^2}\,\mathrm{d}u\\[3pt] &=\frac{\mu}{\rho A^2}\left(x^3-1\right)\left[\frac{2x^3-1}{x^3+x^2+x} -\frac{2\frac{x^3+\gamma}{1+\gamma}-1}{\frac{x^3+\gamma}{1+\gamma}+\left(\frac{x^3+\gamma}{1+\gamma}\right)^{2/3}+\left(\frac{x^3+\gamma}{1+\gamma}\right)^{1/3}}\right].\end{split}\end{equation*}$$

  For the thin-walled shell (Beatty, 2011; Verron et al., 1999; Wang, 1965), |$0<\gamma \ll 1$|⁠, and

  $$\begin{equation*}\begin{split} H(x,\gamma)&=\gamma\frac{4\mu}{3\rho A^2}\int_{1}^{x}\frac{u^6-1}{u^5}\,\mathrm{d}u\\ &=\gamma\frac{\mu}{\rho A^2}\frac{\left(x+1\right)\left(2x^4-x^2-1\right)}{x^3\left(x^3+x^2+x\right)}.\end{split}\end{equation*}$$

  For the spherical cavity (Balakrishnan & Shahinpoor, 1978; Knowles & Jakub, 1965), |$\gamma \to \infty$|⁠; hence

  $$\begin{equation*}H(x,\gamma)=\frac{\mu}{3\rho A^2}\left(x^3-1\right)\frac{5x^3-x^2-x-3}{x^3+x^2+x}.\end{equation*}$$