## Abstract

Cylindrical forms are among one of Nature’s fundamental building blocks. They serve many different purposes, from sustaining body weight to carrying flows. Their mechanical properties are generated through the often complex arrangements of the walls. In particular, in many structures that have elastic responses, such as stems and arteries, the walls are in a state of tension generated by differential growth. Here, the effect of differential growth and residual stress on the overall mechanical response of the cylindrical structure is studied within the framework of morpho-elasticity.

## 1. Introduction

The mechanical description of growth is potentially of great interest in many different fields, as growth plays a fundamental role in both normal development processes, including morphogenesis and homeostasis, and many pathological disorders. Growth, the increase of mass in a system, may occur through different mechanisms. For example, the growth of hard tissue, such as teeth, seashells and horns, is appositional, i.e. growth occurs by accretion and deposition on an existing surface (Skalak & Hoger 1997). Volumetric growth, however, takes place in the bulk of the material and is the generic mechanism for growth in many soft tissues (Taber 1995; Humphrey 2002).

Differential growth within biological tissues may produce residual stress, i.e. stress remains in a body when all external loads are removed. This stress typically arises from an incompatibility of growth. Skalak *et al.* (1996) put forth the key ideas to develop a theory of kinematic growth. They suggested that if the growth of cells compromises continuity of the body, an elastic deformation may be required in order to maintain integrity and compatibility. These ideas were further refined by Rodriguez *et al.* (1994), who formulated a general three-dimensional theory of volumetric growth for soft tissues undergoing large deformations. The basic idea is to decompose the total deformation into a sequence of two mappings: a local addition of material represented by a growth tensor and an elastic response represented by an elastic tensor. A similar decomposition found in elasto-plasticity was introduced by Kroner (1960) and Lee (1969), who separated the finite elastic and plastic deformations. Rachev *et al.* (1998) and Taber (1998*a*,*b*), among others, have applied the theory of kinematic growth presented by Rodriguez *et al*. to study problems in vascular mechanics.

Residual stress, arising from differential growth, is observed in many biological systems and may be associated with important mechanical effects. By making cuts in a tissue, a change in shape may occur, which reveals the presence of residual stress in the body. The existence of residual stress in the trachea (Han & Fung 1991), ventricular wall (Taber *et al.* 1993) and arteries (Vaishnav & Vossoughi 1983, 1987; Liu & Fung 1988) has been demonstrated by performing radial cuts on unloaded rings. Thin rings of each of these tissues will open up into sectors, indicating the presence of internal stresses in the intact structures. The residual stresses in the trachea are known to affect the acoustic characteristics of the trachea and promote elastic instability or collapse of the trachea when it is subjected to compression, as in sneezing or coughing (Han & Fung 1991). Residual stress may also play a role in cardiac morphogenesis by affecting ventricular wall stress during cardiac development (Taber *et al.* 1993). Similarly, residual stresses in arteries are known to play a fundamental role in regulating the peak stress and homogenizing the stress field across the arterial wall (Chuong & Fung 1986). While the role of transverse differential growth in cylindrical structures is well appreciated, the overall effect of axial residual stress on the mechanics of a structure is not completely understood.

