A computational framework for polyconvex large strain elasticity for geometrically exact beam theory

In this paper, a new computational framework is presented for the analysis of nonlinear beam finite elements subjected to large strains. Specifically, the methodology recently introduced in Bonet et al. (Comput Methods Appl Mech Eng 283:1061–1094, 2015) in the context of three dimensional polyconvex elasticity is extended to the geometrically exact beam model of Simo (Comput Methods Appl Mech Eng 49:55–70, 1985), the starting point of so many other finite element beam type formulations. This new variational framework can be viewed as a continuum degenerate formulation which, moreover, is enhanced by three key novelties. First, in order to facilitate the implementation of the sophisticated polyconvex constitutive laws particularly associated with beams undergoing large strains, a novel tensor cross product algebra by Bonet et al. (Comput Methods Appl Mech Eng 283:1061–1094, 2015) is adopted, leading to an elegant and physically meaningful representation of an otherwise complex computational framework. Second, the paper shows how the novel algebra facilitates the re-expression of any invariant of the deformation gradient, its cofactor and its determinant in terms of the classical beam strain measures. The latter being very useful whenever a classical beam implementation is preferred. This is particularised for the case of a Mooney–Rivlin model although the technique can be straightforwardly generalised to other more complex isotropic and anisotropic polyconvex models. Third, the connection between the two most accepted restrictions for the definition of constitutive models in three dimensional elasticity and beams is shown, bridging the gap between the continuum and its degenerate beam description. This is carried out via a novel insightful representation of the tangent operator.


Introduction
Most classical beam theories [1][2][3][4][5][6][7][8] are based on the definition of so called beam strain measures, namely the axial-shear and torsional-bending strain vectors. Work conjugates to those variables (typically denoted as resultant contact forces and resultant contact couples [9]) can be derived from a hyperelastic energy functional which is defined in terms of the beam strain measures [1][2][3][4][5][6][7][8]. Other authors prefer an alternative continuum degenerate approach [10,11] where, typically, the Lagrangian strain tensor is retained as the main strain measure and particularised for the kinematic description of the beam. In these latter formulations, the second Piola-Kirchhoff stress tensor emerges as the work conjugate variable and can be derived from a hyperelastic energy functional expressed as a function of the Lagrangian strain tensor.
The present manuscript aims to extend the variational and computational framework recently introduced in Bonet et al. [12] in the context of three dimensional elasticity to the geometrically exact Simo [13] beam model, the starting point of so many other finite element beam type formulations [1][2][3][4][5][6][7]. The variational approach proposed herein can therefore be viewed as a continuum degenerate formulation emerging from that presented in [12].
The formulation by Bonet et al. [12] has been particularly formulated for use in large strain scenarios [14,15], where appropriate convexity criteria [16][17][18] are necessary to ensure the well posedness of the problem. The most wellestablished of these criteria is the concept of polyconvexity [19][20][21][22][23][24][25] whereby the strain energy is expressed as a convex multi-value function of an extended kinematic set defined by the deformation gradient (fibre map), its determinant (volume map) and its cofactor (area map).
Numerous authors have previously incorporated the concept of polyconvexity into computational models for both isotropic and non-isotropic materials for a variety of applications [26][27][28][29][30][31]. However, the classical approach consists of ensuring that the strain energy satisfies the polyconvexity condition first but then proceeds towards a computational solution by re-expressing the energy in terms of the deformation gradient alone. Pioneered in [29] and [12,32,33], an alternative framework is proposed based on maintaining as independent variables the extended kinematic set on which the strain energy is expressed as a convex function, namely, the deformation gradient, its cofactor and its determinant.
In this paper, the latter approach is particularised for a degenerate beam description. The paper also presents the work conjugates of the extended kinematic set in the context of beam theory as well as novel expressions for the tangent operators. The present formulation utilises a novel algebra based on a tensor cross product operation pioneered in [34] and reintroduced and exploited for the first time in [12,32,33] in the context of solid mechanics. The new tensor cross product operation is particularly helpful when dealing with polyconvex constitutive laws, where invariants of the cofactor and the determinant of the deformation feature heavily in the representation of the strain energy functional.
In addition, the paper shows how the novel algebra facilitates the re-expression of any invariant of the deformation gradient, its cofactor and its determinant in terms of the classical beam strain measures. The latter being very useful whenever a classical beam implementation is preferred. As an example, and to the best of the authors knowledge, this is the first time that the strain energy function for a polyconvex Mooney-Rivlin (or Neo-Hookean) material is expressed entirely in terms of the classical beam strain measures.
Convexity of a beam constitutive model with respect to the classical beam strain measures is a well accepted condition [35] that must be satisfied by admissible strain energy functionals in the large strain regime. The present manuscript shows the relationship between a polyconvex constitutive model defined in terms of the deformation gradient, the cofactor and the determinant of the mapping in the continuum and its degenerate beam counterpart defined in terms of the beam strain measures. Hence, the connection between the two most accepted restrictions for the definition of constitutive models in three dimensional elasticity and beams is shown, bridg-ing the gap between the continuum and its degenerate beam description.
The paper is organised as follows. Section 2 briefly revises the computational framework developed in [12] for large strain scenarios and especially tailored for polyconvex constitutive models. Section 3 presents the new degenerate beam formulation which extends the variational framework presented in [12] to the case of geometrically exact beam theory where the deformation gradient, its cofactor and its determinant are retained as the main strain variables. Section 4 presents the classical framework for the geometrically exact beam theory where the physically meaningful axial-shear and torsional-bending strain vectors are used as the main strain measures. Additionally, a link is established between the convexity criteria required for the constitutive model in both continuum and beam descriptions. Section 5 briefly presents the variational principle associated with the proposed continuum degenerate formulation. Section 6 presents the Finite Element discretisation, where the use of the novel tensor cross product algebra leads to alternative tangent operator representations. Section 7 includes representative numerical examples in large strain scenarios, including a comparison against a recently published mixed continuum based formulation [12]. Finally, Sect. 8 provides some concluding remarks and a summary of the key contributions of this paper.

