The elastic theory of shells using geometric algebra

We present a novel derivation of the elastic theory of shells. We use the language of geometric algebra, which allows us to express the fundamental laws in component-free form, thus aiding physical interpretation. It also provides the tools to express equations in an arbitrary coordinate system, which enhances their usefulness. The role of moments and angular velocity, and the apparent use by previous authors of an unphysical angular velocity, has been clarified through the use of a bivector representation. In the linearized theory, clarification of previous coordinate conventions which have been the cause of confusion is provided, and the introduction of prior strain into the linearized theory of shells is made possible.


Introduction
Thin shells have been an active subject of research for some considerable time, however, in attempting to understand the selfexcited oscillations of flexible tubes, we have had difficulties finding a complete and rational theory in which the underlying physical principles are clear, and which is easy to apply to the practical problem at hand. Specifically, it was found that in order to have a full understanding of the assumptions of various shell theories, it was necessary to derive our own from first principles. We found that in doing this we were able to produce a theory with improved clarity, brevity and with explicit results for linearization about a deformed state. For reference, our complete nomenclature is given in table 1.
We also have an interest in applying geometric algebra (GA) [1][2][3] to new areas of the physical sciences. GA provides the tools to formulate physical laws with as little reference to coordinate systems as possible, which helps with the first aim of clarifying the physical meaning of the equations produced, but it also provides simple tools to allow these equations to be represented in arbitrary coordinate systems, which ensures practical use. This article aims to provide, to our knowledge for the first time in this area, an introduction to shell theory using GA. While in this article we restrict the introduction of GA to the use of bivectors to represent torques and angular velocities, we hope that this will pave the way for more radical developments, such as those completed for the theory of rods [4].              There are a large number (at least 10) of linearized shell theories [5]. The derivations of these theories use a wide variety of notations, coordinate systems and conventions, making it very difficult to compare the assumptions made. In addition, none of the theories reviewed by Leissa [5] allow for prior strain of the shell, which we wish to include for our own analysis. More general shell theories have also been produced, with the most extensive probably that by Naghdi [6], which provides the basis for more modern studies such as [7][8][9], though the theory of Koiter [10] has also been popular with some authors. While rigorous, these theories have limited practical use. They generally require the use of differential geometry [7,11] whose indicial expressions often hide much of the physical meaning of the equations. The general theory presented by Antman [8] is relegated to a final chapter that does not stand alone, meaning that the entire book must be read to use the shell theory. Naghdi [6] discusses in detail the different advantages of developing a shell theory directly from three-dimensional elasticity or by considering two-dimensional surfaces from the start. Antman [8] restricts his development to the former, but we feel that a more concise and lucid theory can be obtained from the latter.
We aim to use GA to develop an accessible, concise, rational shell theory that can be easily linearized to include pre-strain. In doing this, we will provide new developments in the representation of moments and angular velocities with bivectors, and in the representation of bending, which is where most disagreements occur in linearized shell theories.