Residual stress is also observed in many plants, including wheat roots, fennel leaves, rhubarb stalks, sweetgum trees and hypocotyls of cucumber, sunflower, cantaloupe and squash (Burstrom 1971; Kutschera 1989; Brown *et al.* 1995*a*,*b*; Hejnowicz & Sievers 1995, 1996; Hejnowicz 1997; Peters & Tomos 2000). If these cylindrical structures are cut along the longitudinal axes, they will tend to bend outwards, as part of the elastic stress is relieved when the inner tissues elongate and the outer tissues shorten by curving. The existence of tissue tension in plants has been known since as early as the nineteenth century when Hofmeister (1859) observed the outer and inner tissues of *Vitis vinifera* (common grape vine) elastically contract and extend upon separation. This simple experiment can be easily reproduced on a stalk of rhubarb (Atkinson 1900, p. 47). A stalk of rhubarb is made up of an outer layer, consisting of the epidermal tissue and the collenchyma layers, and an inner layer consisting chiefly of parenchyma (figure 1). When a strip of the stalk’s outer layer is peeled, the strip shortens in length (1.5–2.5% in figure 1). When the outer layer is completely removed, the inner part (the pith) extends in length by a noticeable amount (6% in figure 1). This simple experiment shows that the outer wall is in a state of axial tension, while the pith is in a state of axial compression. The inner tissue grows at a greater rate than the outer tissue, which creates tissue tension in the intact plant. Tissue tension in plants played an important role in the discovery of plant growth hormones such as auxins. The demonstration of the existence of hormones rely on the differential effects that a hormone would have on different parts of a growing tissues, resulting in a different shape when the tissue is split. This is the basis of the famous ‘pea test’(Bonner 1933; Thimann & Schneider 1938; Schneider 1942; Masuda 1990; Peters & Tomos 1996).

The possible mechanical role of the residual stresses and the combination of tissues can be appreciated by realizing that the peeled rhubarb has lost most of its rigidity, so much so that it now buckles under its own weight. The possible role of tissue tension in plant mechanics was observed by Sachs: ‘We have here the case of an elastic stiff body consisting of two parts, each in a high degree flexible and by no means stiff; only in their natural connection do the epidermal tissue and internal tissues together form an elastic rigid body’ (Sachs 1882). While it is evident that many plant stems develop axial tissue tension, the overall effect on the mechanics is not yet understood. Similarly, many other cylindrical structures, such as tree trunks, roots and arteries (Holzapfel *et al.* 2000), also develop axial tissue tension and the mechanical role of these residual stresses is not known.

Our goal is to study the development of residual stress in differential growth of biological cylindrical structures and to elucidate its possible mechanical role in modifying material properties. We start by giving a general formulation of growth within a three-dimensional nonlinear elastic body and apply to a cylindrical geometry, which is relevant in many physiological and biological systems (such as the growth of arteries and stems). We model the biological structure as a cylindrical shell composed of two material layers with different growth and elastic properties. Because soft tissues undergo large elastic deformations in response to physiological loading, their physical properties are best described within the theory of finite elasticity.

The specific goal of this paper is to determine the overall rigidity of a grown cylinder, and this is done by using the Euler buckling formula to obtain an effective Young’s modulus. However, in order to obtain an estimate of the cylinder rigidity through an effective Young’s modulus, one needs to obtain the critical buckling load by performing a bifurcation analysis. This analysis is done through a perturbation expansion in which an incremental deformation is superimposed on some finite deformation. If a solution is found to the incremental equation with the appropriate boundary conditions, then the possibility of instability exists and the critical load is determined. This method allows us to understand how residual stresses affect the overall rigidity of the system. We apply this method to cylinders made of isotropic materials, as well as anisotropic materials.

## 2. Kinematics

Consider a continuous body whose *reference configuration* is defined by . Let ** X** denote the position vectors in . Now suppose the body is deformed to a new configuration, . We refer to as the

*current configuration*, where the position of a material point

**is now defined by**

*x***=**

*x***(**

*χ***,**

*X**t*). The deformation gradient, F(

**,**

*X**t*)=grad

**, relates a material segment in the reference configuration to the same segment in the current configuration. The main assumption is that the total deformation can be decomposed into a growth tensor G(**

*χ***,**

*X**t*) and an elastic tensor A(

**,**

*X**t*). This key idea of decomposing the total deformation was introduced by Rodriguez

*et al.*(1994) and takes the form 2.1 That is, the growth deformation tensor may not result in a continuous change from point to point and may not be compatible. However, if we require continuity as the body grows, then an elastic deformation is introduced to maintain compatibility. As shown in figure 2, the growth tensor G(

**,**

*X**t*) maps to the virtual configuration , which is locally stress free. In order to maintain continuity of the body, the elastic deformation is introduced, which in turn causes residual stress in the grown body .

## 3. Equations of motion and Cauchy stress