Continuum mechanics 2.1 Continuum kinematics
Consider the three dimensional deformation of an elastic body from its initial configuration occupying a volume V , of boundary ∂ V , into a final configuration at volume v, of boundary ∂v, where x represents the current position of a particle originally at X (see Fig. 1). Virtual and incremental variations of x will be denoted by δu and u, respectively. It will be assumed that x, δu and u satisfy appropriate essential (displacement) boundary conditions on ∂ u V . Additionally, the body is under the action of certain body forces per unit undeformed volume f 0 and traction per unit undeformed The deformation gradient tensor F (or fibre map) and its determinant J (or volume map) are defined as [16] where ∇ 0 denotes the gradient with respect to material coordinates and dv and dV represent elemental volumes in the initial and final configurations, respectively. The elemental area vector is mapped from initial d A to final d a configurations by means of the cofactor tensor H (or area map), which Fig. 1 Deformation mapping of a continuum and associated kinematics variables: F, H, J is related to the deformation gradient F and its determinant J via the Nanson's rule [16] as In references [12,32], the authors employ an alternative definition of the cofactor H, namely H = 1 2 F × F, which simplifies considerably the algebra [33]. The new definition of the area map H is based on a tensor cross product × introduced in [34] and applied within the context of nonlinear elasticity for the first time in [12]. The properties of this tensor cross product developed in [12] have been included in Appendix 1 for completeness.
Crucially, the first and second directional derivatives of H with respect to geometry changes are now easily evaluated as 1 Similarly, the directional derivatives of the volume map J are easily expressed with the help of (139) and (142) as

Polyconvex elasticity
Polyconvexity is now well accepted as a fundamental mathematical requirement that must be satisfied by admissible strain energy functionals used to describe elastic materials in the large strain regime [18][19][20][21][22]24,25,36,37]. The strain 1 The first and second directional derivatives of the deformation gradient tensor F are well known to be D F[δu] = ∇ 0 δu and D 2 F[δu; u] = 0 [16]. energy per unit undeformed volume must be a function of the deformation gradient ∇ 0 x via a convex multi-valued function W as Crucially, frame invariance (objectivity) implies that W must be independent of the rotational components of F and H, which is typically achieved by ensuring that W depends on F and H via the symmetric tensors F T F and H T H, respectively. The three strain measures F, H, and J have work conjugate stresses F , H , and J defined by [12,32] It is then possible [12,32,33] A tangent elasticity operator D 2 [δu; u] is usually needed [16] in order to ensure quadratic convergence of a Newton-Raphson type of solution process, derived in [12,32,33] as follows where the Hessian operator [H W ] denotes the symmetric positive semi-definite operator containing the second derivatives of W (F, H, J ) as Note that the first term on the right hand side of (8) is necessarily positive for δu = u and, therefore, buckling can only be induced by the second "initial stress" term ( H + J F) : (∇ 0 δu × ∇ 0 u). Finally, it is also possible to re-write the tangent elasticity operator (8) as where use of the tensor cross product operations (133) and (134) yields

A polyconvex constitutive model: Mooney-Rivlin material
The well known Mooney-Rivlin constitutive model is polyconvex. An elegant representation of its strain energy functional in terms of the extended kinematic set {F, H, J } [12,32,33] is with α and β suitable non-negative material parameters and f (J ) a convex function of J (refer to remark 1). The conjugate stresses { F , H , J } (6) for this model are obtained as and the first Piola-Kirchhoff (7) results in The Hessian operator [H W ] (9) adopts the following simple expression where I i I j J = δ i j δ I J and H eq W (11) is expressed as The clear benefits of employing the new tensor cross product algebra in the context of polyconvex elasticity have been presented in [12,32] and thoroughly detailed in [33].

Remark 1
An alternative function f (J ) to that used in (13) is Following [12], the material parameters α, β and λ used in (12) and (18) can be related to the classical Lamé parameters μ lin and λ lin in linearised elasticity (in the reference configuration) through The Poisson ratio ν in the reference configuration is related to material parameters λ lin and μ lin as Substitution of Eq. (19) into (20) leads to the following expression for the Poisson ratio in terms of the material parameters α, β and λ Classical beam theories do not consider dilatations or contractions of the beam section. Hence, a Poisson ratio of value ν = 0 is consistent with such kinematical description of the beam. However, from Eq. (21), it can be inferred that it is not possible to find any combination of non-negative (and hence, satisfying polyconvexity) values for the material parameters α, β and λ compatible with ν = 0. On the contrary, if the volumetric function f (J ) in (13) is used instead, the relation between the material parameters α, β and λ relate to μ lin and λ lin as which leads to an expression for the Poison ratio ν lin The definition of f (J ) in (13) allows to model the particular case ν lin = 0 without violating polyconvexity (non-negative values for α and β with λ = 0).