Geometry of surfaces
Let B and S be two-dimensional surfaces embedded in three-dimensional Euclidean space E 3 . B is the reference configuration of the surface, and S is the spatial configuration, and the two are related by the motion φ t . At time t, the point X ∈ B is moved to φ t (X) ∈ S. Let {X i } be coordinates over B, and {x i } be coordinates over S. We follow the convention that the indices i, j, k, . . . run over 1, 2, and the indices a, b, c, . . . run over 1, 2, 3. We restrict {x i } to be convected coordinates such that x i (x) = X i (φ −1 t (x)), where x ∈ S. We denote the frame associated with {X i } by E i = ∂X/∂X i , and similarly, e i = ∂x/∂x i . The reciprocal frames are denoted by {E i } and {e i }, and are defined to satisfy E i · E j = e i · e j = δ i j . The frames on each configuration are illustrated in figure 1.  Figure 1. Surface geometry.
The local pseudoscalars in the reference and spatial configurations are I = E 1 ∧ E 2 /|E 1 ∧ E 2 | and i = e 1 ∧ e 2 /|e 1 ∧ e 2 |, which satisfy I 2 = i 2 = −1. We denote the pseudoscalar of E 3 by I 3 . We have defined orientations of both configurations and E 3 with these pseudoscalars, which allows us to define the normal vectors to the surfaces, E 3 = −I 3 I and e 3 = −I 3 i. E 3 and e 3 are unit vectors perpendicular to the other frame vectors, so E 3 = E 3 and e 3 = e 3 . {E a } and {e a } now both form a basis of E 3 . The (scalar) volume forms V and v are defined to satisfy VI = E 1 ∧ E 2 and vi = e 1 ∧ e 2 .
The vector derivative of E 3 is denoted by ∇, and the projection of this derivative operator onto either B or S is denoted by ∂. ∂ can be written locally on B as ∂ = E i (∂/∂X i ), and on S as ∂ = e i (∂/∂x i ). For convenience, we define the notation ∂ i = (∂/∂X i ) and ∂ i = (∂/∂x i ). It will be clear from context whether differentiation is on the reference or spatial configuration.
Let G(Y) = Y and g(y) = y be identity functions, where Y is a vector on B, and y is a vector on S. The reason we distinguish these apparently identical linear functions is that they are the metrics of the two surfaces, also called the first fundamental forms. In component form, we have g ab = e a · g(e b ) = e a · e b and g ab = e a · e b . The properties of the reciprocal frame imply that g a b = g b a = δ a b . Analogous results hold for G. The determinant of a function is defined in a coordinate free way by g(i) = (det g)i, from which it is clear that det G = det g = 1. However, it is common to define det(g ij ) = g 11 g 22 − g 12 g 21 , which is not equal to 1, and in fact encodes important geometric information about the manifold. This is possible because the coordinate free definition of det g corresponds to det(g i j ), and not det(g ij ). In fact, we can show that det(g ij ) = v 2 . Recalling the definition of v, this demonstrates in a very obvious way that det(g ij ) is a measure of the 'volume' spanned by the parallelepiped formed from the basis vectors. GA in this instance provides clarification over the fact that g is simply the identity function, and provides a definition of v = det(g ij ) that makes its geometric significance immediately obvious.
We denote the second fundamental forms on B and S by B and b. b is defined to satisfy b(y) = −∂(y ·ė 3 ). In component form, we have b ij = −e j · ∂e 3 /∂x i = e 3 · ∂e j /∂x i , which follows from the fact that e j · e 3 = 0 ⇒ ∂ i (e j · e 3 ) = 0. From this, it is clear that b ij = b ji , and hence b is symmetric, i.e. b(y) =b(y) = −y · ∂e 3 . The eigenvalues of b are the principal curvatures of the surface, denoted by c 1 and c 2 . Analogous results hold for B, whose eigenvalues are denoted by C 1 and C 2 . We define the Christoffel coefficients γ a ib = e a · ∂e b /∂x i = −e b · ∂e a /∂x i , which follows from the fact that e a · e b = δ a b ⇒ ∂ i (e a · e b ) = 0. γ i jk are the usual coefficients associated with a frame on a manifold. The remaining coefficients are closely related to the second fundamental form by γ 3 ij = b ij , γ i j3 = −b i j and γ 3 i3 = 0, as e 3 is a unit vector. We similarly define Γ a ib = E a · ∂E b /∂X i . When considering angular momentum, we will make use of bivectors. To do this, we first introduce some notation. The space of all bivectors in  A, B, C, . . . run over (1, 3), (2, 3), (1, 2). Hence, the space of bivectors is spanned by {e A } and {e A }. Defined in this way, these basis bivectors satisfy e A · e B = δ A B . In an analogous way to vectors, the general bivector ω can be written in component form as ω A = ω · e A (here we follow the conventions of [1, eqn. 1-3.18]). We can also define the bivector Christoffel coefficients γ A iB = e A · ∂e B /∂x i . Given the surface, we have already defined, for which e 3 = e 3 , these satisfy γ (1,2) i (1,2) and γ (1,2) i (2,3) Let m be a bivector valued function of a vector. m(y) can be written as m(y) = m Aa y a e A , where y a = y · e a and m Aa = e A · m(e a ), for example, m (1,2)2 = (e 2 ∧ e 1 ) · m(e 2 ).