The forces distributed on a body include a contact-force density *t*_{n} and a body-force density ** b**. In accordance with Euler’s laws of motion, the balance of linear momentum is given by
3.1
where

*ρ*is the mass density of the material of body and is the time derivative of the velocity vector such that . According to Cauchy’s theorem, the contact-force density depends linearly on the unit normal

**, given by**

*n*

*t*_{n}=T

**, where T is referred to as the Cauchy stress tensor. Using Cauchy’s theorem and applying the divergence theorem to equation (3.1) yields the equilibrium equation 3.2 Balancing the moments of the forces acting on the body reveals that the stress tensor is symmetric, i.e. T**

*n*^{T}=T. If the body is at rest and the body forces are assumed to be negligible, the equation for mechanical equilibrium becomes 3.3 where the divergence is taken with respect to

**in the current configuration. The solutions to the equilibrium equations must satisfy the conditions imposed on the boundary, which can be in the form of dead loading, rigid loading or mixed loading. Dead loading prescribes the normal components of the stresses at the boundary, rigid loading imposes a fixed deformation at the boundary and mixed loading imposes fixed deformations on some part of the body and stresses on the remaining boundaries.**

*x*We assume that the body is hyperelastic. That is, the material can be described by a strain-energy function *W*=*W*(A), and the Cauchy stress is related to the elastic deformation by
3.4
where *p* is a Lagrange multiplier associated with the internal constraint of incompressibility and *J*=det(A) represents the change of volume owing to the elastic deformation. For an incompressible material, *J*=1 and *p* is the hydrostatic pressure. If the material is compressible, then *p*=0.

## 4. Growth formulation for a two-layered model

We consider a two-layered cylinder, with an inner cylinder of initial radii *A* and *B* and an outer layer of initial radii *B* and *C* (shown in figure 3). Initially, we consider a simple finite deformation in which the cylinder is allowed to grow and deform while remaining cylindrical. We consider the finite growth and elastic deformation up to a critical value of the parameters where the linearized incremental problem around the finite deformation supports solutions that do not respect the cylindrical symmetry. That is, the finite deformation of each cylindrical shell is given by
4.1
where the subscripts ‘i’ and ‘o’ correspond to the inner and outer layer, respectively, and
4.2

The geometric deformation in cylindrical coordinates is
4.3
4.4
where the gradient is taken in the reference configuration, the prime denotes differentiation with respect to *R*_{i,o} and *λ* is the axial stretch of the cylinder such that *l*=*λ**L*. The growth tensor is given by
4.5
and the elastic tensor is given by
4.6
where we have used the incompressibility condition det(A_{i,o})=1 and *α*_{i,o}=*r*_{i,o}/(*R*_{i,o}*γ*_{θi,o}). The only non-vanishing equation in (3.3) is
4.7
The strain-energy function *W*_{i,o} relates the elastic deformation to the Cauchy stress tensor T=diag(*t*_{r},*t*_{θ},*t*_{z}) by
4.8
which can be used in equation (4.7) to obtain
4.9
where . Assuming the boundary conditions of zero normal traction on the outer boundary and an internal pressure, *P*, at the inner boundary (*t*_{r}(*a*)=−*P* and *t*_{r}(*c*)=0), we integrate equation (4.9) to obtain an equation for the radial stress
4.10
The axial stress is given by the *z**z* component of equation (4.8), i.e.
4.11
where *W*_{1} and *W*_{3} are the derivatives of *W* with respect to its first and third variables. At the surface of stress discontinuity at *r*=*b*, the associated tractions need to be equal but opposite, which requires the radial stress to be continuous across *r*, i.e.
4.12
The two free parameters *a* and *λ* are obtained by solving simultaneously the boundary condition (4.12) for a given resultant axial load *N*,
4.13
and the deformation is completely determined. The reduced axial force, *F*, is given by *F*=*N*−*π**a*^{2}*P*, which gives the force measured on the cylinder.

## 5. Rigidity of the cylinder