A non-polyconvex constitutive model: Saint Venant-Kirchhoff material
A popular constitutive model for beams is that of a Saint Venant-Kirchhoff material [1,3,8,10,38], where the strain energy functional per unit undeformed volume is defined in terms of the Green-Lagrange strain tensor E as where C is the right Cauchy-Green deformation tensor, λ and μ the Lamé coefficients and I the identity tensor. Alternatively, SV K can be expressed in terms of C as As shown in [30], it is possible to re-express the term tr (CC) in (25) as Introduction of identity (26) into (25) enables to rewrite the strain energy functional in terms of F and J as The conjugate stresses { F , H , J } (6) for this model are obtained as and the first Piola-Kirchhoff (7) results in The Hessian operator [H W ] (9) adopts the following simple expression Unfortunately this model is not polyconvex as can be observed from the lack of positive (semi)-definiteness of the Hessian operator in Eq. (30). Finally, the tensor H eq W (11) is expressed as

Continuum degenerate polyconvex beam formulation
In this section, a continuum degenerate beam model is presented. Contrary to the work in [10,11], where the strain energy functional is presented in terms of the Green-Lagrange strain tensor, the formulation hereby presented exploits the extended set of kinematic strain measures {F, H, J } along with the computational framework outlined in [12,32,33]. In addition, a co-rotational formulation in terms of an extended alternative set of co-rotational kinematic entities is also presented. This co-rotational approach will prove to be a natural link between the continuum degenerate approach and a more classical beam description where well-known engineering strain measures feature as the main kinematic entities.

Beam kinematics
Let us now consider that the elastic body introduced in the previous section is a beam structural element. In particular, it is assumed that the beam in the reference configuration has a straight axis of length L which is completely characterised by the position of a centre line X 0 (s), parametrised in terms of s ∈ [0, L], and an orthonormal triad { D 1 , D 2 , D 3 } where the vectors D 1 and D 2 lay parallel to the cross sectional area A(s) of the beam (with boundary ∂ A(s)) and D 3 is aligned along the undeformed centre line, see Fig. 2. In the following, summation over Greek indexes α ranges from 1 to 2 whereas summation over Latin indexes i ranges from 1 to 3 and Einstein's convention is assumed for repeated indices unless otherwise stated. The beam reference configuration is Fig. 2 Reference and current configuration of a beam defined in terms of the convective coordinates {η α , α = 1, 2} [8] by For the current configuration, an orthonormal triad {d 1 (s), d 2 (s), d 3 (s)} can be introduced with d 1 (s) and d 2 (s) laying on the cross sectional area a(s) (with boundary ∂a(s)) and with d 3 (s) = d 1 (s)×d 2 (s), see Fig. 2. Hence, the motion of the beam can be defined as [8,39] x(η α , s) = x 0 (s) +x(s);x(s) = η α d α (s) , where the orthonormal triads defined above are related via a rotation tensor R ∈ SO(3) defined as R = d i ⊗ D i [1]. It is known that the rotation tensor R can be defined in terms of the Rodrigues parametrisation of a rotation vector θ [40], with associated skew symmetric second order tensor θ = I × θ, 3 as The deformation gradient tensor F (1) derived from the mapping (33) is obtained as [8] where (•) := d(•)(s) ds . Above Eq. (35) in conjunction with (34) leads to a description of the deformation gradient tensor F in terms of the centre line x 0 and the rotation vector θ, typically used in continuum degenerate descriptions [10,11] of beam structural models. Notice that a computational framework built upon the representation (35) of the deformation gradient tensor is not restricted to beams with a straight axis in the reference undeformed configuration.
It can be interesting to re-write F in a more meaningful manner from the physical standpoint. With that in mind, following the algebraic manipulations in [8], the deformation gradient tensor F (35) can be re-expressed as where the beam strain measures andK are defined as withK a skew-symmetric tensor that can be re-written in terms of a vector K asK = I × K . 4 The beam strain measure is known as the axial-shear strain vector whilst K is known as the torsional-bending strain vector. Thus, an alternative representation of (36) expressed in terms of the beam strain Above Eq. (38) is known as the right extended polar decomposition of the deformation gradient F in terms of the beam strain measures { , K } as presented in [8]. Notice that whilst R is a rotation tensor, the non-symmetric strain tensorŪ is not a pure stretch tensor (as in the classical polar decomposition theorem [16]). Note thatŪ is only symmetric in the case that B is colinear with D 3 , which only happens in the absence of shear strain and torsion.
Formulae (38) for the representation of the deformation gradient tensor F are very useful in the context of geometrically exact beam models when the strain energy functional is defined in terms of the classical beam strain measures { , K } [1, 3,41]. Notice that a computational framework built upon the representation (36) of the deformation gradient tensor is restricted to beams with a straight axis in the reference undeformed configuration. A more general expression for this representation of the deformation gradient tensor can be found in [1].

Linearisation of the beam kinematics
As stated above (33), the mapping of the beam is defined in terms of the centre line x 0 and the triad {d i } (or the rota-tion vector θ ). Virtual or incremental variations of the centre line x 0 will be denoted by δu 0 or u 0 , respectively. In addition, virtual or incremental variations of the rotation vector θ will be denoted by δθ or θ , respectively. All of these fields must satisfy appropriate essential and natural boundary conditions at s = 0 and/or s = L [42]. Thus, the virtual or incremental variations of the mapping (33) are then defined as The first and second directional derivatives of the triad with respect to geometry changes are defined by 5 Analogously, for the derivative of the triad with respect to the beam axis, the first and second directional derivatives are The first and second directional derivatives of the deformation gradient F can be obtained as where the directional derivatives of the triad were previously computed in (40)- (41). Notice how the kinematics of the beam introduces, in contrast to the continuum formulation, additional non-linearities into the problem through the nonvanishing second directional derivative of F (recall that this term vanishes in the more general continuum description). The first and second directional derivatives of the cofactor H are computed as 5 These directional derivatives need to comply with the orthogonality Finally, the first and second directional derivatives of the determinant J are computed as Due to the kinematics introduced by the beam problem, further nonlinearities are observed in (43) and (44) with respect to (3) and (4), respectively. Alternatively, a formulation exclusively in terms of the centre line x 0 and the triad {d i } has been presented in references [1,3]. This formulation avoids relating directional derivatives of the triad {d i } with respect to the rotation vector θ as in Eq. (40).