Kinematics
We define X(η) to be a path over B parametrized by the scalar η. dX/dη is then a tangent vector to B, and we can also obtain a tangent vector to S, ∂φ t (X)/∂η. The map between these tangent vectors is denoted by F, and is called the deformation gradient. This encodes stretching information for the surface, but also rigid body rotations. Rigid body rotations are not expected to influence constitutive theory, so we construct the Cauchy-Green tensor C(Y) =FF(Y), which is symmetric. We restrict ourselves to deformations that have an inverse and leave the orientation of B unchanged, which means that the eigenvalues of C will be real and positive. It is therefore meaningful to define λ i as the square roots of the eigenvalues of C. These are the principal stretches of the surface. Using the Cauchy-Green tensor we construct the Green-Lagrange strain tensor, , that is only non-zero when the material is locally stretched. Given that {x i } are convected coordinates, e i = F(E i ). This allows us to obtain the component expressions Hence, we see that using convected coordinates, the metric can be used to encode stretching information. However, our definition is coordinate free.
In three-dimensional elasticity, the strain tensor is sufficient to characterize linear constitutive theory. When dealing with shells we must also consider the bending of the shell, or more precisely, the change of curvature from the reference to the spatial configuration. Hence, we define the change of curvature tensor H(Y) =FbF(Y) − B(Y). Using convected coordinates, we obtain the component expression We are also interested in the strain rate, and to this end we consider the rate of change of a tangent vector as it is convected with the surface, ∂F(Y)/∂t = ∂ 2 φ t (X)/∂t∂η = ∂/∂η(∂φ t (X)/∂t) = F(Y) · ∂v = Y · ∂V, where v and V are the velocities referred to the spatial and reference configurations, respectively (figure 1). Using the fact that e 3 is always normal to {e i }, we can write ∂e 3 /∂t = −e 3 · (∂e i /∂t)e i = −e 3 · (e i · ∂v)e i = −v 3|i e i , where v a|i is defined by v a|i = e a · (e i · ∂v). Combining these we can now construct a function that returns the rate of change of a vector, that need not be tangential to S, as it is convected with the motion φ t . We denote this function l(y) = ∂y/∂t = y · ∂v + y · e 3 (∂e 3 /∂t) (note that y need not be tangential to S in this expression). It is useful to decompose this into its symmetric and antisymmetric parts n(y) = 1 2 (l(y) +l(y)) and w(y) = 1 2 (l(y) −l(y)). After some manipulation, the components of n and w are given, in terms of convected coordinates, by n ij = 1 The symmetric tensor n(y) is closely related to E. We define the rate of change of the strain tensor E with time byĖ(Y) = ∂E(Y)/∂t. The components of this tensor are given byĖ ij = n ij , and hence we see thatĖ(Y) = FnF(Y). This will be important in constitutive theory.
The rate of change of the change of curvature tensorḢ(Y) = ∂H(Y)/∂t can be expressed in component form asḢ ij = ∂l 3j /∂x i − γ k ij l 3k − γ a i3 l aj = l 3j|i = e 3 · (e i ·∂l(e j )) (for details see appendix A). We see from this that the rate of change of H with time is the e 3 component of the second spatial derivative of velocity. Note that bothĖ andḢ are symmetric.
w is an antisymmetric function mapping vectors on S into vectors in E 3 . Hence, it has a single characteristic eigenbivector ω such that w(y) = y · ω, which we can extract as ω = 1 2 e a ∧ w(e a ) [1, §3.4]. Defined in this way ω is the local angular velocity of the shell material, represented as a bivector. The vector representation of angular velocity is given by −I 3 ω. If we consider e 3 being convected with a material point on the surface, then the fact that it is defined to be a unit vector allows us to use ω to write ∂e 3 /∂t = e 3 · ω. We need a representation of angular velocity in shell theory because it is not possible to assume, as it is in three-dimensional elasticity, that couple-stresses are negligible. The bivector representation of angular velocity allows for a much more physical representation of the governing laws of shells than that suggested by Naghdi [6], who requires the use of a rotated angular velocity with components normal to the shell removed.