The Euler buckling formula gives the critical axial load for the buckling of a column, i.e.
5.1
where *F*_{crit} is the critical force, *E* the modulus of elasticity, *I* the area moment of inertia and *l*_{0} the unloaded length of the column. For a cylinder of inner radius *a*_{0} and outer radius *b*_{0}, the area moment of inertia is , and therefore the critical load becomes
5.2

Since the cylinder is in large deformation and supports residual stress, the traditional notion of the Young modulus does not apply directly, and an explicit computation of the rigidity is not possible. However, we can perform the following thought experiment to obtain an effective Young’s modulus. Consider the grown cylinder, measure all its geometric parameters (inner and outer radii *a*_{0} and *b*_{0} and length *l*_{0} being the geometric parameters for zero load) and subject it to a normal load *N* until it buckles at *F*_{crit}. The critical buckling pressure can be used to compute the effective Young modulus, defined as the Young modulus of an equivalent homogeneous elastic cylinder with the same geometry, but with no residual stress and buckling at the same pressure, i.e.
5.3
where is the inverse of the slenderness ratio. The Euler buckling formula is valid for cylinders that are sufficiently slender. Therefore, in order to obtain an estimate of the cylinder rigidity through an effective Young modulus, one needs to obtain the critical buckling load, i.e. to compute the stability of the grown cylinder under applied loads. The effective Young modulus, *E*_{eff}, can then be obtained from the slope of the graph *P*_{crit}=*π*^{2}σ^{2}*E*_{eff}/4.

## 6. Incremental deformation

To perform the stability of the grown cylinder under applied loads, we consider an incremental deformation superimposed on the finite deformation. A deformation ** χ** is introduced as
6.1
where

*χ*^{(0)}is the finite deformation and

*χ*^{(1)}is the incremental deformation that is now a general mapping with no prescribed symmetry. The small parameter

*ε*is introduced to characterize the size of the perturbation. It follows that the total deformation tensor and the elastic tensor are given by where

**A**

^{(0)}=

**F**

^{(0)}·

**G**

^{−1}and

**A**

^{(1)}=

**F**

^{(1)}. The Cauchy stress tensor is expanded around the finite state in

*ε*given by 6.2 and the constitutive relationship to zeroth and first orders reads 6.3 and 6.4 where

*p*=

*p*

^{(0)}+

*p*

^{(1)}

*ε*, the derivatives of

*W*are with respect to

**A**evaluated at

**A**

^{(0)}and is the fourth-order tensor of instantaneous elastic moduli defined by 6.5 In Cartesian components, the instantaneous elastic moduli are 6.6 where is the elastic moduli defined by . The non-zero components of are given by (Ogden 1984) 6.7 The equilibrium equation to first order is 6.8 where div is the divergence in the current configuration. Substituting equation (6.4) into equation (6.8) and using (

**T**

^{(i)})

^{T}=(

**T**

^{(i)}), div(

**T**

^{(0)})=0 and div(

**F**

^{(1)})=0, the incremental equations for an incompressible material are given by 6.9 Explicitly, these three equations are given by 6.10 6.11 6.12 where .

Assuming no extra loading is imposed on the surface of the cylinder, the incremental boundary conditions are given by (**T**^{(1)})^{T}·** N**=0. For given boundary conditions, bifurcation occurs if there exists a non-trivial solution to the system of equations given in (6.10)–(6.12).

## 7. Bifurcation analysis

The general form of the incremental deformation is given by the vector
7.1
The incremental deformation tensor **F**^{(1)} is then given by
7.2
where (*r*,*θ*,*z*) subscripts denote partial derivatives. The incompressibility constraint is
7.3
Substituting equation (7.2) into the incremental equations given in equation (6.9) results in
7.4
7.5
and
7.6
where the prime denotes differentiation with respect to *r*. The set of partial differential equations for the variable *u*,*v*,*w*,*p*^{(1)} can be further simplified to a boundary-value problem for a set of differential equations in the variable *r* by Fourier expanding the dependence in *θ* and *z*, i.e.
7.7
In order that the incremental displacements be single valued, we take the mode number *m* to be an integer *m*≥0. As the incremental end displacement *w* must vanish on the ends, we have
7.8
where *n* is set to unity so that η is inversely proportional to the length of the cylinder. The expressions in equation (7.7) are substituted in equations (7.4)–(7.6), and *h*(*r*) can be eliminated by means of the incompressibility constraint (7.3). The resulting incremental equilibrium equations are
7.9
7.10
and
7.11
The corresponding incremental boundary condition on the inner and outer wall *r*=*a*,*c* is **T**^{(1)}·** n**=