Transition from polyconvex continuum model to polyconvex beam theory
For the case of a polyconvex constitutive model satisfying Eq. (5), it is possible to express the internal virtual work as follows which leads to an identical expression for the first Piola-Kirchhoff stress tensor P as in (7). The tangent elasticity operator is obtained as where where Note that the first term on the right hand side of (47) is necessarily positive for δu 0 = u 0 and δθ = θ and therefore buckling can only be induced by the "initial stress" term associated with the last two terms on the right hand side of (47).
In comparison with Eq. (8), an extra nonlinearity present in above Eq. (47) can be observed, namely the last term on the right hand side. Finally, proceeding as in Sect. 2.2, it is possible to reexpress Eq. (47) in terms of the alternative Hessian operator

Polyconvex co-rotational beam formulation
In this Section, an alternative polyconvex co-rotational formulation is presented. This new approach will prove to be very useful establishing a link between the continuum degenerate description and a classical beam description in terms of the beam strain measures and K defined above. Taking into consideration the right extended decomposition theorem (38) and the necessary objectivity requirement, an equivalent energy functional to W (F, H, J ) can be defined in terms of the strain tensorŪ = R T F (38), its cofactorW = R T H and its determinant J , 6 namelyŴ (Ū,W , J ). Notice that both representations W (F, H, J ) andŴ (Ū,W , J ) have identical expressions in their own arguments. Analogously to (6), a new set of stress variables { Ū , W , J } conjugate to {Ū,W , J } can be defined as 6 Notice that the determinants of F andŪ are equal, namely J = detF = detŪ.
Considering the right extended decomposition theorem (38), the internal virtual work can be written as Conservation of angular momentum in the continuum implies where S is the second Piola-Kirchhoff stress tensor. Hence, an alternative co-rotational representation for the internal virtual work (51) is as follows which results in which leads to the following expression for the co-rotational stress tensorP Following a similar procedure to that of Eq. (47), an alternative expression for the tangent operator of the energy can be obtained in terms of {Ū,W , J } as follows where and with the Hessian operator [HŴ ] denoting the symmetric positive semi-definite operator containing the second derivatives ofŴ Ū ,W , J as As mentioned above, the identical representation of the energy functionals W (F, H, J ) andŴ (Ū,W , J ) with respect to their respective arguments results in identical expressions of the Hessian operators H W and HŴ in their respective arguments. Therefore, both Hessian operators H W and HŴ will be positive semi-definite, provided that W (F, H, J ) and hence,Ŵ (Ū,W , J ) satisfy polyconvexity. The latter representation of the strain energy functional will be further pursued in the next Section.

Relationship with the classical beam formulation
Alternatively to continuum degenerate based formulations, as that presented in Sect. 3, classical beam formulations employ a kinematical description based on the strain measures { , K } by means of the definition of the deformation gradient tensor presented in Eq. (38).
The objective of this Section is twofold. On the one hand, the strain energy per unit undeformed volume (i.e. W (F, H, J ) =Ŵ (Ū,W , J )) will be re-written in terms of the alternative beam strain vector B defined in (38). Further manipulations will lead to the introduction of the strain energy per unit undeformed length in terms of the beam strain measures and K . On the second hand, relationships will be established between the tangent operators of the different energy representations to relate the concept of polyconvexity at a continuum and beam level.

The classical beam strain measures
Manipulation of equations in (37) enable the strain measures { , K } to be expressed in terms of the centre line and the triads [1] as From equations in (59), the first and second directional derivatives of the beam strain measures { , K } follow as and

The strain vector B
The directional derivative of the strain vector B defined in (38) can be written as where use of the tensor cross product formula (130) has been made above.

Co-rotational strain measures
It is possible to obtain explicit expressions of the strain tensor U, its cofactorW and its determinant J in terms of the strain vector B (linearly related to and K ). Recalling equation (38) which gives an explicit representation of the strain tensor U = B ⊗ D 3 + I in terms of the strain vector B, the cofactor W can be evaluated making use of property (143) as where use of the properties (141) and (144) has been made. Similarly, the determinant J can be computed by means of the property (142) and using the above result (63) where use of properties (139), (141) and (144) has been made. The directional derivatives of {Ū,W , J } are obtained from above equations as and

Polyconvex beam theory
It is customary in classical beam formulations to define a strain energy functional per unit undeformed length W b ( , K ) dependent on the classical beam strain measures and K . However, it is insightful to obtain the relationship between the strain energy functional per unit undeformed length W b ( , K ) and the strain energy functional per unit undeformed volume W (F, H, J ) (5) (or its alternative equivalent representationŴ (Ū,W , J )) defined previously.
Moreover, by use of Eq. (38), where the strain tensorŪ is expressed in terms of the strain vector B, in conjunction with the definitions of the cofactor (143) and the determinant of a tensor (142), it is feasible to obtain an explicit expression of the strain energy functional in terms of the strain vector B (linearly related to and K ), namelyW (B( , K )). The relationship between the alternative strain energy functionals can be summarised as