Stress
We consider an arbitrary region of the shell defined by U ⊂ B, which under the motion φ t moves to φ t (U) ⊂ S. In continuum mechanics, it is standard to assume that all the forces on the region U can be described by either body forces or boundary forces, to which we must add body and boundary moments in shell theory. Body forces are expressed in terms of the body force per unit mass b(x, t). The force acting on the region U owing to body forces is given by where ρ is the mass per unit area of the shell, b(X, t) = b(φ t (X), t) and dx, dX are directed volume elements on the spatial and reference configurations. Directed integration theory is introduced in [3, §6.4]. In shell theory, we must also consider body moments. We define the body moment per unit mass c such that the moment acting on the region U owing to body moments is given by Next, we consider boundary forces and moments. We denote a small portion of the boundary ∂φ t (U) by s, with normal vector n. We assume that the material on the outside of s exerts a force f , and moment m, on the material inside. The stress principle of Euler and Cauchy, adapted for a shell, states, as the length s tends to zero, the ratios f / s and m/ s tend to definite limits. Moreover, if two paths passing through a point x have the same normal n, then f / s and m/ s tend to the same value for both of these paths [12].
Using arguments outlined by, among others, [12], we can show that the limits described in this principle can be expressed as linear functions of the normal vector n at each point x ∈ S, given a particular time. This allows us to define the Cauchy stress tensor σ (n) and the couple-stress tensor m(n). We can then write the force on a portion of the shell owing to boundary forces, and the moment on a portion of the shell owing to couple-stresses as where ds is a directed boundary element, related to the normal vector by n|ds| = dsi −1 . σ and m are both tensors on the spatial configuration. We wish to express balance laws on the reference configuration, so we construct the first Piola-Kirchhoff stress tensor T, and the first reference couple-stress tensor M, by T(N) = det F σF −1 (N) and M(N) = det F mF −1 (N). Using these we can write The domain of σ is vectors tangential to S, but its range is E 3 , and similarly the domains of T and S are vectors tangential to B, while their ranges are E 3 . m is not vector valued, but bivector valued, as it represents a moment. Its domain is vectors tangential to S, and its range is the space of all bivectors in E 3 . However, we know that the moments represented by m are because of the stress distribution through the thickness of the shell, and this means that its range is more restricted. More precisely, we can say that m (1,2)i = (e 2 ∧ e 1 ) · m(e i ) = 0. Similarly, we assume that c is due only to shear stresses acting on the upper and lower surfaces of the shell, meaning that its e 1 ∧ e 2 component is zero. Given the use of convected coordinates, the following coordinate expressions hold: N is the tangent space of B. The symmetry of the three-dimensional Cauchy stress tensor implies that N is symmetric in the plane tangent to B. These are convenient when expressing conservation of angular momentum and constitutive laws. Note that M and N are not the physical vector representations of the torque. Their natural emergence in conservation of angular momentum and energy (see §5) explains why [6] was able to make use of an apparently unphysical rotated angular velocity in his formulation, and also justifies the rather strange definition of the vector components of the couple stress given in [5, §1.6.1,eq.1.113]. M and N satisfy the following coordinate expressions: and

