On deforming a sector of a circular cylindrical tube into an intact tube: existence, uniqueness, and stability
OOn deforming a sector of a circular cylindricaltube into an intact tube:existence, uniqueness, and stability
M. Destrade a,b , J.G. Murphy c , R.W. Ogden b,d,ea School of Mathematics, Statistics and Applied Mathematics,National University of Ireland Galway,University Road, Galway, Ireland b School of Electrical, Electronic, and Mechanical Engineering,University College Dublin,Belfield, Dublin 4, Ireland c Department of Mechanical Engineering,Dublin City University,Glasnevin, Dublin 9, Ireland d School of Engineering,University of Aberdeen,King’s College, Aberdeen AB24 3UE, United Kingdom e School of Mathematics and Statistics,University of Glasgow,University Gardens, Glasgow G12 8QW, Scotland, UK a r X i v : . [ n li n . PS ] J a n bstract Within the context of finite deformation elasticity theory the prob-lem of deforming an open sector of a thick-walled circular cylindricaltube into a complete circular cylindrical tube is analyzed. The analysisprovides a means of estimating the radial and circumferential residualstress present in an intact tube, which is a problem of particular con-cern in dealing with the mechanical response of arteries. The initialsector is assumed to be unstressed and the stress distribution result-ing from the closure of the sector is then calculated in the absence ofloads on the cylindrical surfaces. Conditions on the form of the elas-tic strain-energy function required for existence and uniqueness of thedeformed configuration are then examined. Finally, stability of the re-sulting finite deformation is analyzed using the theory of incrementaldeformations superimposed on the finite deformation, implemented interms of the Stroh formulation. The main results are that convex-ity of the strain energy as a function of a certain deformation vari-able ensures existence and uniqueness of the residually-stressed intacttube, and that bifurcation can occur in the closing of thick, widelyopened sectors, depending on the values of geometrical and physicalparameters. The results are illustrated for particular choices of theseparameters, based on data available in the biomechanics literature.
Keywords : nonlinear elasticity; finite elasticity; residual stress; elastic tube;existence and uniqueness; stability
Residual stresses have a very important role in the mechanical functioning ofarteries and the existence of residual stresses is well documented, as in, forexample, Chuong and Fung (1986), Vaishnav and Vossoughi (1983) and Fung(1993); see also the review by Humphrey (2002). They are demonstrated bythe simple device of cutting radially a short ring of an artery, the result of thecut being the springing open of the ring into a sector, thereby releasing someresidual stresses; in general, some residual stress remains, however, as shownby Vossoughi et al. (1993) and Greenwald et al. (1997); see also the reviewby Rachev and Greenwald (2003). A crude estimate of the circumferentialresidual stresses is obtained by measuring the resulting angle of the sectorinto which the ring deforms, although it has to be said that reliable quan-titative data remain elusive. This so-called opening angle method has been2nalyzed in some detail by several authors, including, for example, Delfinoet al. (1997), Zidi et al. (1999), Rachev and Hayashi (1999), Holzapfel etal. (2000), Ogden and Schulze-Bauer (2000), Matsumoto and Sato (2002),Ogden (2003), Raghavan et al. (2004) and Olsson et al. (2006). The typicalapproach is to assume that the (unloaded) opened sector is free of stressand circular cylindrical and then to construct a deformation that brings thesector into an intact tube and to calculate the resulting radial and circum-ferential (residual) stress distributions. This has been done for a single layerand for two-layered tubes. In general this method does not account for anystress or deformation in the axial direction. However, in reality, not onlydoes the ring open into a sector, but the length of the arterial segment tendsto change and there are also residual stresses in the axial direction, althoughquantitative information about axial residual stresses is distinctly lacking. Arecent 3D analysis by Holzapfel and Ogden (2010) is the first attempt to cal-culate the combination of radial, circumferential and axial residual stresses.The purpose of the present paper is to provide an analysis of the openingangle method, with particular reference to the questions of existence anduniqueness of the residually-stressed configuration and its stability.In Section 2 we review the description of the geometry of the deformationwhich first takes the sector into an intact cylinder, allowing for a possiblelength change, and then allows inflation under internal pressure combinedwith axial load while maintaining circular cylindrical symmetry. The corre-sponding stress components and equilibrium equations are then summarizedin respect of an isotropic elastic material whose properties are described interms of a strain-energy function. The existence and uniqueness of the re-sulting configuration is then examined in Section 3, under the assumptionthat the strain energy is convex as a function of a suitably chosen deforma-tion variable. The Mooney-Rivlin and the Fung strain energies obey thisassumption, and we use them, as well as some experimental data on arteries,to illustrate the analysis. Section 4 is devoted to an analysis of the stabilityof the residually stressed tube (in the absence of lateral loads) on the basis ofthe theory of small incremental deformations superimposed on a finite defor-mation, and for this purpose an appropriate version of the Stroh formulationis adopted. For simplicity of illustration, attention is restricted to prismaticincremental deformations. We then treat numerically the case of tubes madeof solids with the Mooney-Rivlin strain energy, and compare the predictionsto some simple experiments we performed with silicone rubber.3
Problem formulation
We consider an annular sector of a circular cylindrical tube with geometrydefined by A ≤ R ≤ B, − (2 π − α ) / ≤ Θ ≤ +(2 π − α ) / , ≤ Z ≤ L, (1)where A and B are the radii of the inner and outer faces, respectively, of thesector, L is its length, and ( R, Θ , Z ) denote the cylindrical coordinates ofa point in the material, with corresponding orthonormal basis vectors ( E R , E Θ , E Z ). Here α ∈ (0 , π ) is the so-called opening angle . The sector isassumed to be unstressed in this reference configuration, which we denote by B . The sector is then deformed into an intact (circular cylindrical) tube sothat the plane faces originally at Θ = ± ( π − α/
2) are joined perfectly to-gether and there is an accompanying uniform axial stretch λ z = l/L > l is the deformed length of the tube. We refer to this new configura-tion as the residually-stressed configuration, which we denote by B r . In thisconfiguration there is no traction on the curved surfaces of the tube, but ingeneral axial loads are required to maintain the deformation since it turnsout that the axial stress depends on the radius r and cannot therefore beset to zero pointwise on the ends of the tube. Finally, the tube is subjectedto a uniform internal pressure of magnitude P per unit deformed area andis maintained in a new circular cylindrical configuration, which we refer toas the loaded configuration , denoted B l ; the three separate configurations aredepicted in Figure 1.Let ( r, θ, z ) denote cylindrical coordinates in either the residually-stressedconfiguration B r or the pressurized configuration B l , with corresponding or-thonormal basis vectors ( e r , e θ , e z ). The deformation may be described bythe equations r = r ( R ) , θ = k Θ , z = λ z Z, (2)from which we may calculate the deformation gradient tensor F = Grad x as F = r (cid:48) e r ⊗ E R + ( kr/R ) e θ ⊗ E Θ + λ z e z ⊗ E Z , (3)where the prime denotes differentiation and k ≡ π π − α , k > , (4)4 Figure 1: Depiction of the three configurations: (a) unstressed open sector(configuration B ), with the opening angle α shown as π/
2; (b) residuallystressed and stretched tube (configuration B r ); (c) tube in (b) subject tointernal pressure P (configuration B l ).is a measure of the opening angle . The left Cauchy–Green tensor is thencalculated as B = F F T = r (cid:48) e r ⊗ e r + ( kr/R ) e θ ⊗ e θ + λ z e z ⊗ e z , (5)and it follows immediately that the principal axes of B are e r , e θ , and e z .The corresponding principal stretches are λ = r (cid:48) , λ = krR , λ = λ z . (6)5rterial wall tissue is generally considered to be essentially incompressible ,so that the constraint det F = 1 is enforced. For the present deformationthis constraint yields r = R − A kλ z + a , (7)where a ≡ r ( A ) is the radius of the inner surface of the tube in configuration B l . Similarly, we define b ≡ r ( B ) for the outer surface. It then follows fromfrom (6) that the principal stretches become λ = Rkλ z r , λ = krR , λ = λ z . (8)Notice from (7) that if a kλ z = A then the principal stretches are con-stants. In other words, the deformation is homogeneous , of the form r = R/ (cid:112) kλ z , θ = k Θ , z = λ z Z, (9)with principal stretches λ = 1 / √ kλ z , λ = (cid:112) k/λ z , λ = λ z . Henceforth we restrict the analysis to homogeneous, incompressible, isotropichyperelastic solids. The strain-energy function is then a symmetric functionof the principal stretches, which we write as W ( λ , λ , λ ), subject to theconstraint λ λ λ = 1. For such materials the Cauchy stress tensor σ iscoaxial with B and for the considered deformation the principal Cauchystress components are σ rr = − q + λ ∂W∂λ , σ θθ = − q + λ ∂W∂λ , σ zz = − q + λ ∂W∂λ , (10)where q is a Lagrange multiplier introduced by the constraint of incompress-ibility, and λ , λ , λ are given by (8).Because the stretches are independent of θ , and z the θ and z componentsof the equilibrium equation ensure that q is also independent of θ and z . Theonly non-trivial equation remaining is the radial equationdd r σ rr + 1 r ( σ rr − σ θθ ) = 0 . (11)6ith this equation we associate the boundary conditions σ rr ( a ) = − P, σ rr ( b ) = 0 , (12)where P is the internal pressure in configuration B l and the outer boundaryis free of traction.At this point it is convenient to introduce the notations x a ≡ kλ z a A , x ≡ kλ z r R , x b ≡ kλ z b B , ε ≡ B A − > , (13) ε being a measure of the thickness of the initial annular sector. Then itfollows easily that x b = ε + x a ε + 1 . (14)For any fixed value of the axial stretch λ z we may now regard W as afunction of x through the stretches, which are given in terms of x by λ = 1 √ kλ z x , λ = (cid:114) kxλ z , (15)together with λ = λ z , and we write (cid:99) W ( x ) = W (1 / √ kλ z x, (cid:112) kx/λ z , λ z ).From (10) it then follows that σ θθ − σ rr = λ ∂W∂λ − λ ∂W∂λ = 2 x (cid:99) W (cid:48) ( x ) . (16)Noting that r d x d r = 2 x (1 − x ) , (17)we may integrate equation (11) and apply the boundary conditions (12), andthus conclude that P = (cid:90) x b x a (cid:99) W (cid:48) ( x )1 − x d x. (18)This equation, together with (14), enables the location of the deformed innerradius to be determined for given values of the pressure P , the thickness (asmeasured by ε ), the opening angle (as measured via k ), and the axial stretch7 z . For completeness we compute the corresponding stress distribution, givenby σ rr = (cid:90) xx b (cid:99) W (cid:48) ( s )1 − s d s, σ θθ = σ rr + 2 x (cid:99) W (cid:48) ( x ) , σ zz = σ rr + λ ∂W∂λ − λ ∂W∂λ . (19)In the special case of the homogeneous deformation (9), we have x ≡ σ θθ − σ rr = 2 (cid:99) W (cid:48) (1) = constant , (20)so that we can integrate the equation of equilibrium (11) to give σ rr = 2 (cid:99) W (cid:48) (1) ln ( r/b ) , (21)where the boundary condition (12) has been used. The boundary condition(12) then yields P = (cid:99) W (cid:48) (1) ln ( b/a ) = (cid:99) W (cid:48) (1) ln ( B/A ) . (22)This is the pressure that needs to be applied inside the tube in order toproduce a homogeneous deformation through the tube wall. We note that inthis case, both σ rr , given by (21), and σ θθ , given by σ θθ = P [1 + ln( r/b )] / ln( b/a ) , (23)have logarithmic variations. Hence, in the case of homogeneous deformation,the stress is a slowly varying function of the radial coordinate.It is clear that in general σ zz depends on r . In the case of homogeneousdeformations this follows immediately from (21) and (19) . Thus, we em-phasize that the pointwise end condition σ zz = 0 cannot be adopted and theresidually stressed configuration B r is not load free. In some situations it maybe possible to enforce the condition that the resultant axial load vanishes,but here we merely consider that the axial stretch is fixed, so that the axialload has to be adjusted accordingly to accommodate this. Because we are considering an incompressible isotropic elastic material, thestrain-energy function can also be written as W = (cid:102) W ( I , I ), where I , I are8he first and second principal invariants of the Cauchy–Green deformationtensor B . Expressed in terms of the principal stretches these are I = λ + λ + λ , I = λ λ + λ λ + λ λ . (24)From (8) we then have explicitly I = (cid:18) Rkλ z r (cid:19) + (cid:18) krR (cid:19) + λ z = 1 kλ z x + kxλ z + λ z ,I = 1 λ z + (cid:18) Rkr (cid:19) + (cid:18) kλ z rR (cid:19) = λ z kx + kλ z x + λ − z . (25)The function (cid:99) W (cid:48) ( x ), occurring in the integrand of (18), can therefore bewritten as (cid:99) W (cid:48) ( x ) = k x − kλ z x (cid:32) ∂ (cid:102) W∂I + λ z ∂ (cid:102) W∂I (cid:33) . (26)Assuming that the standard empirical inequalities ∂ (cid:102) W /∂I > ∂ (cid:102) W /∂I ≥ (cid:102) W ( I , I ) (cid:99) W (cid:48) ( x ) = 0 ⇔ x = x m ≡ /k = (2 π − α ) / (2 π ) . (27)We assume here and henceforth that (cid:99) W ( x ) is a strictly convex function ,i.e. (cid:99) W (cid:48)(cid:48) ( x ) >
0. It follows that for any given value of λ z , x = x m ≡ /k yields the unique minimum of the strain energy function. Further, by (24)we have I ( x m ) = 2 λ − z + λ z , I ( x m ) = 2 λ z + λ − z , (28)and the minimum value of the strain energy is therefore independent of themeasure of the opening angle k ≡ π/ (2 π − α ).From equations (11), (16), and (17) we haved σ rr d x = (cid:99) W (cid:48) ( x )1 − x , (29)and consequently, the radial stress σ rr also has a unique extremum, a min-imum, at x = x m . The corresponding relation for the hoop stress σ θθ isderivable from (19) and has the formd σ θθ d x = (cid:99) W (cid:48) ( x )1 − x + 2 (cid:99) W (cid:48) ( x ) + 2 x (cid:99) W (cid:48)(cid:48) ( x ) . (30)9n contrast to the radial stress, the determination of its singular point(s)depends on the form of the strain energy function. Clearly, the hoop stress isincreasing at x = x m , the radius where the radial stress is at a minimum. Atthat point, according to (19) and (27), σ θθ = σ rr and, by (8), λ = λ ; thusthe in-plane principal stresses and the stretches are equal there. We refer tothe surface defined by x = x m as the neutral surface . The associated radius( r or R ) depends, in general, on both λ z and k .Another consequence of (26) and convexity is that (cid:40) (cid:99) W (cid:48) ( x ) < < x < x m ≡ /k, (cid:99) W (cid:48) ( x ) > x > x m ≡ /k. (31)Application of these conditions to the tube, however, depends on whether x a > x b or x b > x a . It is easy to show from (14) that either (i) x a > x b > > x b > x a . In case (i) we have x m = 1 /k < x m does not liewithin the required range and there is no neutral surface in this case. In case(ii) we require x b > x m > x a . Thus, (31) can be made more explicit: (cid:40) (cid:99) W (cid:48) ( x ) < x a < x < x m ≡ /k, (cid:99) W (cid:48) ( x ) > > x b > x > x m ≡ /k. (32)This will be used in the next section to determine existence and uniquenessof the unloaded, residually-stressed configuration.In closing this section, we note that convexity of (cid:99) W ( x ) is a common fea-ture of many standard strain-energy functions, including the Mooney-Rivlinmodel , (cid:102) W MR = C I −
3) + D I − , C > , D ≥ , (33)where the material constants C and D have dimensions of stress, the Gentmodel (cid:102) W G = − µJ m (cid:18) − I − J m (cid:19) , µ > , J m > , (34)and the Fung model (cid:102) W F = µ c (cid:2) e c ( I − − (cid:3) , µ > , c > , (35)where µ is the shear modulus in the unstressed configuration and J m and c are dimensionless constants. Convexity may easily be checked. For instance, (cid:99) W (cid:48)(cid:48) MR ( x ) = ( C + Dλ z ) / ( kλ z x ) which is clearly positive. Similar, but longer,expressions can be found for (cid:99) W (cid:48)(cid:48) G and (cid:99) W (cid:48)(cid:48) F , which are also strictly positive.10 Existence and uniqueness of the residually-stressed configuration
Here and henceforth, we study the residually-stressed tube in configuration B r . Setting P = 0 in (18) yields0 = (cid:90) x b x a (cid:99) W (cid:48) ( x )1 − x d x, x b = ε + x a ε + 1 . (36)This equation determines an inner radius of the unloaded residually-stressedtube, measured by x a , for a given measure k of the opening angle and givenaxial stretch λ z .Two questions need to be answered.(i) Does equation (36) have a positive solution x a ? In other words, is italways possible to close a cylindrical sector in the way indicated whenno load is applied to its curved faces?(ii) If a solution exists, is it unique?First consider the possible existence of a homogeneous residually-stressedconfiguration, when the deformation is in the form (9), which is equivalent toassuming that x a = 1. Then, by (22), vanishing of P is equivalent to (cid:99) W (1) =0. According to (27) this leads to k = 1, which is not admissible since k > there cannot be a homogeneous residually-stressed state of deformation ,i.e. a residually-stressed configuration is necessarily inhomogeneous. This,of course, is well known in general, as shown by, for example, Hoger (1985);see also the discussion in Ogden (2003).Assume now that x a >
1. Then, as indicated following (31), we musthave x a > x b > > x m , and hence, according to (31), the integrand in (36)has the same sign over the entire range of integration, and there is thereforeno solution to (36) for x a > x a of (36) to exist, it must lie in the range 0 < x a <
1. Ifeither x a < x b < x m or x m < x a < x b then there is again no solution to (36),because the integrand is one-signed over the range of integration. Thus, wemust have the ordering x a < x m < x b and (32) is applicable. Expressing x b in terms of x a by means of (36) in this latter inequality yields the followingrange of values for x a : [1 − ( k − ε ] x m < x a < x m . (37)11hen, to ensure that the solution x a to (36) is positive, we enforce the in-equality ε < / ( k − B A < kk − πα , (38)which prescribes a maximum ratio of outer to inner radius for the originalannular sector for a given opening angle, or equivalently a maximum valueof the opening angle for a given radius ratio.Now, addressing question (i) above is equivalent to determining whetherthe function f defined over the range (37) by f ( y ) = (cid:90) ( ε + y ) / ( ε +1) y (cid:99) W (cid:48) ( x )1 − x d x, [1 − ( k − ε ] x m < y < x m , (39)has a zero. First we examine the value of f ( y ) at the upper bound of therange, which is f ( x m ) = (cid:90) ( ε + x m ) / ( ε +1) x m (cid:99) W (cid:48) ( x )1 − x d x, (40)and this is positive because x m < ε + x m ε + 1 < , (41)so that the integrand in (40) is always positive over the range of integration.Next, we examine f ( y ) at the lower bound of the range (37); after somesimplification, we find that f ([1 − ( k − ε ] x m ) = (cid:90) x m [1 − ( k − ε ] x m (cid:99) W (cid:48) ( x )1 − x d x, (42)which is strictly negative because[1 − ( k − ε ] x m < x m < , (43)so that the integrand in (40) is negative over the whole range of integration.Because f is continuous we conclude that there exists at least one solution y = x a of the equation f ( y ) = 0 in the range (37).The answer to question (ii) can be found by computing f (cid:48) ( y ), which yields f (cid:48) ( y ) = (cid:99) W (cid:48) (( ε + y ) / ( ε + 1)) − (cid:99) W (cid:48) ( y )1 − y . (44)12ecause (cid:99) W (cid:48) is a monotonic increasing function, and because of the ordering y < ε + yε + 1 < , (45)it follows that f (cid:48) ( y ) > y in the range (37). Thus the answer toquestion (ii) is also positive.In conclusion, an annular sector can always be closed into a tube withan accompanying axial stretch if its strain-energy function (cid:102) W ( x ) is strictlyconvex and its thickness is small enough, as prescribed by (38). We have shown in the previous section that if the strain-energy function isstrictly convex, then a unique residually-stressed configuration exists for allsuch energy functions and for any given opening angle and axial stretch. Thedimensions of this residually-stressed tube depend on the form of the strain-energy function, the initial geometry of the unstressed sector, and the valueof the axial stretch. For thin tubes, some general statements can be madeabout the form of the residually-stressed configuration, as follows.Expanding f ( y ) in (39) as a Maclaurin series in the thickness parameter ε , substituting into the equation f ( x a ) = 0, and retaining only the leading-order term, yields the equation (cid:99) W (cid:48) ( x a ) = 0 , (46)which, by (27), has the unique solution x a = x m = 1 /k . Thus, thin tubeshave the nice property that the residually-stressed configuration B r corre-sponds to the minimum of the strain energy.Moving on to thicker tubes, we expand equation (36) to the next orderin ε and this yields (cid:99) W (cid:48) ( x a ) + 12 ε (cid:104) (1 − x a ) (cid:99) W (cid:48)(cid:48) ( x a ) − (cid:99) W (cid:48) ( x a ) (cid:105) = 0 . (47)Then we look for x a in the neighborhood of x m , setting x a = x m [1 + εδ + O ( ε )],say. Substitution of this into the equation above yields (cid:20) x m δ + 12 (cid:18) − k (cid:19)(cid:21) (cid:99) W (cid:48)(cid:48) ( x m ) = 0 , (48)13ut (cid:99) W (cid:48)(cid:48) >
0, so that x m δ = − ( k − / (2 k ), and hence x a = x m (cid:2) − ε ( k − / O ( ε ) (cid:3) , (49)showing that a thicker initial sector leads to a smaller inner radius of theresidually-stressed tube (a result independent of the particular form of strain-energy function within the considered class). We now return to the case of tubes with finite thickness. In the precedinganalysis, ε ≡ ( B/A ) − B . It fol-lows easily from simple kinematics that the thickness ratio of the residually-stressed tube in configuration B r , ˆ ε say, is related to ε byˆ ε ≡ (cid:18) ba (cid:19) − εx a > ε. (50)To make further progress, the strain-energy function has to be specified.First, we consider the Mooney-Rivlin model (33). For this material, it isfound that the defining equation (36) for the residually-stressed state is in-dependent of the material constants C and D and does not depend explicitlyon the axial stretch λ z ; it is given byln (cid:18) ε + x a x a (cid:19) − k ln ( ε + 1) + ε (1 − x a ) x a ( ε + x a ) = 0 , (51)or, equivalently, in terms of ˆ ε ,ln (ˆ ε + 1) − k ln (1 + ˆ εx a ) + ˆ ε ˆ ε + 1 (cid:18) − x a x a (cid:19) = 0 . (52)The constants in these equations now need to be specified. Matsumoto andSato (2002) measured the inner and outer radii of five tubular segments ofbovine thoracic aortas, with mean values of a = 11 .
48 mm and b = 17 . ε = 1 . α = 139 ◦ , which yields k = 1 . x a = 0 . then yield ε = 0 . =kλ ( r/R ) z kλ σ z C+Dλ z x a x m x b θθ rr Figure 2: Variation of the non-dimensional radial stress σ rr and hoop stress σ θθ through an arterial wall modeled as a Mooney-Rivlin material. Thetube is in a residually-stressed configuration. The opening angle is taken as α = 139 ◦ .which, incidentally, is indeed smaller than 1 / ( k −
1) = 1 . x is determined as x b = 0 . , and the radial stress reaches its minimum at x = x m = 0 . σ rr and σ θθ , divided by the quantity ( C + Dλ z ) / ( kλ z ), as a function of x . Thesefigures are entirely consistent with those produced by other authors, as in,for example, Ogden (2003), Raghavan et al. (2004) and Olsson et al. (2006).For our second example, we use data from an article by Delfino et al.(1997), who consider measurements made on human carotid arteries. Intheir experimental, theoretical, and numerical investigation of the mechan-ical behavior of arteries, Delfino et al. (1997) model the arterial wall asa homogeneous, hyperelastic, isotropic solid, exhibiting residual stress andstrain-stiffening effects. They use the opening angle method to account forthe former and the Fung model for the latter. They report the following15ypical values: a = 3 . b = 4 . α = 100 ◦ . These give here:ˆ ε = 0 . k = 1 . µ = 44 . c = 8 .
35. It follows that x a can be determined by solving (36), which reads here as0 = (cid:90) x b x a (cid:99) W (cid:48) F ( x )1 − x d x, x b = (ˆ ε + 1) x a ˆ εx a + 1 . (53)Numerically, we find that x a = 0 . A = a (cid:112) kλ z /x a = 4 .
673 mm.This comes within less than 5% of the experimentally measured radius, whichis 4 .
41 mm (Delfino et al., 1997).
In the absence of body forces the equilibrium equation can be written quitegenerally as Div S = , (54)where S is the nominal stress tensor and Div is the divergence operator withrespect to B . For an elastic material with strain-energy function W ( F ) andsubject to the incompressibility constraint det F ≡
1, the nominal stress isgiven by S = ∂W∂ F − q F − , (55)where q is again the Lagrange multiplier associated with the constraint. Interms of the Cauchy stress tensor σ the equilibrium equation has the equiv-alent form div σ = and we note the connection σ = F S . It was thespecialization of the latter form of the equilibrium equation to the cylindri-cal geometry that was used in Section 2.2. In considering the stability of theconfiguration B r it is convenient to work in terms of the incremental form ofthe equilibrium equation (54) and the incremental form of the constitutiveequation (55), which we now summarize.16 .1 Incremental equations Let u be a small displacement from the configuration B r and let ˙ S be theassociated increment in S . The (linearized) incremental form of the consti-tutive equation (55) is written as ˙ S = A ˙ F + q ˙ F − ˙ q F − , (56)where ˙ F and ˙ q are the increments in F and q , respectively, and A is thefourth-order tensor of fixed-reference elastic moduli, defined by A = ∂ W∂ F ∂ F , A αiβj = ∂ W∂F iα ∂F jβ . (57)The incremental form of the equilibrium equation is thenDiv ˙ S = . (58)Let x denote the position vector of a point in B r . Then, by updatingthe reference configuration from B to B r , we may write the incrementalequilibrium equation in Eulerian form asdiv ˙ S = , (59)where ˙ S = A Γ + q Γ − ˙ q I , (60)is the push forward F ˙ S of ˙ S , Γ is the displacement gradient grad u , I is theidentity tensor, and A , the push-forward of A , is the tensor of instantaneousmoduli, whose components are given by A piqj = F pα F qβ A αiβj . For full detailsof these connections, see (Ogden, 1984). The incremental incompressibilitycondition is then tr Γ ≡ div u = 0 . (61)In the configuration B r it is appropriate to work in terms of cylindricalpolar coordinates r, θ, z and their associated basis vectors e r , e θ , e z . Let( u, v, w ) be the components of u . Then the matrix of components of Γ withrespect to this basis is u ,r ( u ,θ − v ) /r u ,z v ,r ( u + v ,θ ) /r v ,z w ,r w ,θ /r w ,z (62)17nd the incremental incompressibility condition becomes u ,r + u + v ,θ r + w ,z = 0 . (63)Referred to principal axes of B the non-zero components of A are (Og-den, 1984): A iijj = λ i λ j W ij , A ijij = ( λ i W i − λ j W j ) λ i / ( λ i − λ j ) , i (cid:54) = j, A ijji = A jiij = A ijij − λ i W i , i (cid:54) = j, (64)(no sums on repeated indexes here), where W j ≡ ∂W/∂λ j and W ij ≡ ∂ W/∂λ i ∂λ j .In the present situation the cylindrical axes are principal axes of B andthe above apply with i and j running over values r , θ , and z . The compo-nents ˙ S ij are then easily obtained from (60) and (62). Explicit expressionsfor the incremental equations of equilibrium for a thick-walled tube in cylin-drical polar coordinates can be found in Haughton and Ogden (1979). Here,however, for purposes of illustration and for simplicity we now specialize theincremental equations so that w = 0 and u and v are independent of z , so thatwe consider prismatic incremental deformations. The equilibrium equation(59) reduces to two component equations, which can be written as( r ˙ S rr ) ,r + ˙ S θr,θ − ˙ S θθ = 0 , ( r ˙ S rθ ) ,r + ˙ S θθ,θ + ˙ S θr = 0 , (65)and from (60) we obtain˙ S rr = ( A rrrr + q ) u ,r + A rrθθ u + v ,θ r − ˙ q, ˙ S rθ = A rθrθ v ,r + ( A rθrθ − σ rr ) u ,θ − vr , ˙ S θr = A θrθr u ,θ − vr + ( A θrθr − σ θθ ) v ,r , ˙ S θθ = ( A θθθθ + q ) u + v ,θ r + A rrθθ u ,r − ˙ q, (66)and the incremental incompressibility equation (63) reduces to u ,r + u + v ,θ r = 0 . (67)18ote that use has been made of the connections A rθθr + q = A rθrθ − σ rr = A θrθr − σ θθ , (68)for which equations (10) and (64) have been utilized.Note that the components of incremental traction on a surface r =constant are ˙ S rr and ˙ S rθ . We now consider solutions of the form { u, v, ˙ q } = { U ( r ) , V ( r ) , Q ( r ) } e i nθ , (69)where U , V , and Q are functions of r only, so that the solutions are singlevalued, and n = 0 , , , , . . . is the circumferential number. It follows thatthe components ˙ S ij have a similar structure and we therefore write˙ S ij = Σ ij ( r )e i nθ , (70)say, where Σ ij , i, j ∈ { r, θ } are functions of r only. Then we obtainΣ rr = ( A rrrr + p ) U (cid:48) + 1 r A rrθθ ( U + i nV ) − Q, Σ rθ = A rθrθ V (cid:48) + 1 r ( A rθrθ − σ rr )(i nU − V ) , Σ θr = 1 r A θrθr (i nU − V ) + ( A θrθr − σ θθ ) V (cid:48) , Σ θθ = 1 r ( A θθθθ + p )( U + i nV ) + A rrθθ U (cid:48) − Q, (71)and the incremental incompressibility condition (67) yields U (cid:48) = − r U − i nr V. (72)From (71) we obtain V (cid:48) = − i nr (1 − σ ) U + 1 r (1 − σ ) V + Σ rθ α , (73)where we have introduced the notations α = A rθrθ , σ = σ rr /α. (74)19y using the above and then eliminating Σ θr , Σ θθ , and Q from the equilibriumequations in favor of U , V , Σ rr , and Σ rθ , the two equilibrium equations yieldexpressions for ( r Σ rr ) (cid:48) and ( r Σ rθ ) (cid:48) in terms of U , V , Σ rr , and Σ rθ . Then, byintroducing the four-component displacement-traction vector η defined by η = [ U, V, i r Σ rr , i r Σ rθ ] T , (75)we may cast the governing equations in the Stroh form (Shuvalov, 2003)dd r η ( r ) = i r G ( r ) η ( r ) , (76)where the matrix G has the Stroh structure, G = i − n − n (1 − σ ) − i(1 − σ ) 0 − /ακ i κ − i − n (1 − σ ) − i κ κ − n i(1 − σ ) (77)where κ = 2 β + 2 α (1 − σ ) + n (cid:2) γ − α (1 − σ ) (cid:3) ,κ = n (cid:2) β + γ + α (1 − σ ) (cid:3) ,κ = γ − α (1 − σ ) + 2 n [ β + α (1 − σ )] , (78)and the notations γ and β are defined by γ = A θrθr , β = A rrrr + A θθθθ − A rrθθ − A rθθr . (79)Note that κ + κ = ( n + 1) κ /n . Note also that this form of the Strohformulation is entirely consistent with previous forms obtained in the bendingof blocks (Destrade et al., 2009, 2010), or the compression of tubes (Goriely etal., 2008), the only differences being in the actual expressions of the moduli.The above coefficients may also be conveniently expressed in terms of (cid:99) W ( x ) through the connections α = 2 x (cid:99) W (cid:48) ( x ) k x − , γ = k x α, β = 2 x (cid:99) W (cid:48)(cid:48) ( x ) + x (cid:99) W (cid:48) ( x ) − α, (80)and σ is given by (19) and (74). 20e may now devise a numerical strategy to find out whether a given opensector remains stable once closed into a tube, as follows. The prescribedgeometrical quantities are the measure of the initial sector thickness ε , themeasure of the opening angle k , and the axial stretch λ z . The prescribedphysical quantity is the strain energy density (cid:99) W ( x ). These quantities arerequired to determine x a and x b from (36). Then we integrate numericallythe incremental equations (76), written asdd x η ( x ) = i2 x (1 − x ) G ( x ) η ( x ) , (81)(using (17)), and we check whether (stability) or not (instability) the innerand outer faces of the tube are left free of incremental traction, i.e. whetheror not η ( x a ) = [ U ( x a ) , V ( x a ) , , T , η ( x b ) = [ U ( x b ) , V ( x b ) , , T . (82)This two-point boundary-value problem may run into numerical stiffnessproblems, especially for thick sectors, and we use the compound matrixmethod (Haughton and Orr, 1997; Destrade et al., 2009) and the impedancematrix method (Destrade et al., 2009) to take care of these. Motivated by simple experimental observations on sectors made of siliconerubbers, see Figure 2, we now focus on tubes made of solids with the Mooney-Rivlin strain-energy density (33). We find that α = ( Cλ − z + Dλ z ) 1 kx , γ = ( Cλ − z + Dλ z ) kx, β = α + γ, (83)and σ = x (cid:20) (1 − k ) ln (cid:18) x − x b − (cid:19) − ln (cid:18) xx b (cid:19) + 1 x − x b (cid:21) . (84)In the latter equation, x b is determined by solving (51) for x a , a calculationwhich is done independently of C , D , and λ z .It is now clear upon inspection of the incremental equations (81) thatonce we normalize the third and fourth components of η with respect to thequantity ( Cλ − z + Dλ z ), we end up with a completely non-dimensional systemof differential equations which is independent of the material parameters C A = 3 . B = 7 . α = π/ α = 2 π/ D and of the axial stretch λ z . Only three parameters are then requiredto determine the stability of closed-up Mooney-Rivlin tubes, namely thefollowing geometrical quantities: B/A , the ratio of the initial outer to innerradii; α , the opening angle; and n , the mode number.Figure 4 shows the variations of the critical opening angle α cr (in degrees)with B/A , for different values of n . The lowest curve represents the bifurca-tion curve , below which all tubes are stable. In our numerical investigation,we found that this corresponds to the n = 9 curve, although the n = 8 and n = 10 curves are almost indistinguishable from it. It follows that when22 closed-up Mooney-Rivlin tube buckles, it presents at least 8 wrinkles onits inner curved face. When the initial sector is thin, B/A →
0, it can bebent into a tube without buckling, so that α cr → ◦ in that limit, for allmode numbers. This is consistent with the observation that thin rectangularplates with a neo-Hookean strain energy can be bent into full circles withoutencountering bifurcation (Haughton, 1999).The radii ratio of the silicone sectors in Figure 3 is B/A = 2 .
06. We seethat an opening angle of α = 120 ◦ is below the bifurcation curve, so thatthe corresponding sector can be closed into a tube in a stable manner; seePoint (a) in Figure 4. For an opening angle of α = 240 ◦ , we are above thebifurcation curve (see Point (b) in Figure 4), and wrinkles are thus predicted.Finally, the data of Matsumoto and Sato (2002), fitted to the Mooney-Rivlinstrain energy, give B/A = √ x a ˆ ε = 1 .
29; see Section 3.2. In that case,their measured opening angle of α = 139 ◦ is clearly below the bifurcationcurve, see Point (c) in Figure 4, and the artery can be “closed” in a stablemanner. We have proved that the opening angle method gives a unique solution tothe equations of equilibrium provided the strain energy density of the solidis a convex function of the geometrical quantity x = 2 πλ z ( r /R ) / (2 π − α ),where λ z is the axial pre-stretch, α is the opening angle, and r , R are theinitial and current radial coordinates, respectively. The convexity is indeedmet by many standard strain energy densities including the Mooney-Rivlin,Fung, and Gent models. We have also provided an analysis of incrementalprismatic deformations superimposed on the finite deformation in order toassess the stability of the intact tube formed by the closing of the sector.The results are valid for homogeneous isotropic solids and are thus highlyrelevant to the mechanical modeling of the closing of annular sectors intotubes made of rubber-like materials. This correlation was confirmed by ourtoy experiments with sectors made of soft silicone. With respect to biome-chanics, where the opening angle method is a tool widely used for the esti-mation of residual stresses in arteries, we note that there are many articlesin the literature where arterial walls are modeled as being homogeneous andisotropic, as, for example, in Halperin et al. (1991) or Delfino et al. (1997).However, it is now well established and accepted that arteries are inhomo-23 /A cr α n= n= n= a )( c ) ( b ) Figure 4: Instability of a closed tube with Mooney-Rivlin strain energy:critical opening angle α cr versus initial radii ratio B/A , for different modenumbers n . The lowest curve –the bifurcation curve– corresponds to modes n = 8, 9, 10 (almost indistinguishable one from another). The curves ofall other modes are situated above this ( n = 2, 3 are shown). The annularsilicone sector on the left in Figure 3 and the opened artery of Section 3.2correspond to stable geometries (Points (a) and (c) ). The annular siliconesector on the right in Figure 3 gives a point above the bifurcation curve,Point (b) .geneous and anisotropic solids, due, in particular, to a layered structure andthe presence of collagen fiber bundles. The impact of such inhomogeneityand anisotropy on our results remains to be assessed.24 cknowledgements This work is supported by a Senior Marie Curie Fellowship awarded by theSeventh Framework Programme of the European Commission to the firstauthor, and by an E.T.S. Walton Award given to the third author by ScienceFoundation Ireland. This material is based upon works supported by theScience Foundation Ireland under Grant No. SFI 08/W.1/B2580.We are most grateful to Stephen Kiernan for taking the pictures of Fig-ure 3.
References
Chuong, C.J., Fung, Y.C., 1986. On residual stress in arteries, J. Biomech.Eng. 108, 189-192.Delfino, A., Stergiopulos, N., Moore, J.E., Meister, J.-J., 1997. Residualstrain effects on the stress field in a thick wall finite element model of thehuman carotid bifurcation, J. Biomech. 30, 777-786.Destrade, M, Gilchrist, M.D., Motherway, J.A., Murphy, J.G., 2010. Bimod-ular rubber buckles early in bending, Mech. Materials 42, 469–476.Destrade, M., N`ı Annaidh, A., Coman, C.D., 2009. Bending instabilities ofsoft biological tissues, Int. J. Solids Struct. 46, 4322–4330.Fung, Y.C., 1993.
Biomechanics: Mechanical Properties of Living Tissue ,2nd Edition. New York: Springer.Goriely, A., Vandiver, R., Destrade, M., 2008. Nonlinear Euler buckling,Proc. Roy. Soc. A, 464, 3003-3019.Greenwald, S.E., Moore, J.E., Rachev, A., Kane, T.P.C., Meister, J.-J., 1997.Experimental determination of the distribution of residual strains in theartery wall, J. Biomech. Eng. 119, 438-444.Humphrey, J.D., Halperin, H.R., Yin, F.C.P., 1991. Small indentation super-imposed on a finite equibiaxial stretch: Implications for cardiac mechanics,J. Appl. Mech. 58, 1108-1111. 25aughton, D.M., 1999. Flexure and compression of incompressible elasticplates, Int. J. Eng. Sc. 37, 1693–1708.Haughton, D.M., Ogden, R.W., 1979. Bifurcation of inflated circular cylin-ders of elastic material under axial loading II. Exact theory for thick-walledtubes, J. Mech. Phys. Solids 27, 489-512.Haughton, D.M., Orr, A., 1997. On the eversion of compressible elastic cylin-ders, Int. J. Solids Struct. 34, 1893-1914.Hoger, A., 1985. On the residual stress possible in an elastic body withmaterial symmetry, Arch. Rational Mech. Anal. 88, 271-290.Holzapfel, G.A., Gasser, T.C., Ogden, R.W., 2000. A new constitutive frame-work for arterial wall mechanics and a comparative study of material mod-els, J. Elasticity 61, 1-48.Holzapfel, G.A., Ogden, R.W., 2010. Modelling the layer-specific 3D residualstresses in arteries, with an application to the human aorta, J. R. Soc.Interface 7, 787-799.Humphrey, J.D., 2002.
Cardiovascular Solid Mechanics. Cells, Tissues, andOrgans , New York: Springer-Verlag.Matsumoto, T., Sato, M., 2002. Analysis of stress and strain distributionin the artery wall consisted of layers with different elastic modulus andopening angle, JSME Int. J. C45, 906-912.Ogden, R.W. 1984.
Non-Linear Elastic Deformations , Ellis Horwood, Chich-ester. Reprinted by Dover, New York, 1997.Ogden, R.W., 2003. Nonlinear elasticity, anisotropy and residual stresses insoft tissue (lecture notes, CISM Course on the
Biomechanics of Soft Tissuein Cardiovasular Systems ), pp. 65–108. CISM Courses and Lectures Series , Springer, Wien.Ogden, R.W., Schulze-Bauer, C.A.J., 2000. Phenomenological and structuralaspects of the mechanical response of arteries, in J. Casey and G. Bao (eds),
Mechanics in Biology , AMD-Vol. 242/BED-Vol. 46, The American Societyof Mechanical Engineers, New York, pp. 125-140.26lsson, T., St˚alhand, J., Klarbring, A., 2006. Modeling initial strain dis-tribution in soft tissues with application to arteries, Biomech. Model.Mechanobiol. 5, 27-38.Rachev, A., Greenwald, S.E., 2003. Residual strains in conduit arteries, J.Biomech. 36, 661-670.Rachev, A., Hayashi, K., 1999. Theoretical study of the effects of vascularsmooth muscle contraction on strain and stress distributions in arteries,Ann. Biomed. Eng. 27, 459-468.Raghavan, M.L., Trivedi, S., Nagaraj, A., McPherson, D.D., Chandran, K.B.,2004. Three-dimensional finite element analysis of residual stress in arter-ies, Ann. Biomed. Eng. 32, 257-263.Shuvalov, A.L., 2003. A sextic formalism for three-dimensional elastody-namics of cylindrically anisotropic radially inhomogeneous materials, Proc.Roy. Soc. Lond. A459, 1611–1639.Vaishnav, R.N., Vossoughi, J., 1983. Estimation of residual strains in aorticsegments, in C.W. Hall (ed.),
Biomedical Engineering II: Recent Develop-ments , Pergamon Press, New York, pp. 330-333.Vossoughi, J., Hedjazi, Z., Boriss, F.S., 1993. Intimal residual stress andstrain in large arteries, in