Convexity in terms of the strain vector B
In the context of beam mechanics, an appropriate restriction for a constitutive model is convexity with respect to the beam strain vector B. Essentially, the strain energy per unit undeformed volume of the beam must be a function of the deformation gradient via a convex functionW as whereW is convex with respect to its three variables, namely the 3 × 1 components of B. Following a similar approach to that in Sect. 2.2, it is possible to define a new work conjugate stress vector B to the beam strain vector B as This set of conjugate measures enables the directional derivative of the strain energyW to be expressed as Recalling that the directional derivative of the strain energy fulfils  (56), the tangent operator for the strain energy can be re-expressed in terms of the directional derivatives of the strain vector B as where the Hessian operator [HW ] denotes the symmetric positive semi-definite operator containing the second derivatives ofW (B) Note that the first term on the last right hand side of Eq. (74) is necessarily positive for δu 0 = u 0 and δθ = θ. Therefore, buckling can only be induced by the "initial stress" term in (74) defined by B · D 2 B[δu 0 , δθ ; u 0 , θ ], where the second directional derivatives of the beam strain vector B were previously evaluated in (62).

Polyconvexity in terms of the beam strain measures and K
Equivalently, an appropriate restriction for a constitutive model is polyconvexity with respect to the beam strain measures { , K } [35]. It is well known that this condition ensures ellipticity of the static problem and, equivalently, hyperbolicity of its dynamic counterpart [35]. Essentially, the strain energy b 7 per unit undeformed length of the beam must be a function of the deformation gradient via a convex multivalued function where W b is convex with respect to its 6 variables, namely, the 3 × 1 components of and K . Moreover, invariance with respect to rotations is automatically satisfied since 7 It is possible to relate the strain energy per unit undeformed length b with the strain energy per unit of undeformed volume as b := A(s) d A, where A(s) is a generic beam cross sectional area. and K are defined in the reference configuration. Following a similar approach to that in Sect. 2.2, it is possible to define work conjugates and K to the beam strain measures and K , respectively, as It is easy to realise that represents the axial-shear force vector rotated back to the reference undeformed configuration, whilst K is the torsional-bending moment rotated back to the reference undeformed configuration, work conjugates of the axial-shear strain vector and the torsionalbending strain vector K , respectively. This set of conjugate measures enables the directional derivative of the strain energy W b to be expressed as Recalling the relation between the strain energies W b and W in (68) and hence, between b and and Eq. (71), it is possible to obtain the following expression for the directional The directional derivative of B can be expressed in terms of the directional derivatives of { , K } through Eq. (62), leading to the following modified version of Eq. (79) Note that the tangent operator of the strain energy W b naturally emerges as and the Hessian operator

Relationship between tangent elasticity operators
Use of Eqs. (62) and (81) in the second term on the right hand side of Eq. (82) enables both geometric terms of the tangent operators for W b ( , K ) andW (B) to be related as Hence, the constitutive term for the tangent operators for both W b ( , K ) andW (B) must coincide, namely Use of Eqs. (83) and (62) into (86) finally yields Equation (87) establishes a relationship between the tangent operators H W b and HW for the internal energies W b ( , K ) andW (B), respectively. Notice that these two energies are defined in terms of the classical beam strain measures. In addition, it is also possible to relate these two Hessian operators with that of the energyŴ (Ū,W , J ) (56) defined in terms of the more general continuum strain measures. This can be obtained by re-expressing the tangent operator for the energyŴ (Ū,W , J ) (56) in terms of the strain vector B. Notice first that consideration of Eq. (65) enables the simplification of the second term in the right hand side of Eq. (56) as Regarding the third term on the right-hand side in Eq. and [H W b ] (87) will be positive semi-definite too. Therefore, as expected, satisfaction of the ellipticity condition [19] in the continuum through a polyconvex constitutive law implies satisfaction of the ellipticity condition [35] in the consistently derived beam problem. The opposite cannot be inferred. Further inclusion of suitable growth conditions into the constitutive model guarantees well posedness of the overall problem [19].

Polyconvex constitutive models in classical beam theory
In general, given an arbitrary objective energy functional expressed as a function of the extended kinematic set {Ū,W , J }, it is possible to obtain an equivalent expression in terms of the strain vector B (and hence in terms of { , K }) after some algebraic manipulations. 8

A Mooney-Rivlin constitutive model
As presented in Appendix 2, an expression for the strain energy functionalW (B) can be obtained for a Mooney-Rivlin model (154) as (where use has been made of the function f as defined in (18)) As shown in Appendix 2, the stress vector B is obtained as and the Hessian operator [HW ] as Hence, convexity is guaranteed provided that the material parameters α, β and λ are non-negative. Therefore, conforming to the conclusions obtained in Sect. 4 (96)

A Saint-Venant constitutive model
The Saint-Venant-Kirchhoff strain energy in Eq. (24) can be expressed in terms of the strain vector B using Eq. (38) Introduction of Eq. (97) into (24) leads to a final expression of the strain energyW SV K in terms of B. Successive derivations of this expression yield the following Hessian operator HW Notice that for the particular case of B = −D 3 which is clearly non positive semi-definite. Since positive semi-definiteness of HW is not guaranteed for all values of the strain vector B, it can be concluded that the constitutive model is not convex in B. In order to overcome this shortcoming, it is customary [8] to neglect the higher order terms in (97) (last term on the right hand side of (97)) resulting in a linearised version of the Saint Venant-Kirchhoff model as The Hessian operator HW (75) for the linearised Saint-Venant constitutive model is finally computed as: Therefore, this linearised version of the Saint-Venant-Kirchhoff constitutive model is convex with respect to B and hence, according to Sect. 4.3, W b ( , K ) will be convex with respect to { , K }. However, on the contrary to the Mooney-Rivlin and Neo-Hookean models, the energy functional lin SV K still fails to satisfy adequate growth conditions, namely the coercivity condition 9 [30,36] lim Hence, under high compression scenarios, the Saint Venant-Kirchoff constitutive model does not perform well. Given the limitations of this model in large strain scenarios, more suitable constitutive laws need to be considered.