*P*

**F**

^{(1)}η

**−**

*n**P*

^{(1)}

**F**

^{(0)}

**, which takes the form 7.12 where the incremental internal pressure,**

*n**P*

^{(1)}, is taken to be

*P*

^{(1)}=0. The three boundary-condition equations on

*r*=

*a*,

*c*are given by 7.13 7.14 7.15 One must now consider the boundary conditions at the interface of the two layers. It is assumed that there is no penetration at the interface and, therefore,

*f*

_{i}=

*f*

_{o},

*g*

_{i}=

*g*

_{o}and

*h*

_{i}=

*h*

_{o}. Recall that

*h*can be eliminated by using the incompressible condition in equation (7.3) and, therefore, the boundary condition

*h*

_{i}=

*h*

_{o}is equivalent to requiring

*f*′

_{i}=

*f*′

_{o}. Lastly, the radial stress at the interface is required to be equal. The other boundary conditions are 7.16 7.17 7.18 7.19 7.20 7.21 The growth dependence comes about through the coefficients of the system of equations above. The instantaneous elastic moduli depend on

*α*and its derivatives through the relation 7.22

### (a) Numerical integration

The coupled equations for *f*(*r*), *g*(*r*) and *k*(*r*) can be written in the form
7.23
where *y*=(*f*,*f*′,*f*′′,*g*,*g*′,*k*). The boundary conditions (7.13)–(7.21) are given by a set of 12 linear functions *c*_{1,…,12}(*y*(*r*);*r*),
7.24
and
7.25
The determinant method (Haughton & Ogden 1979; Haughton & Orr 1995; Ben Amar & Goriely 2005) is used in order to find values of the parameters for which a solution exists. Consider three copies of the system (7.23), i.e.
7.26
with linearly independent conditions , and that satisfy the conditions at *r*=*a*, such that . Integrate this system up to *r*=*b*, then obtain the initial conditions , and from the end conditions , and using the boundary conditions at the interface in equation (7.25). Integrate this system up to *r*=*c*, then evaluate the 3×3 determinant at *r*=*c* given by
7.27
Bifurcation occurs when there is a linear combination of , and that satisfies the system (7.23) with boundary conditions (7.24) and (7.25). For some initial length of the cylinder, we iterate on *λ*_{z} until the determinant vanishes.

## 8. Results

For each mode *m*, the solution of the system of ordinary differential equations in (7.9)–(7.11) with boundary conditions (7.13)–(7.21) is possible only for a particular combination of parameters *m*, *a*, *b* and *c*. When the determinant vanishes in equation (7.27), the critical stretch, *λ*_{zcrit}, is known and the corresponding critical load *F*_{crit} is now determined. Therefore, for the buckling mode *m*=1, the corresponding overall rigidity of a grown cylinder can be computed and compared with a cylinder without growth.

