Likely oscillatory motions of stochastic hyperelastic solids

# Likely oscillatory motions of stochastic hyperelastic solids

L. Angela Mihai111School of Mathematics, Cardiff University, Senghennydd Road, Cardiff, CF24 4AG, UK, Email: MihaiLA@cardiff.ac.uk   Danielle Fitt222School of Mathematics, Cardiff University, Senghennydd Road, Cardiff, CF24 4AG, UK, Email: FittD@cardiff.ac.uk   Thomas E. Woolley333School of Mathematics, Cardiff University, Senghennydd Road, Cardiff, CF24 4AG, UK, Email: WoolleyT1@cardiff.ac.uk   Alain Goriely444Mathematical Institute, University of Oxford, Woodstock Road, Oxford, OX2 6GG, UK, Email: goriely@maths.ox.ac.uk
###### Abstract

Stochastic homogeneous hyperelastic solids are characterised by strain-energy functions where the parameters are random variables satisfying probability density functions. These models can propagate 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. 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”.

Key words: stochastic hyperelastic models, finite dynamic deformations, quasi-equilibrated motion, finite amplitude oscillations, incompressibility, applied probability.

“Denominetur motus talis, qualis omni momento temporis 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.” - C. Truesdell (1962) [101]

## 1 Introduction

Motivated by numerous long-standing and modern engineering problems, oscillatory motions of cylindrical and spherical shells made of linear elastic material [55, 56, 53, 76] have generated a wide range of experimental, theoretical, and computational studies [5, 6, 7, 17, 27]. 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 [3, 9, 26, 38, 39, 40, 54], have been less investigated, and much of the work in finite nonlinear elasticity has focused on the static stability of pressurised shells [2, 16, 20, 21, 23, 29, 32, 34, 36, 57, 68, 79, 88, 109], or on wave-type solutions in infinite media [44, 75].

More recently, the dynamic response of incompressible hyperelastic cylindrical and spherical shells subjected to periodic loading was discussed in [77, 78]. Radial oscillations of cylindrical tubes and spherical shells of neo-Hookean [100], Mooney-Rivlin [67, 80], and Gent [30] hyperelastic materials were analysed in [14, 15], where it was concluded 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 [8, 10, 82, 108], where the results for Yeoh [107] and Mooney-Rivlin material models were compared. In [18], the static and dynamic behaviour of circular cylindrical shells of homogeneous isotropic incompressible hyperelastic material modelling arterial walls were considered.

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 [31, 42, 49, 72, 74, 99]. 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 characterised by probability density functions, as constructed in [94, 95, 96, 97, 98] and [66]. These models rely on the notion of entropy (or uncertainty) [87, 89] and on the maximum entropy principle for a discrete probability distribution [45, 46, 47], and are able to propagate uncertainties from input data to output quantities of interest [92]. They are also suitable for incorporation into Bayesian methodologies [13] (see also [60]) for models selection and updates [66, 72, 81].

To study the effect of probabilistic model parameters on predicted mechanical responses, in [62, 63, 64, 65], 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. Specific case studies included the classical problems of the Rivlin cube [65], the cavitation of a sphere under uniform tensile dead load [63], the inflation of pressurised spherical and cylindrical shells [62], and the rotation and perversion of anisotropic hyperelastic cylindrical tubes [64].

In this paper, we extend the stochastic framework developed in [63, 62, 64, 65] to the study of 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, and also of 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. 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”.

In Section 2, we provide a summary of the stochastic elasticity prerequisites. Section 3 is devoted to the analysis of the oscillatory motion of a stochastic hyperelastic cuboid under dynamic generalised shear. This is followed in Sections 4 and 5 by the treatment of the radial oscillatory motions of stochastic cylindrical tubes and spherical shells with bounded wall thickness respectively. The limiting cases of thin- and infinitely thick-walled structures are also discussed. Less straight-forward detailed calculations inherent for these problems are deferred to the Appendix A. Concluding remarks are drawn in Section 6.

## 2 Prerequisites

In this Section, we first recall the notion of (universal) quasi-equilibrated motion in finite elasticity, introduced in [101] and reviewed in [102], then summarise the stochastic finite elasticity framework developed in [66] and applied to static stability problems in [63, 62, 64, 65].

### 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 [102, p. 40],

 ρ¨x=div T+ρb, (1) T=TT, (2)

where is the motion of the given elastic solid, is the material density, which is assumed constant, is the body force, is the Cauchy stress tensor, and the superscript defines the transpose. To obtain possible dynamical solutions, one can solve Cauchy’s equation for particular motion or generalise known static solutions to the dynamical solution, using the notion of quasi-equilibrated motion:

###### Definition 2.1

[102, p. 208] A quasi-equilibrated motion is a motion, , of an incompressible homogeneous elastic solid subject to a given body force, , such that, for each value of , defines a static deformation that satisfies the equilibrium conditions under the body force .

###### Theorem 2.2

[102, p. 208] A quasi-equilibrated motion, , of an incompressible homogeneous elastic solid subject to a given body force, , is dynamically possible, subject to the same body force, if and only if the motion is circulation preserving with a single-valued acceleration potential , i.e.,

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

 curl ¨x=0. (4)

Then, the Cauchy stress tensor takes the form

 T=−ρξI+T(0), (5)

where is the Cauchy stress for the equilibrium state at time and 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 .

### 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 that satisfy standard probability laws [66, 94, 95, 96]. In this case, each model parameter is described in terms of its mean value and its variance, which contains information about the range of values about the mean value. As obtaining complete information about a random quantity in an elastic sample of material is rarely possible, the partial information provided by the mean value and the variance is commonly used in practice [19, 24, 42, 59, 69]. Here, we combine finite elasticity [33, 73, 102] and probability theory [37, 47], and rely on the following general hypotheses [66, 65, 63, 62]:

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

• Material isotropy, requiring that the strain-energy function is unaffected by a superimposed rigid-body transformation prior to deformation, i.e., , where . For isotropic materials, the strain-energy function is a symmetric function of the principal stretches of F, i.e., .

• Baker-Ericksen (BE) inequalities, which state that the greater principal (Cauchy) stress occurs in the direction of the greater principal stretch, are [11]

 (Ti−Tj)(λi−λj)>0ifλi≠λj,i,j=1,2,3, (6)

where and denote the principal stretches and the principal Cauchy stresses, respectively. The BE inequalities (6) take the equivalent form

 (λ1∂W∂λ1−λ2∂W∂λ2)(λ1−λ2)>0ifλi≠λj,i,j=1,2,3. (7)

In (6)-(7), the strict inequality “” is replaced by “” if any two principal stretches are equal.

• Finite mean and variance for the random shear modulus, i.e., at any given finite strain deformation, the random shear modulus, , and its inverse, , are second-order random variables [94, 95, 96].