Balance laws
We write each balance law as an integral equation expressed on the spatial configuration, and a local equation of motion expressed on the reference configuration. We are able to express all of these in component free form, which is a common advantage of using GA.
Using this we can define the time independent density ρ 0 = ρ det F.
The algebraic manipulations necessary to achieve this expression are given in appendix C. We can split this expression into its e 1 ∧ e 3 , e 2 ∧ e 3 and e 1 ∧ e 2 components. These components involve taking the divergence of a bivector valued function, which is outlined by Hestenes & Sobczyk [1]. However, using the modified first couple-stress tensor, these components can be written in a more familiar form: where, for convenience, we have defined c 1 = c · (e 3 ∧ e 1 ) = c (1,3) and c 2 = c · (e 3 ∧ e 2 ) = c (2,3) . The bivector versions of these expressions are equally valid, and easier to interpret physically, but less familiar because they involve bivector components. To obtain more familiar expressions, we need to use modified tensors such as M, whose physical meaning is less immediately obvious. The use of bivectors to represent angular velocities and torques has illuminated why it was necessary for [6] to use apparently unphysical quantities to develop his shell theory. Conservation of angular momentum has two major implications. The first, from the e i ∧ e 3 components of the expression, is that stresses normal to the tangent plane of the surface are determined if the couple-stress and body moment are known. This means that we do not need a constitutive law for these components of the stress, we only need constitutive laws for the components of stress within the plane of the shell, and for the couple stress. The second implication, from the e 1 ∧ e 2 component, is that the modified second Piola-Kirchhoff stress is symmetric in the tangent space of the reference configuration. This is important when considering conservation of energy and in constitutive theory. Energy In this article, we assume isothermal elasticity. It is uncomplicated to include thermal effects, simply requiring the inclusion of the second law of thermodynamics and additional constitutive laws, but this addition does not contribute to our aim here of introducing GA to shell theory for the first time, so is not included. Conservation of energy is therefore given by where e(x, t) is the internal energy per unit mass. The negative signs before the moment terms is consistent with the use of bivectors to represent moments and angular velocities (see appendix B). After some algebraic manipulation (see appendix D), on the reference configuration we obtain where E(X, t) = e(φ t (X), t) (which is not the same as E, the Green-Lagrange strain tensor). Note the appearance of the modified second Piola-Kirchhoff stress, and the modified second couple stress tensor. We know the first of these is symmetric in the tangent space of B from conservation of angular momentum, and the second must be assumed symmetric in order to obtain a determinate theory (this assumption was first proposed in [6, §15]). This allows us to use this expression to derive the constitutive laws given in §6.

Constitutive theory
Our basic constitutive assumption is that E is a function of the local values of the tensors E and H. Applying the chain rule to equation (5.9), and noting that the equation is valid for arbitrary deformations, we obtain the constitutive relations: For an introduction to tensor derivatives, see [3, §11.1.2]. Koiter [10] proposes the following form for ρ 0 E, which can be regarded as the application of the Saint Venant-Kirchhoff material to shells: where E y is Young's modulus, ν is Poisson's ratio and h is the thickness of the shell. From this we obtain the following relationships: Note that this only provides the part of S that is tangential to B. The non-tangential part (S 3i ) is found using conservation of angular momentum. There is a fundamental contradiction in arriving at the results presented here. To arrive at the presented form of ρ 0 E shown, we must make the following assumptions: -the midsurface in the reference configuration remains the midsurface under the motion; -a material line that is normal to the midsurface in the reference configuration remains normal to the midsurface; -the shell thickness (measured normal to the midsurface) is constant over the surface and does not change with time; -the first and second moments of the density relative to the midsurface are zero; -the shell thickness is small compared with its principal radii of curvature; -strains within the shell are small; and -normal stress in the shell is negligible.
When applied to Hooke's law in three dimensions, these assumptions imply that the e 3 component of σ is zero, but we know that in shell theory these components are required for conservation of angular momentum. This basic contradiction remains unresolved.