To test the method, the bifurcation curve for the buckling mode *m*=1 is computed for a neo-Hookean two-layer cylindrical model. The strain-energy function is given by *W*=*μ*_{i,o}/2(*I*_{1}−3) where *μ*_{i} and *μ*_{o} are the elastic moduli for the inner and outer layers. The elastic moduli, *μ*_{i} and *μ*_{o}, are related to the Young modulus *μ*_{i,o}=*E*_{i,o}/3 for small deformation. Although modes *m*=0,2, which correspond to axisymmetric barrelling solutions, exist for shorter tubes, the focus of our study here will be on the buckling mode *m*=1. When both layers have the same modulus *μ*_{i}=*μ*_{o}, the buckling problem is equivalent to the classical Euler buckling of a cylindrical tube (Goriely *et al.* 2008). In figure 4, the solution to the bifurcation analysis is compared with the Euler buckling formula by using the relationship *E*_{i,o}=3*μ*_{i,o} for small incremental deformation. The Euler buckling formula is shown to be an excellent approximation for small σ. In the case where Euler’s formula cannot be applied (*μ*_{i}≠*μ*_{o}), the graph of *P*_{crit} as a function of σ^{2} becomes linear for σ small enough. Equation (5.3) can then be used to find *E*_{eff} and, therefore, this bifurcation analysis provides a method to define an effective Young modulus. We now apply this method to three cases: an isotropic one-layer cylinder with radial-dependent growth, an isotropic two-layer model with constant growth and an anisotropic two-layer model with constant growth.

### (a) Isotropic, one-layer model with radial dependent growth

We first consider a neo-Hookean one-layer cylindrical model where the strain-energy function is given by . The growth tensor G is a function of the radial position *r* in the current configuration. Assume a simple, linear growth law in the axial direction, *γ*_{z}=1+κ(*r*−*a*), where κ is the growth parameter, and assume no growth occurs in the radial and circumferential directions. It is interesting to consider the residual stresses created in the grown material in the absence of loads, *N*=0, for various values of the growth parameter. As the axial growth is increased, we note the large zone of compressive force along the axial direction in figure 5*a*. We can now determine the overall rigidity of the grown cylinder. We assume that the cylindrical shell has grown axially by a given amount, then we test the stability of the resulting residually stressed cylindrical shell. In figure 5*b*, we observe that the effect of axial growth is to reduce the stability of the structure. This is due to the fact that growth creates a large zone of compression in the cylinder shown in figure 5*a*. The cylinder is precompressed due to growth and, therefore, the cylinder buckles under a smaller load.

### (b) Isotropic, two-layer model with constant growth

The one-layer cylinder is not sufficient to correctly model many biological cylindrical structures. For example, the plant stem consists of two material layers—the inner layer (pith) being a rather compliant material and the thin outer layer (epidermis) a more rigid layer. In this section, we model the biological structure as a cylindrical shell composed of two material layers with different growth and elastic properties. We assume that the material is hyperelastic, incompressible, homogeneous, isotropic and subject to growth, and each material in the two-layered cylindrical tube is characterized by a strain-energy function *W*_{i,o}=*W*_{i,o}(A_{i,o}). Here, following existing models for arteries (Holzapfel *et al.* 2000) and the data available for stem properties (Hejnowicz & Sievers 1996), we adopt a model where the inner cylinder is a neo-Hookean material and the outer layer is a Fung material with strain-stiffening properties
8.1
where *μ*_{i,o} are the elastic material constants and ν controls the strain-stiffening property of the outer layer (with the neo-Hookean limit as ).