Assumptions (A1)-(A3) are well-known principles in isotropic finite elasticity [33, 73, 102]. In particular, regarding the assumption (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 hold [58]. Another important deformation is that of simple shear superposed on axial stretch, defined by

 x1=αX1+kX2α2,x2=X2α2,x3=αX3, (8)

where and are the Cartesian coordinates for the reference (Lagrangian) and the current (Eulerian) configuration, respectively, and and are positive constants representing the shear parameter and the axial stretch ( for axial tension and for axial compression) respectively. For this deformation, the principal stretches, , satisfy

 λ21=α6+k2+1+√(α6+k2+1)2−4α62α4,λ22=α6+k2+1−√(α6+k2+1)2−4α62α4,λ23=α2, (9)

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

 Ti=λi∂W∂λi−p,i=1,2,3, (10)

where is the Lagrange multiplier for the incompressibility constraint. Then, if the BE inequalities hold, the nonlinear shear modulus defined by [61, 66]

 ˜μ=T1−T2λ21−λ22 (11)

is positive, i.e., , for all and . In the linear elastic limit, where and , the nonlinear shear modulus given by (11) converges to the classical shear modulus under infinitesimal deformation, i.e.,

 lima→1,k→0˜μ=μ. (12)

Assumption (A4) then contains physically realistic expectations on the (positive) random shear modulus , which will be characterised by a suitable probability density function.

In the next Sections, we analyse the dynamic generalised 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 known 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 [94, 66],

 W(λ1,λ2,λ3)=μ12(λ21+λ22+λ23−3)+μ22(λ−21+λ−22+λ−23−3), (13)

where and are random variables. The non-deterministic model (13) reduces to a stochastic neo-Hookean model if . If the random parameters and are replaced by their respective mean values, and , then the resulting mean hyperelastic model coincides with the usual deterministic one.

For the stochastic material model described by (13), the shear modulus in infinitesimal deformation is defined as . For this modulus, we set the following mathematical constraints, which are consistent with the assumption (A4) made in Section 2 [66],

 {E[μ]=μ––>0,E[log μ]=ν,such that |ν|<+∞, (14)

i.e., the mean value of the shear modulus is fixed and greater than zero, and the mean value of is fixed and finite. It follows that and are second-order random variables, i.e., they have finite mean and finite variance [90, 91]. Under the constraints (14), follows a Gamma probability distribution with hyperparameters and satisfying

 μ––=ρ1ρ2,Var[μ]=ρ1ρ22, (15)

where is the mean value, is the variance , and is the standard deviation of . The corresponding probability density function takes the form [1, 48]

 g(μ;ρ1,ρ2)=μρ1−1e−μ/ρ2ρρ12Γ(ρ1),for μ>0 and ρ1,ρ2>0, (16)

where is the complete Gamma function

 Γ(z)=∫+∞0tz−1e−tdt. (17)

When and , we can define the auxiliary random variable [66]

 R1=μ1μ, (18)

such that . Then

 μ1=μR1,μ2=μ−μ1=μ(1−R1). (19)

Setting the realistic constraints [94, 95, 96, 66],

 (20)

follows a standard Beta distribution [1, 48], with hyperparameters and satisfying

 R––1=ξ1ξ1+ξ2,Var[R1]=ξ1ξ2(ξ1+ξ2)2(ξ1+ξ2+1). (21)

where is the mean value, is the variance, and is the standard deviation of . The corresponding probability density function takes the form

 β(r;ξ1,ξ2)=rξ1−1(1−r)ξ2−1B(ξ1,ξ2),for r∈(0,1) and ξ1,ξ2>0, (22)

where is the Beta function

 B(x,y)=∫10tx−1(1−t)y−1dt. (23)

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

 μ––1=μ––R––1,μ––2=μ––−μ––1=μ––(1−R––1), (24)

and the variances and covariance are, respectively,

 Var[μ1]=(μ––)2Var[R1]+(R––1)2Var[μ]+Var[μ]Var[R1], (25) Var[μ2]=(μ––)2Var[R1]+(1−R––1)2Var[μ]+Var[μ]Var[R1], (26) Cov[μ1,μ2]=12(Var[μ]−Var[μ1]−Var[μ2]). (27)

Note that the random variables and are independent, depending on parameters and , respectively, whereas and are dependent variables as they both require to be defined. Explicit derivations of stochastic isotropic hyperelastic models calibrated to experimental data are presented in [66].

Numerically, throughout this paper, we assume that the random shear modulus follows the Gamma distribution represented in Figure 1, with the shape and scale parameters and respectively [64, 63, 62].

## 3 Generalised shear motion of a stochastic hyperelastic cuboid

The generalised shear motion of a cuboid of elastic material is described by [26]

 x=X√α,y=Y√α,z=αZ+u(X,Y,t), (28)

where and are the Cartesian coordinates for the reference (Lagrangian, material) and current (Eulerian, spatial) configuration respectively, is a given constant, and , 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.

By the governing equations (28), the condition (4) is valid for if and only if

 0=curl ¨x=⎡⎢⎣∂¨z/∂y−∂¨y/∂z∂¨x/∂z−∂¨z/∂x∂¨y/∂x−∂¨x/∂y⎤⎥⎦=⎡⎢⎣∂¨u/∂Y−∂¨u/∂X0⎤⎥⎦. (29)

This condition imposes very strict constraints on the motion. Yet, we will see that even though the generalised shear motion (28) is not quasi-equilibrated, exact solutions are still available, although these solutions are not universal [70, 104].

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

 F=⎡⎢⎣1/√α0001/√α0uXuYα⎤⎥⎦,

where and denote the partial first derivatives of with respect to and respectively. The corresponding left Cauchy-Green tensor is

 (30)

and has the principal invariants

 I1=tr (B)=u2X+u2Y+2α+α2,I2=12[(tr% B)2−tr(B2)]=u2Xα+u2Yα+1α2+2α,I3=detB=1. (31)

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

 T=−pI+β1B+β−1B−1, (32)

where is the Lagrange multiplier for the incompressibility constraint (), and

 β1=2∂W∂I1,β−1=−2∂W∂I2 (33)

are the nonlinear material parameters, with , given by (31).

### 3.1 Shear oscillations of a cuboid of stochastic neo-Hookean material

We now specialise to the case of a cuboid of stochastic neo-Hookean material, with and in (13), where the non-zero components of the Cauchy stress tensor given by (32) are as follows

 Txx=Tyy=−p+μα,Tzz=−p+μ(u2X+u2Y+α2),Txz=μ√αuX,Tyz=μ√αuY. (34)

Then, by the equation of motion (1),

 ∂p∂x=0,∂p∂y=0,∂p∂z=−ρ¨u+μ(uXX+uYY), (35)

where and represent the second derivatives of with respect to and respectively. Hence, is independent of and .

We consider the undeformed cuboid to be long in the -direction, and impose an initial displacement and velocity . For the boundary condition, we distinguish the following two cases:

#### (i)

If we impose null normal Cauchy stresses, , on the faces perpendicular to the - and -directions, at all time, then is constant and .

#### (ii)

If , as cannot be made point-wise zero, we denote the normal force acting on the cross-sections of area in the -direction at time by

 Nz(t)=∫ATzzdA, (36)

and consider this force to be zero, i.e., at all time. Then, is independent of , and, by (35), it is also independent of and , hence, .

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

 ¨u=μρ(uXX+uYY). (37)

It remains to solve, by standard procedures, the linear wave equation (37), describing the propagation of waves, subject to the given initial and boundary conditions. To solve this equation, we let the shear stresses and , defined by (34), vanish at the sides, i.e.,

 Txz(0,Y,Z,t)=Txz(1,Y,Z,t)=0⟺uX(0,Y,t)=uX(1,Y,t)=0,Tyz(X,0,Z,t)=Tyz(X,1,Z,t)=0⟺uY(X,0,t)=uY(X,1,t)=0. (38)

In this case, the general solution takes the form

 u(X,Y,t)=∞∑m=1∞∑n=1[Amncos(ωmnt)+Bmnsin(ωmnt)]cos(πmX)cos(πnY), (39)

where

 ωmn=π√(m2+n2)μρ, (40)

and

 Amn=4∫10[∫10u0(X,Y)cos(πmX)dX]cos(πnY)dY, (41) Bmn=4ωmn∫10[∫10˙u0(X,Y)cos(πmX)dX]cos(πnY)dY. (42)

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

As is a random variable, it follows that the speed of wave propagation, , is stochastic. Hence, both the period and the amplitude of the oscillations are stochastic. As an example, we consider the initial data and leading to and . In Figure 3, we illustrate the stochastic dynamic displacement on the edges when , , , , and is drawn from the Gamma distribution with hyperparameters and , as represented in Figure 1. We note from Figure 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 in the top of Figure 3 and by the tight distribution around the mean in the bottom of Figure 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 of Figure 3).

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

For a circular cylindrical tube, the combined radial and axial motion is described by

 r2=a2+R2−A2α,θ=Θ,z=αZ, (43)

where and are the cylindrical polar coordinates in the reference and current configuration, respectively, such that , and are the inner and outer radii in the undeformed state respectively, and are the inner and outer radius at time respectively, and is a given constant (when , the tube is everted, so that the inner surface becomes the outer surface). When , the time-dependent deformation (43) simplifies to that studied also in [14, 50, 51]. The case when is time-dependent was considered in [84].

The radial-axial motion (43) of the cylindrical tube is fully determined by the inner radius at time , which in turn is obtained from the initial conditions. Thus, the acceleration can be computed in terms of the acceleration on the inner surface. By the governing equations (43), the condition (4) is valid for , since

 (44)

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

 −∂ξ∂r=¨r=˙a2r+a¨ar−a2˙a2r3, (45)

and, by integrating (45), the acceleration potential, , is given by [102, p. 215]

 −ξ=˙a2logr+a¨alogr+a2˙a22r2=˙r2logr+r¨rlogr+12˙r2. (46)

The deformation gradient of (43), with respect to the polar coordinates , is equal to

 F=diag(Rαr,rR,α), (47)

the Cauchy-Green deformation tensor is

 B=F2=diag(R2α2r2,r2R2,α2), (48)

and the principal invariants take the form

 I1=tr (B)=R2α2r2+r2R2+α2,I2=12[(tr% B)2−tr(B2)]=α2r2R2+R2r2+1α2,I3=detB=1. (49)

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

 T(0)rr=−p(0)+β1R2α2r2+β−1α2r2R2,T(0)θθ=T(0)rr+(β1−β−1α2)(r2R2−R2α2r2),T(0)zz=T(0)rr+(β1−β−1r2R2)(α2−R2α2r2), (50)

where is the Lagrangian multiplier for the incompressibility constraint (), and

 β1=2∂W∂I1,β−1=−2∂W∂I2 (51)

are the nonlinear material parameters, with and given by (49).

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

 ∂T(0)rr∂r=T(0)θθ−T(0)rrr. (52)

Hence, by (50) and (52), the radial Cauchy stress for the equilibrium state at time is equal to

 T(0)rr(r,t)=∫(β1−β−1α2)(r2R2−R2α2r2)drr+ψ(t), (53)

where is an arbitrary function of time. Substitution of (46) and (53) into (5) then gives the principal Cauchy stress components at time as follows,

 (54)

In (54), the function can be interpreted as the following nonlinear shear modulus [61]

 ˜μ=β1−β−1α2, (55)

corresponding to the combined deformation of simple shear superposed on axial stretch, described by (8), with shear parameter and stretch parameter . As shown in [61], this modulus is positive if the BE inequalities hold. In this case, the integrand is negative for and positive for . Using the first equation in (43), it is straightforward to show that (respectively, ) is equivalent to (respectively, ). When , the modulus defined by (55) coincides with the generalised shear modulus defined in [102, p. 174], and also in [14].

In the limiting case when and , the nonlinear shear modulus defined by (55) converges to the classical shear modulus from the infinitesimal theory [61],

 μ=limα→1limk→0˜μ. (56)

In this case, as , the three stress components defined by (54) are equal.

Next, for the cylindrical tube deforming by (43), we set the inner and outer radial pressures acting on the curvilinear surfaces and at time (measured per unit area in the present configuration), as and respectively [102, pp. 214-217]. Evaluating and , using (54), with and respectively, then subtracting the results, then gives

 (57)

Setting the notation

 u=r2R2=r2α(r2−a2)+A2,x=aA,γ=B2A2−1, (58)

we can rewrite

 (aA¨aA+˙a2A2)logb2a2+˙a2A2(a2b2−1)=(¨xx+˙x2)log(1+γαx2)−˙x2γαx21+γαx2=12xddx[˙x2x2log(1+γαx2)]

and

 ∫ba˜μ(r2R2−R2α2r2)drr=∫ba˜μ[r2α(r2−a2)+A2−α(r2−a2)+A2α2r2]drr=12∫x2x2+γα1+γ˜μ1+αuα2u2du.

Then, we can express the equation (57) equivalently as follows,

 2xT1(t)−T2(t)ρA2=12ddx[˙x2x2log(1+γαx2)]+xρA2∫x2x2+γα1+γ˜μ1+αuα2u2du. (59)

Note that, when the BE inequalities hold, , and the integral in (57), or equivalently in (59), is negative if (i.e., if ) and positive if (i.e., if ).

In the static case, where and , (57) becomes

 T1(t)−T2(t)=∫ba˜μ(r2R2−R2α2r2)drr, (60)

and (59) reduces to

 2T1(t)−T2(t)ρA2=1ρA2∫x2x2+γα1+γ˜μ1+αuα2u2du. (61)

For the cylindrical tube in finite dynamic deformation, we set

 G(x,γ)=1ρA2∫x1/√α⎡⎣ζ∫ζ2ζ2+γα1+γ˜μ1+αuα2u2du⎤⎦dζ, (62)

and find that is monotonically decreasing when and increasing when . 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)

 2αT1(t)−T2(t)ρA2={0% if t≤0,p0if t>0, (63)

where is constant in time. Then, integrating (59) once gives

 12˙x2x2log(1+γαx2)+G(x,γ)=p02α(x2−1α)+C, (64)

with defined by (62) and

 C=12˙x20x20log(1+γαx20)+G(x0,γ)−p02α(x20−1α), (65)

where and are the initial conditions. By (64),

 ˙x=±   ⎷p0α(x2−1α)+2C−2G(x,γ)x2log(1+γαx2). (66)

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

 E=12m(x)˙x2+V(x), (67)

The energy is , the potential is given by and the position-dependent mass is . Due to the constraints on the function , this system has simple dynamics. Depending on the constant , the system may have a static state or periodic motion. Indeed, the radial motion is periodic if and only if the following equation,

 G(x,γ)=p02α(x2−1α)+C, (68)

has exactly two distinct solutions, representing the amplitudes of the oscillation, and , such that . Then, by (58), the minimum and maximum radii of the inner surface in the oscillation are equal to and respectively, and by (66), the period of oscillation is equal to

 T=2∣∣∣∫x2x1dx˙x∣∣∣=2∣∣ ∣ ∣ ∣∣∫x2x1   ⎷x2log(1+γαx2)p0α(x2−1α)+2C−2G(x,γ)dx∣∣ ∣ ∣ ∣∣. (69)

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

### 4.1 Radial oscillations of a cylindrical tube of stochastic Mooney-Rivlin material

For cylindrical tubes of stochastic Mooney-Rivlin material defined by (13), with , evaluating the integral in (62) gives (see Appendix A for detailed calculations)

 G(x,γ)=˜μ2αρA2(x2−1α)log1+γ1+γαx2, (70)

where . In this case, assuming that the nonlinear shear modulus has a uniform lower bound, i.e.,

 μ>η, (71)

for some constant , it follows that

 limx→0G(x,γ)=limx→∞G(x,γ)=∞. (72)

#### (i)

If and , then equation (68) has exactly two solutions, and , satisfying , for any positive constant . It should be noted that, by (54), if at and , so that , then and at and , unless and