Variational formulation
The objective of this section is to succinctly present the variational formulation of the geometrically exact beam theory following a continuum degenerate approach. This will be used as a starting point when describing the discretisation strategy employed. As a starting point, let us consider the following standard total energy minimisation variational principle

P D α d A ds
where the first Piola-Kirchhoff stress tensor P is evaluated via Eq. (7) in terms of the conjugate stresses F , H and J . An iterative 11 Newton-Raphson process is applied to obtain the solution. This is usually achieved by solving a linearised system for the increments {u 0 , θ } as [1] where the update of the position of the center line x 0 and of the rotation matrix R is obtained via and with the rotation matrix R ( θ ) obtained particularising Eq. (34) to the incremental rotation vector θ . Finally, in the absence of follower loads, the second derivative of the total energy functional is given by where the tangent operator is evaluated using Eq. (49).

Finite element discretisation
Traditional finite element discretisations based on rotations do not satisfy objectivity [1,2,4,6]. In order to circumvent this shortcoming, a formulation in which the nodal degrees of freedom correspond to displacements of the centre line x 0 and the director triad {d 1 , d 2 , d 3 } is considered following reference [1,6]. In this approach, x 0 and the director vectors are interpolated as where N a u and N a d are standard Finite Element shape functions. Using a Galerkin approach Particularisation of Eqs. (40) to the triads d k collocated at the nodes of the underlying mesh, it results in where the cross product property a × (b × c) = b (a · c) −  c (a · b) has been used in Eq. (112) in order to obtain the expression for D 2 d a i [δθ a ; θ b ].

Discretised kinematics
Based on the discretisation presented and making use of Eq. (112), the discretised form of the first and second directional derivatives of F in Eq. (42) can be obtained as with

Discretisation of the weak forms
Substitution of the discretisation expressions for x 0 and d k and its virtual variations in (110) and (111) respectively, into (106) leads to the final discretised form of the principle of virtual work as where P is evaluated according to Eq. (7). Notice that the first term in the right hand side of Eq. (115) can be identified as the weak form of conservation of linear momentum whereas the second term can be identified as the weak form of the conservation of angular momentum.

Discretisation of linearisations
Bearing in mind Eq. (49), the discretisation of the second directional derivative of the potential M in Eq. (109) will be presented in this section. For instance, second derivatives with respect to the displacement of the centre line are discretised as where g = H + J F. Equation (116) can be equivalently expressed as where the stiffness matrix K ab uu is obtained as where E is the third order permutation tensor. It is easy to observe that the second (geometric) term in (118) vanishes, leading to the final expression for K ab uu as Similarly, the discretisation of the cross derivatives of the potential M with respect to changes in the displacement of the centre line and rotations would lead to an expression of the form which can be written equivalently as Finally, the discretisation of the second derivatives with respect to changes in rotations is obtained as Equation (123) is equivalent to D 2 M [δθ a ; θ b ] = δθ a · K ab θθ θ b , with the stiffness matrix K ab θθ being defined as (124)

Numerical examples
Note that it is not the primary aim of this paper to propose a new finite element discretisation. Nevertheless, the objective of this section is to present a series of numerical examples in order to show the applicability of the variational and computational frameworks presented in Sects. 3, 5 and 6. For all the examples included, linear shape functions have been used for the interpolation of the centre line and the director triad as presented in [1]. The consideration of dynamic effects for one of the examples presented has been implemented in a straightforward manner [4] and thus, it is not further discussed. Reduced numerical integration [1] of the internal virtual work contributions has been carried out in order to alleviate locking effects, whilst exact integration of the inertial contributions has been implemented. For the last example presented, a comparison of the beam model against a sophisticated mixed continuum polyconvex computational framework as reported in [12] is carried out.