We estimate the material constants *μ*_{i} and *μ*_{o} and the strain-stiffening parameter ν from a study done by Hejnowicz & Sievers (1996), in which they performed uniaxial tests on the inner and outer tissues of the hypocotyl of the sunflower (*Helianthus annuus*). Strain stiffening in the outer layer is evident from the data points (fig. 2, Hejnowicz & Sievers 1996), and therefore the data points were fitted to the axial stress *t*_{z}=*α*_{z}∂*W*_{o}/∂*α*_{z}−*p*, where *W*_{o} is the Fung strain-energy function given in equation (8.1). In a uniaxial test on a rectangular strip of tissue, the elastic tensor is given by A=diag(*α*,*α*,1/*α*^{2}). The vertical surfaces are stress free and, therefore, the boundary condition *t*_{x}=*t*_{y}=0 implies *p*=*α*∂_{x}*W*_{o}. Substituting the strain-energy function *W*_{o} into the stress–strain relationship, the axial stress is given by
8.2
Nonlinear regression can be performed to fit equation (8.2) to the experimental data and the parameters *μ*_{o} and ν can be determined. The results produced values *μ*_{o}=6.7 and ν=195 with an *R*^{2} value of 0.992. The large value of ν is needed to have this type of exponential behaviour for such small strains. Figure 6 compares the approximate experimental data points and the corresponding stress–strain fit. The Fung strain-energy function is shown to be an excellent fit to the data.

A uniaxial test was also performed on the inner tissue (fig. 6, Hejnowicz & Sievers 1996). Unlike that for the outer tissue, the stress–strain curve for the inner tissue was nearly linear within the range of forces applied (0–0.4 N), and therefore we assume that the tissue is composed of a neo-Hookean material. The axial stress for a strip of neo-Hookean material is given by
8.3
A nonlinear regression can be performed to fit equation (8.3) to the experimental data while solving for the parameter *μ*_{i}. In doing so, the value of *μ*_{i} is found to be 1.25 with an *R*^{2} value of 0.991. Figure 6 compares the approximate experimental data points (Hejnowicz & Sievers 1996) and the corresponding stress–strain fit. It should be noted that the outer tissue is a complex tissue containing both the epidermis and one or two layers of cortical parenchyma.

Now that the form of the strain-energy functions are known, it is particularly interesting to consider the residual stresses created in the grown material in the absence of loads, *N*=0. In figure 7, we show the axial stress profile as we vary the ratio of the inner layer stiffness to the outer layer stiffness and the ratio of the inner layer width to the outer layer width (*w*_{i}/*w*_{o}), where *w*_{i}=(*B*−*A*) and *w*_{o}=(*C*−*B*). As the thickness of the inner layer is decreased, we note that large stress gradients are created in the outer layer.

We can now perform a bifurcation analysis and compute the overall rigidity as a function of the parameters. In order to take into account the geometry of the structure, we consider the flexural rigidity, *E*_{eff}*I*. We would like to understand the specific effect of axial growth and therefore set *γ*_{ri,o}=*γ*_{θi,o}=1. Without loss of generality, we take *γ*_{zo}=*μ*_{i}=1. We study the effect of the ratio of the thickness of the inner tissue (*w*_{i}) to the thickness of the outer wall (*w*_{o}). The radii *B* and *C* are set to 0.9 and 1, respectively, and the inner radius *A* is varied. In the absence of growth, the stiffness of the outer layer *μ*_{o} provides a substantial improvement on the overall flexural rigidity of the structure, which is essentially linear in the ratio (figure 8). However, the thickness of the inner layer has very little effect.

We can now fix the values of *μ*_{o} and consider the effect of growth. In figure 9, we observe that, as the thickness of the inner layer increases, even modest differential axial growth has a large effect on the rigidity of the structure. Increasing the thickness of the inner layer creates a large tensile stress in the outer layer (figure 7) and, therefore, a large force must be imposed to incur the stresses needed on the concave side of the cylinder to produce bending. A similar effect due to the increase in the strain-stiffening parameter ν was noticed in Vandiver & Goriely (2008), where an estimate based on a simple physical argument was also provided.

### (c) Anisotropic two-layer model with constant growth

Thus far, we have neglected an important effect necessary to obtain a precise description of the material properties of biological structures, namely the anisotropic response. In particular, tissues such as blood vessels and plant stem walls exhibit anisotropic behaviour due to the fibre orientations that tend to have preferred directions (Holzapfel *et al.* 2000; Baskin 2005). In this section, we study the effect of anisotropy, internal pressure and residual stress on the overall mechanical stability of a two-layer cylinder.

Here, rather than looking at a general situation, we use the multi-layer model for arteries introduced by Holzapfel *et al.* (2000). The arterial wall consists of three layers, each containing specific histological features: intima (the innermost layer), media and adventitia (the outer layer). The intima is known to have negligible (solid) mechanical contributions (Holzapfel *et al.* 2004) and, therefore, the arterial wall is approximated as a two-layer structure containing the media and adventitia. Each layer is reinforced with collagen fibres with preferred directions that render the material anisotropic.

The contribution of collagen fibres in the arterial wall is considered to be negligible at low pressures and therefore the mechanical response is assumed to be isotropic. However, at higher pressures, the collagen fibres are stretched and the resulting mechanical response is anisotropic. Therefore, for each layer, Holzapfel *et al.* (2000) suggest a strain energy of the form
8.4
where the invariants *I*_{1}, *I*_{4} and *I*_{6} of the right Cauchy–Green tensor C=F^{T}·F are given by
8.5
The invariants *I*_{4} and *I*_{6} are the squares of the stretches in the directions of the two families of collagen fibres. The unit vectors ** M** and

**′ describe the orientations of these fibres, and these vectors are defined for each layer in the reference configuration. In a cylindrical coordinate system, the components of the direction vectors are 8.6 where the subscripts ‘i’ and ‘o’ refer to the inner layer (media) and the outer layer (adventitia), respectively, and**

*M**ϕ*

_{j}are the angles between the collagen fibres (figure 10). The isotropic and anisotropic responses in equation (8.4) are given by 8.7 and 8.8 where the material parameter

*k*

_{1}has the dimension of stress and ν is a dimensionless parameter. The anisotropic term contributes to the mechanical response only when the fibres are extended, in other words when either

*I*

_{4}>1 or

*I*

_{6}>1 or both. The cylinder is subject to an internal pressure

*P*and, therefore, the axial equilibrium and the boundary conditions at the inner and outer surfaces yield 8.9 8.10

The pressure/axial stretch plot shown in figure 11, with parameters from Holzapfel (2004), reveals an interesting trend. For a particular axial stretch, the length does not change as the pressure increases. This result is consistent with the inversion line found in Schulze-Bauer *et al.* (2003). We now perform a bifurcation analysis and determine the effective modulus of the grown, pressurized cylinder. It should be noted that the parameter values *a*_{0}, *c*_{0} and *l*_{0} are the geometric parameters of the grown, pressurized cylinder, and these values are used in equation (5.3) to calculate the effective modulus. We consider a constant growth in the radial and circumferential directions in the inner layer such that *γ*_{ri}=*γ*_{θi}, *γ*_{ro}=*γ*_{θo}=1 and *γ*_{zi}=*γ*_{zo}=1. Both growth and an applied internal pressure provide substantial improvements in the overall rigidity of the anisotropic cylindrical shell (figure 12). It is interesting to note the effect of anisotropy on the overall stability of the structure. In figure 12, it is shown that, in an isotropic material (with values taken from Holzapfel *et al.* (2002)), the internal pressure slightly decreases the overall rigidity. However, in an anisotropic material, the internal pressure dramatically increases the effective modulus. We conclude that anisotropy creates an environment in which the internal pressure and differential growth are advantageous.

## 9. Conclusion

In this paper, we have explored the role of differential growth and residual stress in cylindrical elastic structures. We have presented a general formulation of growth and a bifurcation analysis for a three-dimensional nonlinear elastic body and applied it to cylindrical growth, which is relevant in many physiological and biological systems. The analysis presented was performed using the general framework of nonlinear elasticity that can be easily generalized to include various effects such as anisotropic response, as shown here.

The emphasis was on the change of rigidity induced by growth differences between the inner and the outer layer in a cylindrical tube. By itself, inhomogeneous axial differential growth does not notably change the overall rigidity of a cylindrical structure. Similarly, in the absence of growth, a strain-stiffening material in the outer layer will not change the overall structural rigidity. During differential growth, the outer layer is in a state of tensile stress. Therefore, the incremental modulus is correspondingly very large, and the overall rigidity of the structure is greatly enhanced. Similarly, the coupling between radial growth and anisotropic response also greatly enhances the overall flexural rigidity.

## Acknowledgements

This material is based, in part, upon work supported by the National Science Foundation under grant no. DMS-0604704 (A.G.) and an ARCS fellowship (R.V.).

## Footnotes

One contribution of 12 to a Theme Issue ‘Mechanics in biology: cells and tissues’.

- © 2009 The Royal Society