Linearization
We define the displacement U(X, t) = φ t (X) − X, and we assume that it takes the form U = U 0 + U , where is small. Neglecting terms of O( 2 ) we obtain The Green-Lagrange strain tensor can then be written as We also need to write the change of curvature tensor in its perturbed form. To do this, we first express the convected basis vectors {e a } as where × is the vector cross product, defined by a × b = −I 3 a ∧ b. This allows us to writeFbF as which in turn allows us to express H as We can now write the modified second reference couple stress tensor N as To express the second Piola-Kirchhoff stress tensor, we first need to express F −1 bF(Y): This allows us to write S and T as Conservation of mass can be expressed as ∂ ∂t (ρ det F) = ∂ ∂t (ρ(det F 0 + det F 0 tr(F −1 0 F ))) = 0. (7.9) We define ρ 0 = ρ det F 0 and ρ = ρ det F 0 tr(F −1 0 F ) (note the adjustment of the definition of ρ 0 ). We assume that the initial displacement U 0 satisfies the governing equations separately, so both ρ 0 and ρ are independent of time. Conservation of momentum can be expressed as (ρ 0 + ρ ) ∂ 2 ∂t 2 (U 0 + U ) =Ṫ 0 (∂) + Ṫ (∂) + (ρ 0 + ρ )b. (7.10) We denote the body force acting on the body in its initial deformed state (defined by U 0 ) by b 0 , and then decompose b as b = b 0 + b . This includes the assumption that the additional body force acting on the body after the perturbation U is small. Subtracting conservation of momentum for the initial deformation U 0 , we obtain Usually we assume that U 0 is time independent, meaning that we obtain We can write M as = M 0 (y) + M (y). (7.13) If, as with b, we assume that c can be decomposed as c = c 0 + c , then this allows us to write the perturbed part of conservation of angular momentum as F 0 (E i ) ∧ T (E i ) + F (E i ) ∧ T 0 (E i ) +Ṁ (∂) + ρ 0 c + ρ c 0 = 0. (7.14)

Small displacements
If we assume U 0 = 0 (or that it is constant), then we obtain the following simplifications: and Following the method of Ciarlet [7], and using the coordinate independent notation of GA, the change of curvature tensor takes the form H(Y) = E j ((E j ·∂Ḟ (Y)) · E 3 ). (7.16) From this, it is clear how the linearized change of curvature tensor is closely related to the E 3 components of the second derivative of the displacement field U . It is worth pointing out an anomaly at this point, which is made clearer by the use of GA. Much of the previous work done on linearized shell theory (summarized by Leissa [5], with one of the more rigorous derivations provided by Vlasov [13,14]) uses a rather strange coordinate system, which does not aid comprehension. Coordinates are chosen such that the lines X i = constant define lines of principal curvature on the reference configuration. This allows the components of several tensors to be expressed more simply, as the basis vectors are orthogonal, and are eigenvectors of B. However, the coordinate system is not constrained to be orthonormal, meaning that the reciprocal frame and frame do not coincide (though E i is parallel to E i ). This is not a problem in GA, as we can use an arbitrary coordinate system, but the solution adopted by many authors is to create a new normalized frame {Ê i = E i /|E i |}. Differentiation is performed with respect to the coordinates {X i }, but tensor and vector components are expressed relative to the frame {Ê i }. This adds considerable complication to the expressions for strain and change of curvature, which, through the use of GA, we have simplified.

Uniaxial strain of a cylinder
We now consider the case in which {E i } forms an orthonormal basis, and are also the eigenvectors of B. In this case, we have Γ i jk = 0 and Moreover, we take C 1 = 0 and C 2 = C. We take the background deformation to be uniaxial strain such that U 0 = εX 1 E 1 . X 1 is the axial distance along the cylindrical shell, and X 2 is the azimuthal distance around the circumference. We can write F 0 as F 0 (Y) = Y + ε(Y · E 1 )E 1 , but E 1 is a basis vector specific to the tangent space of the reference configuration. For clarity, we therefore define a unit vector aligned with the axis of symmetry of the cylindrical shell e, which is defined everywhere in E 3 . On the reference configuration e = E 1 . Using this we write: F 0 (Y) = Y + ε(Y · e)e, F 0 (y) = y + ε(y · e)e, f −1 0 (y) = y − ε λ (y · e)e, and det F 0 = 1 + ε = λ, ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (7.21)