Bending test
This example consists of a straight cantilever beam of length L with squared cross sectional area of side 1 25 m. This example has already been reported in [1,3] using a linearised Saint Venant-Kirchhoff constitutive model (100) with λ lin = 1 N/m 2 and μ lin = 1 2 N/m 2 the Lamé coefficients of the material defined in the reference configuration. With these parameters, the Young's modulus is computed For this example, four nonlinear constitutive models are compared by employing the computational framework described above. For consistency with the reference [1,3], the material characterisation of the different constitutive models employed is carried out at the origin utilising the Lamé coefficients λ lin and μ lin defined above. Specifically, the four nonlinear constitutive models are: (i) a Saint Venant-Kirchhoff model defined in (27) The strain energy defined defined above in Eq. (125) has been extracted from a modified Mooney-Rivlin model given in [29]. Material characterisation renders the following equivalences between the respective material parameters and those of the linearised model (i.e. λ lin and μ lin ), as In order to close the definition of the material paremeters for the Mooney-Rivlin and modified Mooney-Rivlin As reported in reference [43], for a linearised Saint Venant-Kirchhoff model (100) with Poisson ratio ν = 0 and an applied moment of value M = M max , the beam closes on itself forming a closed loop which can be described with an available analytical solution. In general, for a nonlinear constitutive model (or when the Poisson ratio ν is not equal to zero), the exact closure of the beam configuration cannot be a priori guaranteed. Figure 4a shows the deformed shape of the beam for the general nonlinear Saint Venant-Kirchhoff model (27). As can be observed in this figure, there can be observed a small interpenetration due to the nonlinearity of the considered model.
With respect to the stress distribution, Fig. 4b shows the contour plot of the hydrostatic pressure for the Saint Venant-Kirchhoff constitutive model. For the other three constitutive models, the hydrostatic pressure distribution is almost identical, with hardly any differences, and is thus not displayed. This can be explained due to the moderate strains undergone through the deformation process. For this particular example, the choice of a specific model does not seem to be relevant for the overall solution.
The objective of this example is also to demonstrate the p-order of accuracy of the formulation, as a function of the chosen finite element approximation spaces. For this purpose, and particularising for the Mooney-Rivlin model, the beam is initially discretised with 80 elements and, subsequently, h-refinement is carried out generating a total of 3 discretisations. As a closed form solution is not available for this problem (due to the nonlinearity of the constitutive model adopted), the finest mesh is used to generate numerically the so-called "benchmark" solution for comparison purposes. The two coarsest meshes are compared against the finest mesh. The error between the benchmark solution and the other discretisations is measured in the L 2 norm for all the unknown variables. Let us define for a tensor (e.g. scalar, vector or second order) field, the L 2 norm as associated with the magnitude of the tensor field ζ . Although the integrands associated to the internal virtual work have been underintegrated along the centre line, the integral defined in Eq. (127) for the evaluation of the error is computed exactly. Notice that in the evaluation of the error norm, a reconstruction of the continuum is carried out. In our case, ζ can be any of the kinematic or kinetic variables, namely x, 12 x 0 , R = d i ⊗ D i F , H and J . This enables the definition of the following error norm where ζ b stands for the benchmark solution and ζ i the solution of the i-th mesh, with i = 1, . . . , (b − 1). This can then be used to assess the convergence of the algorithm under h-refinement. Figure 5a shows the order of accuracy of the variables x, x 0 , R, F , H and J . The convergence observed is p + 1 in x, x 0 , R and J , and p for F , H . As expected, primary variables lead to p + 1 order of convergence whereas, with the exception of J for this particular example, derived magnitudes (i.e. stresses) lead to a reduced order of convergence p. For completeness, the quadratic convergence of the Newton-Raphson algorithm is shown in Fig. 5b.

Beam with slope discontinuity with zero Poisson ratio
In this example, a series of interconnected beams with a Young modulus of E = 10 6 N/m 2 and a Poisson ratio of ν = 0 are considered. The geometry and boundary conditions of the problem are shown in Fig. 6a. In order to model the connection between the beams at any angle, continuity of the incremental rotation angle θ is strongly imposed. This example has been presented in [3] using the linearised Saint Venant-Kirchhoff constitutive model (100). In this case, the Mooney-Rivlin constitutive model defined in Eqs. (12) and (13) is considered, with material parameters obtained after material characterisation in the undeformed reference configuration and the consideration of α M R = β M R . 12 Note that x is computed according to Eq. (33) namely x = x 0 +x.  (12) and (13). Two meshes analysed: 120 and 240 linear elements. Reference solution taken from a 480 element mesh The objective of this example is to demonstrate the use of a polyconvex model with null Poisson ratio in the reference configuration for the case of a non-straight beam configuration. Notice the excellent agreement between the results presented in reference [3] for a linearised Saint Venant-Kirchhoff model and those for the current Mooney-Rivlin model in Fig. 6b. Figure 7 shows the contour plot for different representative variables, namely the hydrostatic pressure p, the Jacobian J and different stress components (i.e. σ 11   {ox 1 , ox 2 , ox 3 } are coincident). In order to emphasize the high deformation that the beam is subjected to, a shadowed representation of the beam in its undeformed configuration is included. For visualisation purposes, the different spans of the structural system in the deformed configuration are shown slightly disconnected, enabling a clearer observation of the contour plot of the different variables.
Finally, Fig. 8 shows the expected p + 1 order of accuracy for the primary variables x, x 0 and R, and the reduced p order of accuracy for the derived magnitudes F , H and J , where two meshes of 120 and 240 elements are compared with a reference solution defined from a mesh of 480 elements.

Constrained torsion-compression example: coercivity of the models
In this example, a straight beam is subjected to high compressive and torsional effects. The objective of this example is to show the un-realistic behaviour of the Saint Venant-Kirchhoff (24) model specially under high compression scenarios in contrast to that of the Neo-Hookean, Mooney-Rivlin (12)-(13) and the modified Mooney-Rivlin (125) models.
The beam configuration and geometry and the boundary conditions of the problem are shown in Fig. 9. The free end is subjected to a compressive force F = −9 × 10 −1 N in the O X 3 axis and to a torsional moment M = −2.7 × 10 −2 Nm about the O X 3 axis. The displacement of all the nodes along the centre line is constrained to remain aligned with the straight axis of the undeformed beam (i.e. x 0 · D α = 0). In this way, physical buckling is prevented and the lack of coercivity of the models can be explored. The material properties in the reference configuration are E = 1 N/m 2 and ν = 0.35. Figure 10a shows the contour plot of σ 12 (notice that axis in reference {O X 1 , O X 2 , O X 3 } and deformed {ox 1 , ox 2 , ox 3 } are coincident) and the deformed configuration of the beam for 100 % of the total load for the Neo-Hookean, Mooney-Rivlin and modified Mooney-Rivlin models. For the Saint Venant-Kirchhoff model, only 30 % of the total loading has been applied. Figure 10b-d shows the contour plot for some representative conjugate stresses for the differ-  Fig. 10c) whilst the conjugate stress J vanishes for the Saint-Venant model (see Fig. 10d). Figure 11 displays the curve relating the absolute value of the displacement in the O X 3 axis of the free end of the beam versus the absolute value of the force F applied. For small strains, the constitutive response of the four constitutive models is, as expected, almost identical. However, the lack of the required coercivity and convexity requirements for the Saint Venant-Kirchhoff constitutive model leads to the unphysical buckling observed in Fig. 11.

