A refined dynamic finite-strain shell theory for incompressible hyperelastic materials: equations and two-dimensional shell virtual work principle
AA refined dynamic finite-strain shell theory forincompressible hyperelastic materials: equationsand two-dimensional shell virtual work principle
Xiang Yu , Yibin Fu , Hui-Hui Dai , ˚ Department of Mathematics, City University of Hong Kong, Kowloon Tong, HongKong School of Computing and Mathematics, Keele University, Keele, UK
Abstract
Based on previous work for the static problem, in this paper we first derive one form ofdynamic finite-strain shell equations for incompressible hyperelastic materials that involvethree shell constitutive relations. In order to single out the bending effect as well as toreduce the number of shell constitutive relations, a further refinement is performed, whichleads to a refined dynamic finite-strain shell theory with only two shell constitutive rela-tions (deducible from the given three-dimensional (3D) strain energy function) and somenew insights are also deduced. By using the weak formulation of the shell equations andthe variation of the 3D Lagrange functional, boundary conditions and the two-dimensional(2D) shell virtual work principle are derived. As a benchmark problem, we consider the ex-tension and inflation of an arterial segment. The good agreement between the asymptoticsolution based on the shell equations and that from the 3D exact one gives verification ofthe former. The refined shell theory is also applied to study the plane-strain vibrations ofa pressurized artery, and the effects of the axial pre-stretch, pressure and fibre angle onthe vibration frequencies are investigated in detail.
In recent years, biological materials have attracted a lot of interest; see, for example, thereview article by Holzapfel and Ogden [1] on constitutive modelling of arteries. There aretwo noteworthy properties of biological materials. materials. One is they are very soft andcan undergo large elastic deformations with finite-strain; the other is that the volume is pre-served during the deformation. So, they are normally modelled as incompressible hyperelasticmaterials. Many biological tissues and organs are thin structures. Due to the complexity ofthe 3D formulation and the cost and ineffectiveness of 3D computations (in particular, forpost-bifurcation solutions), often one needs to use a 2D shell model to study their behaviors. ˚ Corresponding author
Email address : [email protected] The work described in this paper is fully supported by a GRF grant (Project No.: CityU 11303718) fromthe Research Grants Council of the Government of HKSAR, P.R. China. a r X i v : . [ c s . C E ] A ug hell theories have a long history, which date back to the pioneering work of Love [2] in1888. Since then, they have been studied extensively during the past 130 years. Numerousworks on shell theories have been done in the framework of linearized elasticity and/or linearconstitutive relation with geometric nonlinearity. Here, the focus is on soft materials modelledby a strain energy function with incompressibility constraint, for which one needs to considermaterial nonlinearity. It is out of the scope of the present study to give an extensive review onlinear shell theories or those with geometric nonlinearity, and for a selected review, we referto Li et al. [3]. Instead, we only give a review on derived shell theories for incompressiblehyperelastic materials , for which, relatively speaking, there are not so many works.In [4], Makowski and Stumpf formulated a finite-strain shell theory for incompressible hy-perelastic materials by assuming the material lines normal to the shell surface remain straightduring the deformation. Itskov [5] assumed that the position vector in the deformed shell islinear in the thickness variable (with six parameters). The incompressibility constraint is usedto eliminate the transverse normal strain, and based on which, a numerical shell theory withfive parameters for a generalized orthotropic incompressible hyperelastic material was devel-oped. In [6], Chapelle et al. examined whether the plane stress assumption or the asymptoticlimits of thickness can commute with the incompressibility constraint, justifying the usagesof classical shell models and a modified 3D shell model in the incompressible conditions. InKiendl et al. [7], a shell theory for compressible and incompressible isotropic hyperelasticmaterials was developed based on the Kirchhoff-Love kinematics which includes the assump-tions of zero transverse normal stress and straight and normal cross-sections, and then anisogeometric discretization was introduced for numerical computation. Recently, Amabili etal. [8], for a tube (a special kind of shells), developed a shell theory for incompressible biolog-ical hyperelastic materials by assuming the in-plane displacement components are third-orderpolynomials of the thickness variable while the out-plane component is a fourth-order polyno-mial. Further simplification in that work includes the dropping of certain nonlinear terms inthe strain-displacement relations and incompressibility condition, which enables one to rep-resent the four coefficients in the out-plane displacement in terms of other unknowns. As aresult, a nine-parameter shell theory was obtained. All the works mentioned above employ ad hoc assumptions and cross-thickness integrations to eliminate the thickness variable. Asa result, one cannot expect that the resulting shell theories are consistent with the 3D fieldequations, top and bottom traction conditions and incompressibility condition in a pointwisemanner. It is difficult to assess the reliability of such inconsistency for general loading. Also,when higher-order expansions are used, higher-order resultants need to be introduced buttheir physical meanings are not clear. Thus, it is more desirable to construct a shell theorywithout these ad hoc assumptions/simplifications, which is consistent with the 3D formula-tion (field equations and top/bottom traction condition and incompressibility constraint) toa proper asymptotic order in a pointwise manner.We also mention that by the Γ-convergence method, Li and Chermisi [9] rigorously derivedthe von K´arm´an shell theory for incompressible hyperelastic materials. However, this kindof approach depends on some a priori scaling assumptions, which cannot yield a shell theorywith both stretching and bending effects.In a recent paper of Dai and Song [10], a dimension-reduction method was proposed to2onstruct a consistent plate theory with both stretching and bending effects via series expan-sions with only smoothness assumption (without any ad hoc kinematic or other assumptions).The idea is to directly work with the 3D field equations and traction conditions on the topand bottom surfaces, and then to establish some recurrence relations for the expansion coeffi-cients. Then, the approach has been used to derive a dynamic plate theory [11], a static shelltheory [12], a static plate theory for incompressible materials [13] and a static shell theory forincompressible materials [3].In this paper, we follow Dai and Song’s approach to first derive one form of dynamic shelltheory for incompressible hyperelastic materials that involves three shell constitutive relationsand six boundary conditions at each edge point. The completely new part is on the furtherrefinement by elaborate calculations (cf. the procedure for a plate in [14]), which reducesthe number of shell constitutive relations to two and singles out the bending term. It turnsout that the refined shell equations alone can reveal a few new insights already. For theforce boundary, in practice one only knows four conditions: the bending moment along theedge tangent direction and the three components of the cross-thickness resultant. To proposeproper boundary conditions, we incorporate the weak form of the refined shell equationsinto the variation of the 3D Lagrange functional δL . By some elaborate calculations, whichprovide guidance on choosing the variation of the displacement vector in the 3D edge term in δL when specializing to a 2D shell theory, suitable shell boundary conditions and the 2D shellvirtual work principle are obtained. A benchmark problem of an artery segment subjected toextension and internal pressure is considered. Finally, as an application of the refined shelltheory, the plane-strain vibrations of a pressurized artery are studied, and the results revealthe influences of the axial pre-stretch, pressure, and fibre angle on the vibration frequencies. Notation.
Throughout this paper, we use boldface letters to denote vectors and second-order tensors; we use curly letters to denote higher-order tensors. The summation conventionfor repeated indices is adopted. In a summation, Greek letters α, β, γ, . . . run from 1 to 2,whereas Latin letters i, j, k, . . . run from 1 to 3. A comma preceding indices means differen-tiation and a dot over variables indicates time derivative. The time argument in variables isusually omitted for brevity.Let R be the three-dimensional Euclidean space with standard basis p e , e , e q . Thesymbol I : “ e i b e i is reserved for the identity tensor of R . The notation ^ means crossproduct. For a scalar-valued function of a tensor W p F q , the derivative of the W with respectto F is defined to be B W B F : “ B W B F ji e i b e j ; higher-order derivatives are defined in a similar way.The divergence of a tensor S is defined by Div p S q : “ B S ij B x i e j . The tensor contractions aredefined by A r B s “ tr p AB q : “ A ji B ij , A r A s : “ A ij(cid:96)k A k(cid:96) e i b e j , A r a , b s : “ Aa ¨ b “ A ij a j b i . (1.1) We consider a thin shell of constant thickness 2 h composed of an incompressible hyperelas-tic material which occupies a region Ω ˆ r , h s in the reference configuration. The thickness2 h of the shell is assumed to be small compared with the length scale of the bottom surface3 and its ratio against the radius of curvature is less than 1. The position of a materialpoint is denoted by X in the reference configuration and by x in the current configuration.The geometric description of a shell has been given in [15] and [16], and here we give a briefsummary.The bottom surface Ω of the shell is parameterized by two curvilinear coordinates θ α , α “ ,
2. The position of a point on Ω is written as r “ r p θ α q . Then the tangent vectorsalong the coordinate lines are given by g α “ B r {B θ α , which form a covariant basis of thetangent plane of the bottom surface. Their contravariant counterparts g α , which satisfy therelations g α ¨ g β “ δ αβ , form a contravariant basis of the same plane. The unit normal vector n to the bottom surface is defined via n “ g ^ g {| g ^ g | , so that by setting g “ g “ n , t g i u and t g i u , i “ , , X “ r p θ α q ` Z n p θ α q , ď Z ď h, (2.1)where Z is the coordinate of the point along the normal direction n . The change of theunit normal vector is captured by the curvature map, which is defined as the negative of thetangent map of the Gauss map n : Ω Ñ S [15], where S denotes the two-dimensional unitsphere; thus we have k “ ´B n {B r “ ´ n ,α b g α . We point out that the curvature tensor k issymmetric in the sense that k “ k T . Associated to k , the mean curvature and the Gaussiancurvature are respectively defined by H “ tr p k q and K “ det p k q .The covariant basis vectors at a point in the shell Ω ˆ r , h s are given by p g α “ B X B θ α “ B r B θ α ` Z B n B r B r B θ α “ p ´ Z k q g α , (2.2)where : “ I ´ n b n “ g α b g α denotes the projection onto the tangent plane of Ω; it is alsothe identity map of the same plane. Setting µ “ ´ Z k , we see from (2.2) that p g α “ µg α and thus p g α “ µ ´ T g α . Note that the previous geometric assumption which asserts | hk βα | ă µ ´ is well-defined. By the change of variables formula, the volumeelement of the shell is computed by dV “ p p g ^ p g q ¨ n dθ dθ dZ “ det p µ qp g ^ g q ¨ n dθ dθ dZ “ µ p Z q dAdZ, (2.3)where µ p Z q “ det p µ q “ ´ HZ ` KZ and dA “ | g ^ g | dθ dθ is the area element on thebottom surface.On the boundary B Ω, let s be the arc length variable, and let τ and ν be respectivelythe unit tangent vector and the unit outward normal vector such that p τ , n , ν q forms a right-handed triple (i.e., ν “ τ ^ n ). Then let N , T and da be respectively the unit outwardnormal vector, unit tangent vector and the area element of the lateral surface such that p T , n , N q forms a right-handed triple. A similar argument as in (2.2) yields T “ p ´ Z k q τ {? g τ , where ? g τ denotes the magnitude of vector p ´ Z k q τ and is given by ? g τ “? ´ Z kτ ¨ τ ` Z kτ ¨ kτ . Using the change of variables formula again, we have N da “ µτ ds ^ n dZ “ p ´ Z k q τ ^ n dsdZ . Then from the equality p kτ q^ n “ tr p k qp τ ^ n q´ k p τ ^ n q ,we deduce that N da “ p ` Z p k ´ H qq ν dsdZ. (2.4)4ince p ´ Z k q τ “ ? g τ T and p T , n , N q forms a right-handed triple of unit vectors, we have da “ ? g τ dsdZ and ? g τ N “ p ` Z p k ´ H qq ν from the above equations.The deformation gradient is then calculated by F “ B x B X “ B x B θ α b p g α ` B x B Z b n “ p ∇ x q µ ´ ` B x B Z b n , (2.5)where ∇ : “ BB θ α g α denotes the 2D gradient operator on the base surface Ω. We remark thatfor the 2D gradient operator, one has the following Stokes’ theorem ż Ω ∇ ¨ p a q dA “ ż B Ω a ¨ ν ds, ż Ω ∇ ¨ p S q dA “ ż B Ω S T ν ds (2.6)for a vector field a and a tensor field S , respectively.For an incompressible material, one has the following incompressibility constraint R p F q “ det p F q ´ “ . (2.7)Assume further that the material is hyperelastic with a strain energy function W p F q . Thenthe associated elastic moduli are defined by A i p F q “ B i ` W B F i ` , i “ , , . . . . The strain energyfunction is assumed to satisfy the strong-ellipticity condition: p A p F qr a b b sqr a b b s ą a b b ‰ q ` and q ´ are the external loads applied on the top and the bottom surfacesof the shell respectively. The boundary B Ω of the bottom surface Ω is divided into two parts:the position boundary B Ω subjected to the prescribed position b and the traction boundaryΩ q subjected to the applied traction q . Then the kinetic energy K , the strain energy S , andthe load potential V of the shell are respectively given by K “ ż Ω ż h ρ x ¨ x µ p Z q dZdA, S “ ż Ω ż h W p F q µ p Z q dZdA, (2.8) V “ ´ ż Ω p q ´ p r q ¨ x p r , q ` q ` p r q ¨ x p r , h q µ p h qq dA ´ ż Ω ż h q b ¨ x µ p Z q dZdA ´ ż B Ω q ż h q p s, Z q ¨ x p s, Z q da, (2.9)where ρ is the mass density of the shell, q b is the body force and da is the area element onthe lateral surface B Ω ˆ r , h s .By Hamilton’s principle, the 3D equations are obtained when the energy functional E “ K ` S ` V attains its minimum under the constraint condition (2.7). Therefore we are led toconsider the Lagrange functional L p x p X q , p p X qq “ K ` S ` V ´ ż Ω ż h p p X q R p F q µ p Z q dZdA, (2.10)where p p X q is the Lagrange multiplier. To attain the minimum, it is necessary that thevariation of L with respect to x is zero, and a direct calculation shows δL “ ż Ω ż h p ρ : x ´ Div p S q ´ q b q ¨ δ x µ p Z q dZdA ´ ż Ω p S T n | Z “ ` q ´ q ¨ δ x p r , q dA ` ż Ω p S T n | Z “ h ´ q ` q ¨ δ x p r , h q µ p h q dA ` ż B Ω q ż h p S T N ´ q q ¨ δ x p s, Z q da “ , (2.11)5here S “ B W B F ´ p F ´ (2.12)is the nominal stress tensor of the incompressible hyperelastic material [17], and we used B R {B F “ det p F q F ´ and det p F q “
1. Since δ x in (2.11) is arbitrary, we obtain the following3D momentum equations together with boundary conditions:Div p S q ` q b “ ρ : x in Ω ˆ r , h s , (2.13) S T n | Z “ “ ´ q ´ , S T n | Z “ h “ q ` in Ω , (2.14) S T N “ q p s, Z q on B Ω q ˆ r , h s , (2.15) x “ b p s, Z q on B Ω ˆ r , h s . (2.16)The above equations together with the incompressibility constraint (2.7) form the 3D dynamicequations for the shell structure, which contain an independent vector variable x and anindependent scalar variable p . In this section, we shall first derive one form of consistent shell equations with threeshell constitutive relations. Here the consistency means each term in (2.11) should be of arequired asymptotic order, separately for the approximation. Then, a refinement is performedto reduce the number of shell constitutive relations from three to two. Also, the bending termis singled out. For the first part, the derivation is similar to that of the static case [3], but tobe self-contained, we present the main steps.
We assume sufficient smoothness for the quantities involved. Then x p X q , p p X q , F p X q and S p X q have Taylor expansions about the bottom surface Z “
0. From (2.5) and the nonlinearrelation (2.12), the following relations among their expansion coefficients can be found: F p q “ ∇ x p q ` x p q b n , F p q “ ∇ x p q k ` ∇ x p q ` x p q b n , (3.1)and S p q “ A ´ p p q F p q´ , S p q “ A r F s ` p p q F p q´ F p q F p q´ ´ p p q F p q´ , (3.2)where the superscript p i q denotes the i th derivative with respect to Z at Z “
0, and A “B W {B F | F “ F p q . From the above expressions, one easily checks that S p i q is linear algebraicin p p i q and x p i ` q , i “ i “
2; for brevity the relations for F p q and S p q areomitted). It is due to this linearity that some recurrence relations can be established for theexpansion coefficients upon further using the field equations in the subsequent derivations. Remark 3.1.
The expressions for S p i q p i “ , , q give three relations between the stresscoefficients and the position vector coefficients. In the sequel, we abuse the terminology alittle and call equations (3.2) and (3.2) and that for S p q to be shell constitutive relations.The reason is that the derived shell equations are represented in terms of S p i q and throughthese relations the unknown in the shell equations is actually the position vector x p q .6ow, we shall proceed to do the dimension reduction process by using the 3D formulation.First, the bottom traction condition (2.14) yields S p q T n “ p A ´ p p q F p q´ q T n “ ´ q ´ . (3.3)To ease notation, we introduce the vector g “ g p x p q q “ det p F p q q F p q´ T n “ x p q , ^ x p q , { a | g ^ g | (see [3]). Then noting (3.1) and det p F p q q “
1, the above equation can be written as p A p ∇ x p q ` x p q b n qq T n “ ´ q ´ ` p p q g , (3.4)Next, substituting the Taylor expansion for S into the field equation (2.13) and equatingthe coefficients of Z i ( i “ , , . . . ) on both sides, we have ∇ ¨ S p q ` S p q T n ` q p q b “ ρ : x p q , (3.5) ∇ ¨ S p q ` S p q T n ` p kg α q ¨ S p q ,α ` q p q b “ ρ : x p q , (3.6)where ∇ ¨ S : “ g α ¨ S ,α denotes the 2D divergence of the tensor S . Then substituting theTaylor expansion for F into the constraint equation (2.7) and equating the coefficients of Z i to be zero, we obtain g ¨ x p q ´ “ , (3.7) g ¨ x p q ` tr p F p q´ p ∇ x p q k ` ∇ x p q qq “ , (3.8)where in (3.7) we have used the equality F p q´ x p q “ n implied by (3.1) . By the way, wepoint out that there is a typo in p q of [3].With the use of (3.2) , equation (3.5) can be simplified into Bx p q ` f ´ p p q g “ ρ : x p q (3.9)by defining Ba “ p A r a b n s ` p p q F p q´ p a b n q F p q´ q T n ðñ B ij “ A i j ` p p q F p q´ i F p q´ j , (3.10) f “ p A r ∇ x p q k ` ∇ x p q s ` p p q F p q´ p ∇ x p q k ` ∇ x p q q F p q´ q T n ` ∇ ¨ S p q ` q p q b . (3.11)From (3.8) and (3.9), we obtain p p q “ g ¨ B ´ g p g ¨ B ´ p f ´ ρ : x p q q ´ tr p F p q´ p ∇ x p q k ` ∇ x p q qqq , (3.12) x p q “ B ´ p p p q g ` ρ : x p q ´ f q . (3.13)Note that the strong-ellipticity condition guarantees that B is positive definite and hence isinvertible. The explicit expressions of x p q and p p q can be obtained similarly, whose expres-sions are omitted. The explicit expressions of x p q and p p q are not needed since they areintermediate variables. The explicit expressions for x p q and p p q are encoded in (3.4) and(3.7), which are nonlinear algebraic equations in general, so they can only be solved when thestrain energy function is specified. Nevertheless, the strong-ellipticity condition together with7he implicit function theorem ensures that x p q and p p q can be uniquely solved in terms of x p q (cf. [13]).Finally, the top traction condition (2.14) states S p q T n ` h S p q T n ` h S p q T n ` h S p q T n ` O p h S p q T n q “ q ` . (3.14)Subtracting (3.14) multiplied by µ p h q “ ´ Hh ` Kh from (3.3) and then simplifying(see [11] for details), we arrive at one form of a 2D dynamic vector shell equation ∇ ¨ r S ` O p h S p q , h k S p q q “ ρ : r x ´ r q ` O p h : x p q , h k : x p i q , h q p q b , h k q p i q b q , (3.15)where i “ , r S “ p ` h p k ´ H qq S p q ` h p ` h p k ´ H qq S p q ` h S p q “ h ż h p ` Z p k ´ H qq S dZ ` O p h S p q , h k S p q q , (3.16) r x “ p ´ hH ` h K q x p q ` h p ´ hH q x p q ` h x p q “ h ż h x µ p Z q dZ ` O p h x p q , h k x p i q q , (3.17) r q “ µ p h q q ` ` q ´ h ` r q b , (3.18)and r q b is defined in the same way as r x . Remark 3.2.
The quantity r S is considered as the averaged stress, and r q the averaged shellbody force due to surface traction and 3D body force. We point out that (3.15) can be alsodeduced by multiplying the field equation (2.13) by µ p Z q and then integrating it with respectto Z from 0 to 2 h followed by applying the equality ż h Div p S q µ p Z q dZ “ ∇ ¨ p ż h p ` Z p k ´ H qq S dZ q ` S T n | Z “ h µ p h q ´ S T n | Z “ , (3.19)which is a consequence of Stokes’ theorem.Similar to [12], suitable edge boundary conditions can be imposed, and then it can beshown that each of the five terms in (2.11) is of O p h q , which satisfies the consistency crite-rion. The details are omitted. Also, it is clear from the derivation process that the bottomtraction condition, the 3D field equations, the incompressibility condition and the top tractioncondition are all satisfied in a pointwise manner (with an error of O p h q , see (3.14)), an impor-tant feature not enjoyed by shell theories based on ad hoc assumptions and/or cross-thicknessintegrations. Although the above-derived shell theory is consistent, there are still a few undesirablefeatures as follows. 1. There are a little too many (three) shell constitutive relations (equations(3.2) and (3.2) and that for S p q ). In particular, the relation between S p q and x p q is8ery complicated and can cause some technical difficulties for implementation in concreteapplications. 2. From the shell equations, one cannot tell clearly which term(s) represents thebending effect. 3. Although the associated weak form can be obtained from the shell equations,physically it does not represent the shell virtual work principle. 4. The shell equations arethree coupled fourth-order PDEs for x p q , which require six boundary conditions at an edgepoint. If one knows the displacement and/or stress distributions, there is no difficulty imposingthem. However, in many practical situations for the traction edge, one only knows fourconditions: the cross-thickness force resultant and the bending moment (with direction alongthe edge tangent), and one does not know how to impose the other two boundary conditions.For a plate theory, these issues were addressed in [14]. Here, with some modifications, thoseideas from this previous work will be used for a shell theory. In this subsection, we shallresolve the first two issues by performing some manipulations to eliminate S p q and to singleout the bending term. As a price to pay, the relative errors for some problems may not beas good as before. We point out that one cannot simply drop h S p q in (3.16), as thebending effect is also dropped. So, one needs to do some elaborate calculations to extract thebending term first and then to drop the relative higher-order terms. The last two issues willbe resolved in the next section.First, we rewrite (3.15) into two parts: ∇ ¨ r S ` O p h S p q , h k S p q q “ ρ : r x t ´ r q t ` O p h : x p i q t , h k : x p i q t , h q p i q bt , h k q p i q bt q , (3.20) p ∇ ¨ r S q ¨ n ` O p h S p q , h k S p q q “ ρ : r x ´ r q ` O p h : x p q , h k : x p i q , h q p q b , h kq p i q b q , (3.21)where “ I ´ n b n “ g α b g α and the subscript t indicates the projection into the tangentplane; thus a t : “ a “ a and S t : “ S for a vector a and a tensor S respectively. Notethat since r S satisfies the equality r S “ r S (see (3.16)), we have r S t “ r S “ r S . (3.22)Next, we want to extract terms related to in-plane stress r S t from the in-plane equation(3.20) in order to gain some insights as well for later use for deriving the 2D shell virtual workprinciple. For this purpose, we need the following two equalities for a tensor field S and avector field a : ∇ ¨ S “ ∇ ¨ p S q ´ k αβ S β g α , (3.23) p ∇ ¨ S q ¨ a “ ∇ ¨ p Sa q ´ tr p ∇ aS q . (3.24)To prove (3.23), it suffices to show that ∇ ¨ p S ´ S q “ ´ k αβ S β g α . (3.25)Since “ I ´ n b n , we have S ´ S “ Sn b n . Further, we have ∇ ¨ p Sn b n q “ p g β ¨ p Sn b n q ,β q “ g β ¨ p Sn q ,β n ` p g β ¨ Sn q n ,β (3.26) “ ´p g β ¨ Sn q kg β “ ´ k αβ S β g α . (3.27)Thus (3.23) follows. Equation (3.24) can be proved by a direct calculation starting from ∇ ¨ p Sa q by using the definition of the 2D divergence.9sing (3.23), (3.24) and (3.22), and noting that ∇ n “ ´ k , (3.20) and (3.21) can berewritten as ∇ ¨ r S t ´ k αβ r S β g α ` O p h S p q , h k S p q q “ ρ : r x t ´ r q t ` O p h : x p q t , h k : x p i q t , h q p q bt , h k q p i q bt q , (3.28) ∇ ¨ p r Sn q ` tr p k r S t q ` O p h S p q , h k S p q q “ ρ r x ´ r q ` O p h : x p q , h k : x p i q , h q p q b , h kq p i q b q . (3.29)Now, we shall manipulate the third equation (3.29) further to single out the bending term.Adding (3.14) multiplied by µ p h q to (3.3), we obtain p ´ hH ` h K q S p q T n ` h p ´ hH q S p q T n ` h S p q T n ` O p h S p q T n , h k S p i q T n q “ m , (3.30)where i “ , m “ p µ p h q q ` ´ q ´ q{
2. To extract the bending term from (3.29), wesubtract the 2D divergence of (3.30) multiplied by from the left from (3.29) (with thesubstitution of (3.16)). Note that the focus for this manipulation is on the S p q terms in thesetwo equations. Then, upon further using (3.3) and (3.6), we obtain ∇ ¨ pp ` h p k ´ H qq S p q n ´ pp ´ hH ` h K q ` h k q S p q T n q` h ∇ ¨ pp ` h p k ´ H qq S p q n ´ p ´ hH q S p q T n q ` h ∇ ¨ p S p q n ´ S p q T n q` tr p k r S t q ` h ∇ ¨ p pp kg α q ¨ S p q ,α qq ` h ∇ ¨ p ∇ ¨ S p q q ` O p h S p q , h k S p i q q“ ρ : r x ´ r q ` h ∇ ¨ p ρ : x p q t ´ q p q bt q ´ ∇ ¨ m t ` h ∇ ¨ p kq ´ t q ` O p h : x p q , h k : x p i q , h q p q b , h kq p i q b q . (3.31)We also want to extract the in-plane stress parts of the last term h ∇ ¨ p ∇ ¨ S q on theleft-hand side. Observe that we have the decomposition S p q “ IS p q I “ p ` n b n q S p q p ` n b n q “ S p q t ` n b S p q T n ` S p q n b n . (3.32)Further, routine calculations show that ∇ ¨ p ∇ ¨ p n b S p q T n qq “ ´ ∇ ¨ p H S p q T n qq , (3.33) ∇ ¨ p ∇ ¨ p S p q n b n qq “ ´ ∇ ¨ pp g α ¨ S p q n q kg α q . (3.34)Upon using the above three equations, (3.31) can be recast as ∇ ¨ pp ` h p k ´ H qq S p q n ´ pp ´ hH ` h K q ` h k q S p q T n q` h ∇ ¨ pp ` h p k ´ H qq S p q n ´ p ´ hH q S p q T n q ` h ∇ ¨ p S p q n ´ S p q T n q` tr p k r S t q ` h ∇ ¨ p ∇ ¨ S p q t q ´ h ∇ ¨ p H S p q T n q ´ h ∇ ¨ pp g α ¨ S p q n q kg α q` h ∇ ¨ p pp kg α q ¨ S p q ,α qq ` O p h S p q , h k S p i q q“ ρ : r x ´ r q ` h ∇ ¨ p ρ : x p q t ´ q p q bt q ´ ∇ ¨ m t ` h ∇ ¨ p kq ´ t q ` O p h : x p q , h k : x p i q , h q p q b , h kq p i q b q . (3.35)10o eliminate S p q terms in a consistent manner, we shall drop any term which is relatively O p h q or O p h q smaller than another term (so that the shell theory yields results with arelative O p h q or O p h q error). It is justified, as shown in the following simple example: for A ` B ` C “
0, if C “ O p h B q or C “ O p hB q , the dropping of C causes at most a relativeerror of O p h q or O p h q , no matter A ą O p B q or A ď O p B q . Any terms which cannot satisfythe above requirement will be kept.We make the following observations. 1. In (3.28), h S p q in r S (cf. (3.16)) is dropped,as it is O p h q smaller than S p q or O p h q smaller than h S p q if S p q “ F p q “ R and thus S p q “
0, where R is a rotation tensor). As it is possible that S p q terms become the leading ones, theyshould be kept. 2. The last three terms on the left-hand side of (3.35), h ∇ ¨ p K S p q T n q , h ∇ ¨ pp k ´ H q S p q n q and h ∇ ¨ p H S p q T n q are dropped as they are O p h q smaller thantr p k r S t q or either O p h q smaller than tr p k r S t q or zero if S p q “
0. 3. The third term on theleft-hand side of (3.35) is dropped as it is O p h q smaller than ∇ ¨ p S p q n ´ S p q T n q or O p h q smaller than h ∇ ¨ p S p q n ´ S p q T n q if S p q “
0. 4. On the right-hand sides, h x p q in r x (cf. (3.17)) is dropped , as it is O p h q smaller than x p q , and a similar treatment is made to r q b . From these observations, we have the refined 2D dynamic shell equations as follows: ∇ ¨ S t ´ k αβ S β g α “ ρ : x t ´ q t , (3.36) ∇ ¨ p S ‹ n ´ S T ‹ n q ` tr p kS t q ` h ∇ ¨ p ∇ ¨ S p q t q“ ρ : x ´ q ` h ∇ ¨ p ρ : x p q t ´ q p q bt q ´ ∇ ¨ m t ` h ∇ ¨ p kq ´ t q , (3.37)where S “ p ` h p k ´ H qq S p q ` h p ` h p k ´ H qq S p q , (3.38) S ‹ “ p ` h p k ´ H qq S p q ` h S p q , (3.39) S T ‹ “ p ` h p k ´ H qq S p q T ` h S p q T , (3.40) x “ p ´ hH ` h K q x p q ` h p ´ hH q x p q , (3.41) q “ µ p h q q ` ` q ´ h ` q b , (3.42)and q b is defined in the same way as x .From the above shell equations, one can observe some important insights. 1. For a plate(or a shell with | k αβ | ď O p h q ) in linear elasticity, the bending term h ∇ ¨ p ∇ ¨ S p q t q becomesthe leading term, so it should be kept although it looks like an O p h q term. 2. For the in-planeequation (3.36), the in-plane forces and inertia effects are resisted by two sources: the in-planestress part (the first term on the left-hand side) and the out-plane shear stresses due to thecurvature effect (the second term). 3. For the out-plane equation (3.37), the out-plane forcesand inertia effects are resisted by three sources: (i) the out-plane shear stresses (the first termon the left-hand side) due to geometric and/or material nonlinearity; (ii) the in-plane stressesdue to the curvature effect (the second term); (iii) bending effect due to the in-plane stresses(the last term). 4. Although the out-plane normal stress does not appear explicitly in theseshell equations, it plays a role in expressing x p q and p p q in terms of x p q (see (3.3) and (3.7)),11o it should not be ignored (as in some ad hoc theories, which assume the out-plane componentof the displacement is independent of Z ). 5. Only two shell constitutive relations are needed,which are provided by (3.2) and (3.2) . 6. These shell equations provide results with at mosta relative O p h q error, although in some cases the error can be O p h q . Note that higher-orderTaylor expansions do not necessarily lead to higher-order correct plate/shell equations.After substitutions of all recurrence relations, the above shell equations become a systemof differential equations involving x p q only. Once it is solved, x p q (with a relative error equalto or smaller than O p h q ) is obtained and the position vector x can then be recovered. Now we shall resolve the last two issues mentioned in the beginning of the previous subsec-tion. Actually, boundary conditions for a derived shell theory can cause considerable difficulty(see Steigmann [18]). Here, we shall use both the variation of the 3D Lagrange functional andthe weak form of the shell equations to get the appropriate boundary conditions and the 2Dshell virtual work principle.For the shell equations, the bottom traction condition (2.14) , and the vanishing coeffi-cients of the field equation (2.13) and the incompressibility constraint (2.7) are used to findthe recurrence relations. As a result, (2.13) (up to required order) and (2.14) can be treatedas identities. To obtain the 2D shell virtual work principle from the vanishing of the variationof 3D Lagrange functional (2.11), we need to specialize it to the 2D case (by using the Taylorexpansions for the quantities involved as in deriving the shell equations). The first two termsin (2.11) can be set to be identically zero because of the above-mentioned two identities.Then, in order to remove δ x p r , h q (we still use x p r , h q for the writing purpose but it meansthe Taylor expansion of the position vector at Z “ h ) in the third integral and to introduce δ x p r , h q to the variation (needed for the 2D shell virtual work principle), we add to δL threeidentically zero terms (the first three terms below) to obtain δL “ h ż Ω A t ¨ p δ x t p r , h q ´ δ x t p r , h qq dA ` h ż Ω A ¨ p δx p r , h q ´ δx p r , h qq dA ` h ż Ω p ∇ ¨ p C qq ¨ δ x p r , h q dA ` ż Ω p S T n | Z “ h ´ q ` q ¨ δ x p r , h q µ p h q dA ` ż B Ω q ż h p S T N ´ q q ¨ δ x p s, Z q da “ , (4.1)where A t “ A “ C “ δ x p r , h q terms canceleach other (upon dropping relatively higher-order terms as in Section 33.2), and thus we have δL “ ´ h ż Ω A t ¨ δ u mt dA ´ h ż Ω A ¨ δu m dA ` ż B Ω q ż h p S T N ´ q q ¨ δ u p s, Z q da “ , (4.2)where we have used the virtual displacement δ u to replace the virtual position vector and thesubscript m denotes the middle surface Z “ h . Actually, the first two terms are just the weak12orm for the shell equations (3.36) and (3.37). We remark that when the boundary conditionsare involved, one can only expect to obtain the leading-order results in general; thus in thesequel, any term, which is relatively smaller than another term, will be dropped.To get the 2D shell virtual work principle, we shall further add two identities to the aboveequation, which are associated with the virtual work due to the moment, which is given by M “ ż B Ω ż h pp x ´ x p r , h qq ˆ S T N q? g τ dZds. (4.3)Then, the twist moment (along N m direction) and the bending moment (along T m direction)per unit arc length of B Ω are given by respectively T “ ż h pp x ´ x p r , h qq ˆ S T N q ¨ N m ? g τ dZ “ h S p q T r ν , ν ˆ x p q s ` h S p q T r ν , ν ˆ x p q s ` O p h , h k q , (4.4) M “ ż h pp x ´ x p r , h qq ˆ S T N q ¨ T m ? g τ dZ “ h S p q T r ν , τ ˆ x p q s ` h S p q T r ν , τ ˆ x p q s ` O p h , h k q . (4.5)It was shown in [19] (Section 2.5; the authors attributed the argument to Kirchhoff) thatthe derivative of the twisting moment with respect to the arc length T ,s is equivalent to adistributed shear force (along the downward thickness direction). Thus, this twist momentgenerates a virtual work per unit arc length: ´ T ,s δu m (the smoothness of B Ω is assumed). Onthe other hand, the bending moment generates a virtual work per unit arc length: ´ M δα m ,where α m is the rotation angle at the edge of the middle surface. It is defined as the change ofthe angle between the vector N m and the projected vector onto the nN m -plane of the tangentvector at an edge point of the intersection curve of the middle surface and the nN m -planeduring the deformation, which is given by (after some calculations) α m “ arctan p ∇ m x p r , h qr N m s ¨ n ∇ m x p r , h qr N m s ¨ N m q ˘ pπ “ arctan p u m ,ν ` ∇ u mt r ν , ν s ` O p k, hk qq ˘ pπ, (4.6)where ∇ m “ BB θ α p g α | Z “ h (see (2.2) for the definition of p g α ) is the gradient operator on themiddle surface and p is a natural number.Now, we add the two identities ´ T ,s δu m ` T ,s δu m “ ´ M δα m ` M δα m “ δL “ ´ h ż Ω A t ¨ δ u mt dA ´ h ż Ω A ¨ δu m dA ` ż B Ω q ż h p S T N ´ q q ¨ δ u p s, Z q da ` ż B Ω T ,s δu m ds ´ ż B Ω T ,s δu m ds ` ż B Ω M δα m ds ´ ż B Ω M δα m ds “ . (4.7)Next, substituting the expressions of A t and A according to the shell equations (3.36)and (3.37) into the above equation and then doing integration by parts by Stokes’ theorem,13e obtain, after dropping O p h , h k q terms,2 h ż Ω p tr p S t ∇ δ u mt q ` k αβ S β g α ¨ δ u mt ` p ρ : x t ´ q t q ¨ δ u mt q dA ` h ż Ω ` p S ‹ n ´ S T ‹ n q ¨ ∇ δu m ´ tr p kS t q δu m ` h ∇ ¨ pp S p q t τ ´ S p q r x p q ˆ ν sq δu m ,s q ´ h ∇ ¨ p S p q r x p q ˆ ν s δu m ,s q` h ∇ ¨ p S p q t ν δu m ,ν ´ S p q r τ ˆ x p q s δα m ‹ q ´ h ∇ ¨ p S p q r τ ˆ x p q s δα m ‹ q ´ h tr p S p q t ∇∇ δu m q´ h p ρ : x t ´ q p q bt q ¨ ∇ δu m ` m t ¨ ∇ δu m ´ h kq ´ t ¨ ∇ δu m ` p ρ : x ´ q q δu m ˘ dA “ h ż B Ω S t ν ¨ δ u mt ds ` h ż B Ω ` p S ‹ n ´ S T ‹ n q ¨ ν ` h p ∇ ¨ S p q t ´ ρ : x p q t ` q p q bt q ¨ ν ´ h p S p q T r ν , ν ˆ x p q sq ,s ´ h p S p q T r ν , ν ˆ x p q sq ,s ` m t ¨ ν ´ h kq ´ t ¨ ν ˘ δu m ds ´ p h ż B Ω S p q T r ν , τ ˆ x p q s δα m ‹ ds ` h ż B Ω S p q T r ν , τ ˆ x p q s δα m ‹ ds q´ ż B Ω q ż h p S T N ´ q q ¨ δ u p s, Z q da, (4.8)where α m ‹ “ arctan p u ,ν {p ` ∇ u mt r ν , ν sqq ˘ pπ . Also, we have used the decomposi-tion ∇ δu m “ δu m ,s τ ` δu m ,ν ν and have transformed the integrals ş B Ω T ,s δu m ds and ş B Ω q M δα m ds into integrals over Ω by Stokes’ theorem. Remark 4.1.
In (4.8), the reason that the O p h k q terms can be dropped is because theyare either relatively O p h q smaller than 2 h ş Ω tr p kS t q δu m dA or relatively O p h q smaller than2 h ş Ω tr p S t ∇ δ u mt q dA (since | hk | ă O p h k q termwill be put into the reminder, which are droppable for the same reasoning. We also pointout that, in order to make the above decomposition of ∇ u m as well as the 2D divergenceof T ,s and M well-defined, the unit vectors τ and ν have to be defined in Ω, which can bedone as follows. The boundary B Ω can be described by an implicit function F p θ α q “ θ α “ θ α , τ can be defined asthe unit tangent vector of the curve F p θ α q “ F p θ α q at the point and ν can then be definedvia the formula ν “ τ ^ n . Note that the variables θ α of τ and ν are changed into θ α in(4.8).Now, we are ready to address the boundary conditions, which should come from the last3D edge term. For the 3D case, the vanishing of this term for any δ u leads to the 3D boundarycondition (2.16) for arbitrary Z , which, obviously, a 2D shell theory cannot satisfy. So, fora 2D shell theory one needs to make some special choices for δ u . Here, the criterion is thatthe lateral force q should generate the virtual work; at the same time for such a choice, theremaining three terms on the right-hand side should give the virtual work done by the external3D force at the edge so that after the vanishing of the last term, (4.8) gives the 2D shell virtualwork principle (that is the main reason that the above calculations are about). According tothis criterion, we choose δ u p s, Z q “ δ u mt ` δu n ` p Z ´ h qp δu m ,s p ν ˆ x p q q ´ δα m p τ ˆ x p q qq` p Z ´ h qp u m ,s p ν ˆ x p q q ´ δα m p τ ˆ x p q qq (4.9)14n B Ω q . Then the vanishing of the last integral of (4.8) leads to ż B Ω q ż h S Tt N ¨ δ u mt da ` ż B Ω q ` ż h S T N ¨ n ? g τ dZ ´ p ż h p Z ´ h q S T N ¨ p ν ˆ x p q q? g τ dZ q ,s ´ p ż h p Z ´ h q S T N ¨ p ν ˆ x p q q? g τ dZ q ,s ˘ δu m ds ´ r ż B Ω q ż h p Z ´ h q S T N ¨ p τ ˆ x p q q δα m da ` ż B Ω q ż h p Z ´ h q S T N ¨ p τ ˆ x p q q δα m da s“ ż B Ω q ż h q t ¨ δ u mt da ` ż B Ω q ` ż h q ? g τ dZ ´ p ż h p Z ´ h q q ¨ p ν ˆ x p q q? g τ dZ q ,s ´ p ż h p Z ´ h q q ¨ p ν ˆ x p q q? g τ dZ q ,s ˘ δu m ds ´ r ż B Ω q ż h p Z ´ h q q ¨ p τ ˆ x p q q δα m da ` ż B Ω q ż h p Z ´ h q q ¨ p τ ˆ x p q q δα m da s (4.10)Next we shall examine each integral on the left-hand side of (4.10) upon using the Taylorexpansions (i.e., specializing to the 2D shell theory) and its counterpart on the right-handside.1. The first integral L on the left-hand side of (4.10) is found to be L “ h ż B Ω q S Tt ν ¨ δ u mt ds ` O p h q , (4.11)which agrees with the first integral on the right-hand side of (4.8) over B Ω q .The applied in-plane force per unit arc length of B Ω q is p q t “ ş h q t ? g τ dZ , so the firstintegral R on the right-hand side of (4.10) can be written as R “ ż B Ω q p ż h q t ? g τ dZ q ¨ δ u mt ds “ ż B Ω q p q t ¨ δ u mt ds, (4.12)which is the virtual work by the applied 3D in-plane force.2. The second integral L on the left-hand side of (4.10) is L “ h ż B Ω q ` p S ‹ n ´ S T ‹ n q ¨ ν ` h p ∇ ¨ S p q t ´ ρ : x p q t ` q p q bt q ¨ ν ´ h p S p q T r ν , ν ˆ x p q sq ,s ´ h p S p q T r ν , ν ˆ x p q sq ,s ` m t ¨ ν ´ h kq ´ t ¨ ν ˘ δu m ds ` O p h , h k q , (4.13)where use has been made of (3.30) and (3.6). We see that L is same as the second integralon the right-hand side of (4.8) over B Ω q .The applied shear force per unit arc length of B Ω q is q s “ ş h q ? g τ dZ . The twistingmoment at the edge about the middle surface due to the applied force q is written as T q “ ż h pp x ´ x p r , h qq ˆ q q ¨ N m ? g τ dZ “ ż h p Z ´ h qp ν ˆ x p q q ¨ q ? g τ dZ ` ż h p Z ´ h qp ν ˆ x p q q ¨ q ? g τ dZ ` O p h , h k q , (4.14)15hose derivative T q,s with respect to the arc length variable is equivalent to a downward shearforce. Then, the second integral R on the right-hand side is R “ ż B Ω q p q s ´ T q,s q δu m ds ` O p h , h k q “ ż B Ω q p q δu m ds ` O p h , h k q , (4.15)where p q is the total effective applied shear force per unit arc length of B Ω q , and one can see R is the virtual work done by the applied 3D force due to the virtual displacement δu m .3. The third term L on the left-hand side of (4.10) is L “ ´ h ż B Ω q S p q T r ν , τ ˆ x p q s δα m ‹ ds ´ h ż B Ω q S p q T r ν , τ ˆ x p q s δα m ‹ ds ` O p h , h k q , (4.16)which is the same as the third term on the right-hand side of (4.8) over B Ω q .The bending moment at the edge point about the middle surface due to the applied force q is p m “ ż h pp x ´ x p r , h q ˆ q q ¨ T m ? g τ dZ “ ż h p Z ´ h qp τ ˆ x p q q ¨ q ? g τ dZ ` ż h p Z ´ h qp τ ˆ x p q q ¨ q ? g τ dZ ` O p h , h k q . (4.17)Then, the third term R on the right-hand side of (4.10) can be written as R “ ´ ż B Ω q p m δα m ‹ ds ` p h , h k q , (4.18)which is the virtual work by the applied 3D force due to the virtual rotation angle.Finally, the equalities L i “ R i p i “ , , q lead to the following boundary conditions onthe traction edge B Ω q :2 h S Tt ν “ p q t , (4.19)2 h ` p S ‹ n ´ S T ‹ n q ¨ ν ` h p ∇ ¨ S p q t ´ ρ : x p q t ` q p q bt q ¨ ν ` h p S p q Tt r ν , τ sq ,s ` m t ¨ ν ´ h kq ´ t ¨ ν ˘ “ p q , (4.20)23 h S p q Tt r ν , ν s ` h S p q T r ν , τ ˆ u p q s ` h S p q T r ν , τ ˆ u p q s “ p m , (4.21)where q t and p q are respectively the applied in-plane force and total effective shear force (perunit arc length of B Ω q ), and p m is the applied bending moment about the middle surface,which are all supposed to be prescribed. In the above equations, we have made use of x p q “ n ` u p q and x p q “ u p q , which result from the relation x “ X ` u , and the two terms ´ h p S p q T r ν , ν ˆ u p q sq ,s and ´ h p S p q T r ν , ν ˆ u p q sq ,s have been dropped in (4.20) for thefollowing reason: for large deformations, they are O p h q smaller than p S ‹ n ´ S T ‹ n q ¨ ν , whilefor small deformations, they are smaller than h p S p q Tt r ν , τ sq ,s . Thus, no matter for largeor small deformations they can be dropped.Based on work conjugates, on the displacement edge B Ω , the boundary conditions are: u mt “ p u mt , u m “ p u m , α m ‹ “ p α m ðñ u m ,ν ` ∇ u mt r ν , ν s “ tan p p α m q , (4.22)16here u m “ u p q ` h u p q ` O p h q , and p u m and p α m are respectively the prescribed displacementand rotation angle of the middle surface.Upon using these boundary conditions for the right-hand side of (4.8), we obtain the 2Dshell virtual work principle (as the right-hand side represents the virtual work done by theapplied effective 3D force at the edge):2 h ż Ω p tr p S t ∇ δ u mt q ` k αβ S β g α ¨ δ u mt ` p ρ : x t ´ q t q ¨ δ u mt q dA ` h ż Ω ` p S ‹ n ´ S T ‹ n q ¨ ∇ δu m ´ tr p kS t q δu m ´ h tr p S p q t ∇∇ δu m q´ h p ρ : x t ´ q p q bt q ¨ ∇ δu m ` m t ¨ ∇ δu m ´ h kq ´ t ¨ ∇ δu m ` p ρ : x ´ q q δu m ˘ dA “ ż B Ω q p q t ¨ δ u mt ds ` ż B Ω q p q δu m ds ´ ż B Ω q p m δα m ‹ ds. (4.23)In obtaining the above equation, the following four terms in (4.8) have been dropped:23 h ∇ ¨ pp S p q t τ ´ S p q r x p q ˆ ν sq δu m ,s q , ´ h ∇ ¨ p S p q r x p q ˆ ν s δu m ,s q (4.24)23 h ∇ ¨ p S p q t ν δu m ,ν ´ S p q r τ ˆ x p q s δα m ‹ q , ´ h ∇ ¨ p S p q r τ ˆ x p q s δα m ‹ q (4.25)which can be justified as follows. From the relations x p q “ n ` u p q and x p q “ u p q , the twoterms in (4.24) can be simplified as ´ h ∇ ¨p S p q u p q q δu m ,s and ´ h ∇ ¨p S p q r u p q ˆ ν s δu m ,s q . From (4.6), the variation of α m ‹ is calculated by δα m ‹ “ p ` ∇ u mt r ν , ν sq δu m ,ν ´ u m ,ν ∇ δ u mt r ν , ν sp ` ∇ u mt r ν , ν sq ` p u m ,ν q “ δu m ,ν ` O p ∇ u δu m ,ν , ∇ u ∇ δ u mt q , (4.26)where the second equality is for small deformations. In (4.25), the terms related to ∇ δ u mt are relatively O p h q smaller than 2 h tr p S t ∇ δ u mt q and can thus be dropped. For the re-maining terms left in (4.25) and the two terms in (4.24), for large deformations, they arerelatively O p h q smaller than 2 h p Sn ´ S T ‹ n q ∇ δu m , while for small deformations they are of O p h S p q u p q δu m ,s , h S p q u p q δu m ,s q and O p h S p q u p q δu m ,ν , h S p q u p q δu m ,ν q , which aresmaller than ´ h tr p S p q t ∇∇ δu m q . Thus they can be dropped no matter the deformation islarge or small.The 2D shell virtual work principle (4.23) supplemented by boundary conditions (4.19)-(4.21) and (4.22) provides a framework for implementing finite element schemes, which willbe left for future investigations. In this section, we apply the previously derived shell theory to study the extension andinflation of an arterial segment, for which the exact solution is available in [20]. We willcompare the asymptotic solution obtained from the shell theory and the exact solution toshow its validity. 17ollowing [1], we consider an artery as a thick-walled circular cylindrical tube, which inits reference configuration has internal and external radii A and B , respectively, and length L . So, its geometry may be described in terms of cylindrical polar coordinates p R, Θ , X q by A ď R ď B, ď Θ ď π, ď X ď L. (5.1)They are related to the Cartesian coordinates p X , X , X q by X “ R cos Θ , X “ R sin Θ , X “ X. (5.2)In the notation of the shell theory, we have the corresponding relations θ “ Θ , θ “ X, Z “ R ´ A, h “ B ´ A. (5.3)We choose the inner surface of the circular cylindrical tube as the base surface. Let p e R , e Θ , e X q denote the standard basis vectors of the cylindrical polar coordinates. A direct calculationusing (5.2) shows g “ A g “ A e Θ , g “ g “ e X , g “ g “ e R “ n . (5.4)Thus the 2D gradient operator is given by ∇ “ A BB Θ e Θ ` BB X e X . The curvature tensor iscalculated by k “ ´ n ,α b g α “ ´ A e Θ b e Θ , which implies that H “ ´ A and K “ p r, θ, z q by a ď r ď b, ď θ ď π, ď z ď l, (5.5)where a, b and l are the deformed counterparts of A, B and L respectively and deformation isgiven by r “ r p R q , θ “ Θ , z “ λ z X, (5.6)where λ z “ l { L is the uniform stretch in the axial direction. Let p e r , e θ , e z q denote thestandard basis vectors of the cylindrical polar coordinates p r, θ, z q which actually agree with p e R , e Θ , e X q . In cylindrical polar coordinates, the shell equations (3.36) and (3.37) take thefollowing form1 A B S Θ θ B Θ ` B S Xθ B X ` A S Θ r “ ρ : x Θ ´ µ p h q q ` Θ ` q ´ Θ h ´ q b Θ , (5.7)1 A B S Θ z B Θ ` B S Xz B X “ ρ : x X ´ µ p h q q ` X ` q ´ X h ´ q bX , (5.8)1 A p B S ‹ Θ r B Θ ´ B S T ‹ Θ r B Θ q ` B S ‹ Xr B X ´ B S T ‹ Xr B X ´ A S Θ θ ` h p A B S p q Θ θ B Θ ` A B S p q Θ z B Θ B X ` A B S p q Xθ B Θ B X ` B S p q Xz B X q“ ρ : x R ´ µ p h q q ` R ` q ´ R h ´ q bR ` h p A BB Θ p ρ : x p q Θ ´ q p q b Θ q ` BB X p ρ : x p q X ´ q p q bX qq´ p A B m Θ B Θ ` B m X B X q ´ hA B q ´ Θ B Θ , (5.9)18here S , S ‹ , S T ‹ , x and q b are defined below (3.37).The deformation gradient arsing from the deformation (5.6) is given by F “ rR e θ b e Θ ` λ z e z b e X ` r e r b e R . (5.10)On the inner and outer surfaces of the circular cylindrical tube, we consider the tractionboundary conditions caused by the internal pressure P q ´ “ P F p q´ T n “ P λ z aA e R , q ` “ . (5.11)On its end surface, we impose a resultant axial force F “ π ż BA S Xz R dR ´ πa P. (5.12)The artery is modelled as an incompressible hyperelastic material reinforced by two sym-metrically disposed families of fibres, which has a strain energy function [21] given by W p I , I , I q “ c p I ´ q ` k k ÿ i “ , p e k p I i ´ q ´ q , (5.13)where I “ tr p C q is the first principal invariant of the right Cauchy-Green tensor C “ F T F ,and I “ M ¨ p CM q and I “ M ¨ p CM q , where the unit vectors M “ cos ϕ e Θ ` sin ϕ e X and M “ ´ cos ϕ e Θ ` sin ϕ e X represent the directions of the two fibres. It follows from(5.10) that I and I are I “ I “ r R cos ϕ ` λ z sin ϕ : “ I. (5.14)For the strain energy function (5.13), the associated nominal stress is given by S “ c F T ` k p I ´ q e k p I ´ q M b F M ` k p I ´ q e k p I ´ q M b F M ´ p F ´ . (5.15)First, substituting (5.10) into (5.15) and doing a Taylor expansion yield S p q “ p c r A ´ p Ar ` k p I ´ q e k p I ´ q r A cos ϕ q e Θ b e θ ` p cλ z ´ p λ z ` k p I ´ q e k p I ´ q λ z sin ϕ q e X b e z ` p cr ´ p r q e R b e r , (5.16) S p q “ p c r A ´ r A ´ p r ´ r Ar ´ p Ar ` k pp ` k p I ´ q q I r A ` p I ´ q r A ´ r A q e k p I ´ q cos ϕ q e Θ b e θ ` p´ p λ z ` k p ` k p I ´ q q e k p I ´ q I λ z sin ϕ q e X b e z ` p r ´ p r ` p r r q e R b e r , (5.17)where r i , p i , I i denote the i th derivatives of r, p, I with respect to Z at Z “
0, respectively; inparticular, we have I “ I | Z “ “ r A cos ϕ ` λ z sin ϕ, I “ B I B Z | Z “ “ r p r A ´ r q A cos ϕ. (5.18)19ext we obtain from (3.3) and (3.7) the recurrence relation for p and r : p “ c A λ z r ` P, r “ Aλ z r , (5.19)and from (3.12) and (3.13) the recurrence relation for p and r : p “ ´ c p λ z r ´ A q λ z Ar ´ k e k p I ´ q p I ´ q cos ϕλ z A , r “ λ z r ´ A λ z r . (5.20)Finally the only nontrivial shell equation (5.9) becomes1 A p S p q Θ θ ` hS p q Θ θ q “ q ´ R h “ P h λ z r A . (5.21)Substituting the recurrence relations (5.19) and (5.20) into the above equation, we obtain anequation involving r only as expected (cid:37) ´ cλ ´ a λ ´ z p λ a λ z ´ q ´ k e k p I ´ q p I ´ q λ z ´ cos ϕ ` h ˚ ` (cid:37)λ ´ a λ ´ z ` c λ ´ a λ ´ z p λ a λ z ´ λ a λ z ` λ a λ z ´ q ` k e k p I ´ q λ ´ a λ ´ z cos ϕ ˆ pp λ a λ z ´ qp I ´ q ` λ a p λ a λ z ´ qp ` k p I ´ q q cos ϕ q ˘ ` h ˚ (cid:37)λ ´ a λ ´ z p λ a λ z ´ q “ , (5.22)where the scales are set as h ˚ “ h { A , P “ (cid:37) h { A and λ a “ r { A “ a { A . We observe from(5.22) that (cid:37) “ cλ ´ a λ ´ z p λ a λ z ´ q ` k p I ´ q e k p I ´ q λ z ´ cos ϕ ` O p h ˚ q . (5.23)Substituting the above equation into the O p h ˚ q term of (5.22), we have P “ (cid:37)h ˚ “ h ˚ p cλ ´ a λ ´ z p λ a λ z ´ q ` k p I ´ q e k p I ´ q λ z ´ cos ϕ q´ h ˚ ` cλ ´ a λ ´ z p λ a λ z ` λ a λ z ´ q ` k e k p I ´ q λ ´ a λ ´ z cos ϕ ˆ p λ a λ z p I ´ q ` λ a p λ a λ z ´ qp ` k p I ´ q q cos ϕ q ˘ ` O p h ˚ q . (5.24)Then according to (4.19), the boundary condition (5.12) gives2 h pp ` hA q S p q Xz ` hS p q Xz q “ F ` πa P πA . (5.25)Substituting (5.24) into above equation, we have F ˚ “ h ˚ p cλ ´ a λ ´ z p λ a λ z ´ λ a λ z ´ q ` k e k p I ´ q λ ´ z p I ´ qp λ z sin ϕ ´ λ a cos ϕ qq` h ˚ ` cλ ´ a λ ´ z p λ a λ z ` λ a λ z ´ λ a λ z ´ λ a λ z ` q` k e k p I ´ q λ ´ z ` p I ´ qpp λ a λ z ´ q cos ϕ ` λ z sin ϕ q` p λ a λ z ´ qp ` k p I ´ q qp λ a cos ϕ ´ λ z sin ϕ cos ϕ q ˘˘ ` O p h ˚ q , (5.26)where F ˚ “ F {p πA q is the normalized resultant axial force. Equations (5.24) and (5.26)form the asymptotic solution of the problem.20n the other hand, the problem has an exact solution of the following form [20]: P “ ż λ a λ b p λ λ z ´ q ´ ψ λ dλ, (5.27) F “ πA p λ a λ z ´ q ż λ a λ b p λ λ z ´ q ´ p λ z ψ λ z ´ λψ λ q λ dλ, (5.28)where λ b “ b { B “ b λ ´ z pp λ a λ z ´ q A { B ` q , ψ λ “ B ψ {B λ , ψ λ z “ B ψ {B λ z , and ψ is givenby ψ p λ, λ z q “ c p λ ` λ z ` λ ´ λ ´ z ´ q ` k k p e k p λ cos ϕ ` λ z sin ϕ ´ q ´ q , (5.29)Doing a routine Taylor expansion, we see that P “ h ˚ λ ´ a λ ´ z ψ λ p λ a , λ z q ´ h ˚ λ ´ a λ ´ z p ψ λ p λ a , λ z q ` λ a p λ a λ z ´ q ψ λλ p λ a , λ z qq ` O p h ˚ q , (5.30) F ˚ “ h ˚ λ ´ z p λ z ψ λ z p λ a , λ z q ´ λ a ψ λ p λ a , λ z qq ` h ˚ λ ´ a λ ´ z ` λ a λ z ψ λ z p λ a , λ z q ´ ψ λ p λ a , λ z q` p λ a λ z ´ qp λ a ψ λλ p λ a , λ z q ´ λ z ψ λλ z p λ a , λ z qq ˘ ` O p h ˚ q , (5.31)where ψ λλ “ B ψ {B λ and ψ λλ z “ B ψ {B λ B λ z . If the expansions are carried out on the middlesurface, then the O p h ˚ q terms are not present, and the errors are of O p h ˚ q as well; seeequations (6.5) and (6.6) in [22]. Using (5.29), it is easy to check that the exact solution(5.30) and (5.31) are the same as the asymptotic solution (5.24) and (5.26), validating theshell equations.To illustrate a numerical example, we set the geometrical and material parameters of theartery as in Table 1; these parameters are cited from [21] and are given for a carotid arteryfrom a rabbit.In Figure 1, we compare the exact solution and the asymptotic solution of the pressure P and the normalized resultant axial force F ˚ for the artery described by the above parameters.It is seen that the asymptotic solution is very close to the exact one, which can be viewed asa numerical validation of the shell equations.Table 1: Geometrical and material data for a carotid artery from a rabbit A (mm) 2 h (mm) c (kPa) k (kPa) k (-) ϕ ρ (g { cm )1 .
43 0 .
26 3 2 . . ˝ . As an application of the derived refined shell theory, we consider the plane-strain vibrationsof an artery superimposed on a pressurized state considered in the previous section. Theresults may be useful in determining the material parameters of an artery. Due to the spacelimit, other vibration modes together with wave propagation will be reported in a separatepaper. The shell equations are three nonlinear PDEs for x p q . For deformations superimposed21 xact solutionasymptotic solution1.00 1.05 1.10 1.15 1.20 1.25 1.30 λ a P ( kPa ) (a) exact solutionasymptotic solution1.0 1.1 1.2 1.3 1.4 1.5 λ z F * ( kPa ) (b)Figure 1: Comparison of the exact solution and the asymptotic solution (a) Variation of theinner pressure P with respect to λ a for fixed λ z “ F ˚ “ F {p πA q with respect to λ z for fixed λ a “ x p q “ x p q b ` δ u p q , where the known vector x p q b is theposition vector of the deformed bottom surface in the base state and δ u p q is the incrementaldisplacement vector. For the pressurized state, we have x p q b “ r e R ` λ z X e X . For theplane-strain vibration modes, we set the components of δ u p q to be δu p q Θ “ U exp p i p n Θ ´ ωt qq , δu p q X “ V exp p i p n Θ ´ ωt qq , δu p q R “ W exp p i p n Θ ´ ωt qq , (6.1)where p U, V, W q are constants, and ω is the angular frequency and n is the circumferentialmode number. Substituting the above two equations into the shell equations in cylindri-cal polar coordinates (5.7)-(5.9) and linearizing, one has three linear algebraic equations for p U, V, W q in the form: ¨˚˝ m m m m m ˛‹‚¨˚˝ UVW ˛‹‚ “ ¨˚˝ ˛‹‚ , (6.2)where the coefficients m , etc. are related to n , ω and the known quantities in the basestate, whose expressions are omitted. For the existence of nontrivial solutions, we need thedeterminant of the coefficient matrix to be zero, which leads to D D “ D “ m and D “ m m ´ m m . We note that this equation gives a relation between thefrequency and the material parameters of an artery; in particular, it may be used to determinethe material parameters of an artery, if the technology is available to measure its vibrationfrequency. The equation D “ δu p q X that is also independent of X , which is thus called the axialmode . The equation D “ X -independent coupled motions with bothcircumferential and radial displacements but without axial displacements, which are calledthe circumferential-radial mode and radial-circumferential mode respectively. This way ofnaming is according to their displacement components when n approaches zero. Precisely,when n “
0, the circumferential-radial mode has the circumferential displacement only andthe radial-circumferential mode has the radial displacement only. Now, we examine the effectsof the axial stretch, pressure and fibre angle on the frequencies for different mode numbers n (with the same material and geometric parameters in the previous section). The numericalresults will be displayed in terms of the non-dimensional frequency ω ˚ : “ ω h { a c { ρ .22e first investigate how the axial pre-stretch affects the frequencies of the plane-strainvibration modes of the pressurized artery. For fixed P “ λ z “ , . , .
6, the frequencies of the plane-strain vibration modes areshown in Table 2. The circumferential-radial mode with n “ n “
0; the circumferential-radial mode with n “ λ z ω ˚ , n “ ω ˚ , n “ ω ˚ , n “
31 0 .
330 0 .
661 0 . . .
399 0 .
799 1 . . .
503 1 .
005 1 .
508 (a) λ z ω ˚ , n “ ω ˚ , n “
31 0 .
354 0 . . .
392 0 . . .
438 0 .
785 (b) λ z ω ˚ , n “ ω ˚ , n “ ω ˚ , n “ ω ˚ , n “
31 0 .
525 0 .
806 1 .
290 1 . . .
511 0 .
799 1 .
288 1 . . .
529 0 .
837 1 .
348 1 .
905 (c)Next we turn to determine the influence of the pressure on the frequencies of the plane-strain vibration modes. For fixed λ z “ P “ , , P ω ˚ , n “ ω ˚ , n “ ω ˚ , n “
30 0 .
242 0 .
484 0 . .
330 0 .
661 0 . .
410 0 .
820 1 .
229 (a)
P ω ˚ , n “ ω ˚ , n “
30 0 .
032 0 . .
354 0 . .
505 1 .
017 (b)
P ω ˚ , n “ ω ˚ , n “ ω ˚ , n “ ω ˚ , n “
30 0 .
463 0 .
644 1 .
006 1 . .
525 0 .
806 1 .
290 1 . .
644 0 .
969 1 .
552 2 .
181 (c)Finally, we check the effect of the fibre angle on the frequencies of the plane-strain vibrationmodes. For fixed λ z “ P “ ϕ “ ˝ , ˝ , ˝ , the frequencies of the plane-strain vibration modes are shown in Table4. It is seen that the frequencies of all vibration modes increase with the mode number. In23ddition, among the three vibration modes, the frequencies of the circumferential-radial modeand radial-circumferential mode decrease with the fibre angle, while frequencies of the axialmode does not always decrease with the fibre angle.Table 4: The frequencies of the plane-strain vibration modes at different fibre angles (a) Axialmode (b) Circumferential-radial mode (c) Radial-circumferential mode ϕ ω ˚ , n “ ω ˚ , n “ ω ˚ , n “ ˝ .
330 0 .
661 0 . ˝ .
364 0 .
728 1 . ˝ .
350 0 .
699 1 .
049 (a) ϕ ω ˚ , n “ ω ˚ , n “ ˝ .
354 0 . ˝ .
344 0 . ˝ .
322 0 .
541 (b) ϕ ω ˚ , n “ ω ˚ , n “ ω ˚ , n “ ω ˚ , n “ ˝ .
525 0 .
806 1 .
290 1 . ˝ .
410 0 .
661 1 .
068 1 . ˝ .
259 0 .
479 0 .
798 1 .
139 (c)
A consistent static finite-strain shell theory is available in the literature (see [3]), whichinvolves three shell constitutive relations (deducible from the 3D constitutive relation) andsix boundary conditions at each edge point. This work first presents a consistent dynamic finite-strain shell theory for incompressible hyperelastic materials in parallel. Novel aspects ofour current study include: 1. The derivation of the refined shell equations through elaboratecalculations which single out the bending effect with only two shell constitutive relations. 2.Many insights can be deduced from the refined shell equations. 3. It is not an easy taskto get the proper number and proper form of physically meaningful boundary conditions ina shell theory. Here, by using the weak form of the shell equations and the variation ofthe 3D Lagrange functional, four shell boundary conditions at each edge point are derived.4. The 2D shell virtual work principle is obtained. A major advantage of this new shelltheory is that its derivation does not involve any ad hoc kinematic or scaling assumptions(as almost all the existing derived shell theories for incompressible hyperelastic materials do).Due to its consistency with the 3D formulation in an asymptotic sense, one does not need toworry about its reliability in predicting the behaviors of incompressible hyperelastic shells forvarious loading conditions. In contrast, for assumptions-based shell theories some defects areevident. For example, some such shell theories involve higher-order stress resultants, whosephysical meanings are not clear, and one does not know how to impose the proper boundaryconditions for them. Another example is the Donnell shell theory, for which the traction fromthe top and bottom surfaces is assumed to be imposed on the middle surface, and if the sheartraction on the top and bottom surfaces has the equal magnitude and opposite sign, that shelltheory does not work. Another simple example is that some shell theories use the assumptionthat the thickness does not change, which is obviously not valid when a large tensile loadis applied at the edge (e.g., large uniform extension of a tube). Due to the simplicity of24ome assumptions-based shell theories, if, for particular applications, experiences/intuitionsindicate that the assumptions involved do not cause a big error, by all means, they can beused. So, at least in theory, there are two differences between the present shell theory andthose assumptions-based ones: prediction reliability (or confidence level) and generality. Thisshell theory is also tested against a benchmark problem: the extension and inflation of anarterial segment. Good agreement with the exact solution to a suitable asymptotic order givesa verification of this shell theory. As an application to a dynamic problem, the plane-strainvibrations in a pressurized artery is considered, and the results reveal the influences of theaxial pre-stretch, pressure and fibre angle on the vibration frequencies, which may be usefulfor determining the artery parameters.Due to the space limit, we only present one application. In subsequent works, we intendto develop a general incremental shell theory by linearizing the present shell theory arounda known base state. Then, we shall study wave propagation in an infinitely-long pressurizedartery and vibrations in all mode types in a finitely-long pressurized artery with suitableedge conditions. Analytical and numerical studies based on this shell theory for determiningsome post-bifurcation behaviors of incompressible hyperelastic shells will be left for futureinvestigations.
References [1] G. A. Holzapfel and R. W. Ogden, “Constitutive modelling of arteries,”
Proceedings ofthe Royal Society A , vol. 466, no. 2118, pp. 1551–1597, 2010.[2] A. E. H. Love, “Xvi. the small free vibrations and deformation of a thin elastic shell,”
Philosophical Transactions of the Royal Society of London.(A.) , no. 179, pp. 491–546,1888.[3] Y. Li, H.-H. Dai, and J. Wang, “On a consistent finite-strain shell theory for incom-pressible hyperelastic materials,”
Mathematics and Mechanics of Solids , vol. 24, no. 5,pp. 1320–1339, 2019.[4] J. Makowski and H. Stumpf, “Finite strains and rotations in shells,” in: PietraszkiewiczW. ed.,
Finite Rotations in Structural Mechanics , pp. 175–194, Springer, Berlin, 1986.[5] M. Itskov, “A generalized orthotropic hyperelastic material model with application toincompressible shells,”
International Journal for Numerical Methods in Engineering ,vol. 50, no. 8, pp. 1777–1799, 2001.[6] D. Chapelle, C. Mardare, and A. M¨unch, “Asymptotic considerations shedding light onincompressible shell models,”
Journal of Elasticity , vol. 76, no. 3, pp. 199–246, 2004.[7] J. Kiendl, M.-C. Hsu, M. C. Wu, and A. Reali, “Isogeometric Kirchhoff–Love shell for-mulations for general hyperelastic materials,”
Computer Methods in Applied Mechanicsand Engineering , vol. 291, pp. 280–303, 2015.258] M. Amabili, I. Breslavsky, and J. Reddy, “Nonlinear higher-order shell theory for incom-pressible biological hyperelastic materials,”
Computer Methods in Applied Mechanics andEngineering , vol. 346, pp. 841–861, 2019.[9] H. Li and M. Chermisi, “The von K´arm´an theory for incompressible elastic shells,”
Calculus of Variations and Partial Differential Equations , vol. 48, no. 1-2, pp. 185–209,2013.[10] H.-H. Dai and Z. Song, “On a consistent finite-strain plate theory based on three-dimensional energy principle,”
Proceedings of the Royal Society A: , vol. 470, no. 2171,p. 20140494, 2014.[11] Z. Song and H.-H. Dai, “On a consistent dynamic finite-strain plate theory and its lin-earization,”
Journal of Elasticity , vol. 125, no. 2, pp. 149–183, 2016.[12] Z. Song and H.-H. Dai, “On a consistent finite-strain shell theory based on 3-d nonlinearelasticity,”
International Journal of Solids and Structures , vol. 97, pp. 137–149, 2016.[13] J. Wang, Z. Song, and H.-H. Dai, “On a consistent finite-strain plate theory for incom-pressible hyperelastic materials,”
International Journal of Solids and Structures , vol. 78,pp. 101–109, 2016.[14] F.-F. Wang, D. J. Steigmann, and H.-H. Dai, “On a uniformly-valid asymptotic platetheory,”
International Journal of Non-Linear Mechanics , vol. 112, pp. 117–125, 2019.[15] P. G. Ciarlet, “An introduction to differential geometry with applications to elasticity,”
Journal of Elasticity , vol. 78, no. 1-3, pp. 1–215, 2005.[16] D. J. Steigmann, “Extension of Koiter’s linear shell theory to materials exhibiting ar-bitrary symmetry,”
International Journal of Engineering Science , vol. 51, pp. 216–232,2012.[17] R. W. Ogden,
Non-linear elastic deformations . Dover, New York, 1997.[18] D. J. Steigmann, “Koiters shell theory from the perspective of three-dimensional nonlin-ear elasticity,”
Journal of Elasticity , vol. 111, no. 1, pp. 91–107, 2013.[19] E. Ventsel and T. Krauthammer,
Thin plates and shells: theory, analysis and applica-tions . Marcel Dekker, New York, 2001.[20] D. Haughton and R. Ogden, “Bifurcation of inflated circular cylinders of elastic materialunder axial loading II. exact theory for thick-walled tubes,”
Journal of the Mechanicsand Physics of Solids , vol. 27, no. 5-6, pp. 489–512, 1979.[21] G. A. Holzapfel, T. C. Gasser, and R. W. Ogden, “A new constitutive framework for ar-terial wall mechanics and a comparative study of material models,”
Journal of Elasticity ,vol. 61, no. 1-3, pp. 1–48, 2000.[22] Y. Fu, J. Liu, and G. Francisco, “Localized bulging in an inflated cylindrical tube ofarbitrary thickness–the effect of bending stiffness,”