Static twisting column
This example includes the twisting of a cantilever beam of length L = 6 m and a squared cross sectional area of side a = 1 m as shown in Fig. 12. The beam is clamped at X 3 = 0 and subjected to a torsion on its free end, namely X 3 = 6. The torsion at the free end is generated through Dirichlet boundary conditions as follows where θ is the angle of rotation. As can be observed, the section is not restricted to in-plane torsion and zero Neumann boundary conditions are imposed normal to the cross sectional area. tively. A similar example has been presented by the authors in previous reference [12]. The objective of this example is to benchmark the current beam formulation testing the capabilities of the beam representation against a very robust and precise continuum formulation. In particular, a mixed formulation associated to a special Hu-Washizu type of variational principle as presented in reference [12] is considered for that purpose. In this mixed formulation, the unknowns are displacements x, the fibre, area and volume maps {F, H, J } and their stress conjugates { F , H , J }. Regarding the selection of functional spaces: continuous quadratic interpolation of the displacement field x, piecewise linear interpolation of the strain and stress fields F, H, F and H and piecewise constant interpolation of the Jacobian J and its associated stress conjugate J are considered. The Mooney-Rivlin model defined in (12)-(13) has been considered. Despite the obvious differences between beam and continuum formulations (i.e. different kinematical description, different interpolation spaces), Fig. 13 shows reasonable agreement between both beam and continuum representations in terms of tangential stresses. This agreement is excellent for stages of the deformation in which the kinematical assumptions of the beam model are applicable. Hence, when warping of the cross section of the column is pronounced, as in Fig. 13c, the comparison is still reasonable but not as accurate, as expected. The warping of the cross sectional area shown by the continuum representation leads eventually to the buckled configuration represented in Fig. 14

Dynamic twisting column
Finally, a beam with geometry depicted in Fig. 12 is considered. In this time dependent problem, an initial angular velocity 0 = 35 sin π X 3 2L D 3 compatible with the boundary conditions is prescribed. The material properties of the beam are compatible with a shear modulus and Poisson ratio in the reference configuration of E = 0.0179 GPa and ν = 0.3, respectively. A similar example has been presented by the authors in reference [32]. The objective of this section is to benchmark the current beam formulation with the continuum formulation described in Sect. 7.4.1 in a time dependent problem. For both the continuum and beam degenerate problems, a generalised alpha method time integrator is employed [44].
The Mooney-Rivlin model defined in (12)-(13) is considered with material parameters defined in (126). Figures  15 and 16 compare well in terms of displacements between the continuum and beam descriptions. In addition, Figs. 15 and 16 also show a good agreement between the contour plots of the tangential stress σ 31 (notice that axis in reference {O X 1 , O X 2 , O X 3 } and deformed {ox 1 , ox 2 , ox 3 } are coincident) and the conjugate stress H 31 , respectively, for the continuum and beam models. However, as expected, the simplifications introduced in the kinematical description of the beam model (lack of contraction/expansion of the cross section of the beam) lead to different results between both continuum and beam models regarding normal stresses.

Concluding remarks
This paper has provided a novel variational and computational approach to formulate polyconvex large strain geometrically exact beam theory, extending the original ideas introduced by Bonet et al. [12]. In addition, three key novel contributions are incorporated in the present work. First, the deformation gradient, its cofactor and its determinant, namely {F, H, J } are used for the first time as the main strain measures in the context of beam theory. Their respective work conjugates, namely { F , H , J } also feature for the first time in a geometrically exact beam formulation. Moreover, their co-rotational strain {Ū,W , J } and stress { Ū , W , J } counterparts are also presented.
For the first time, the strain energy of a Mooney-Rivlin model has been presented in terms of the classical beam strain measures { , K } by taking advantage of the novel algebra associated to the tensor cross product operation presented in [12] and thoroughly detailed in [33]. Notice that this re-expression procedure can be generalised to more complex constitutive models (i.e. anisotropy, higher nonlinearities).
Finally, the authors have shown that polyconvexity of a continuum constitutive model defined via a continuum strain energy functional W (F, H, J ) implies convexity with respect to the classical beam strain measures of the equivalent beam strain energy functional W b ( , K ), stating explicitly the relationship between alternative tangent operators.
Further extension of our work will include multi-physics electro-magneto-mechanical effects.
Remark 2 It is easy to show using simply algebraic manipulations based on the permutation properties of E or the fact that E i jk E kln = δ il δ jn − δ in δ jl , that the above tensor cross products satisfy the following properties (note that v, v 1 , v 2 , w, w 1 and w 2 denote arbitrary vectors and A, A 1 , A 2 , B, B 1 , B 2 and C are second order tensors): where B = + K ×X . The stress vector B can then be obtained following Eq. (70) as where the partial derivative term contained within the second term on the right hand side of above Eq. (155) can be further The Hessian operator HW (75) is then obtained as HW = 2α I + 2β (I + D 3 ⊗ D 3 )