1 Introduction

The geometrically exact rod theory and their finite element implementations have been extensively studied over the last few decades. The pioneering work of Reissner^[1] was the first to provide beam theories capable of describing arbitrary large displacements and the rotation in the plane static problem. A three-dimensional generalization of this kind of beam theory and its finite element formulations were later reformulated in terms of geometric language in the important works of Simo and co-workers^[2-5]. In their theory,a rod is a mechanical model capable of accommodating general three-dimensional motions,and arbitrarily large deformations including extension,flexure,shear,and torsion. The term geometric exact is used to emphasize the exact satisfaction by the model of all the invariance requirements under superposed rigid body motions,and the lack of restrictions on the degree of allowable deformations.

A distinctive feature of this nonlinear rod theory is that kinematic variables include both displacements and rotations. Formally,a configuration of a rod is totally described by a map $S\mapsto {\varphi}(S)\in \mathbb{R}^{3}$ specifying the position of the current line of the centroid and an orthogonal matrix field $S\mapsto{\Lambda}(S)\in {\mathrm{SO}}(3) $,which physically models the cross section. The presence of rotation group has important consequences in the formulation of numerical models. The rotation variables belong to the nonlinear manifold,and thus their treatment requires a special attention. Within the finite element frame work,a particularly important issue is the interpolation of rotation variables,which is needed to provide informations of the rotation field ${\Lambda}$ within each element. It is evident that a standard interpolation of the nodal rotations does not give an orthogonal matrix. In the early work^{[3, 6-7]},approximated rotation fields were obtained by the direct non-linear interpolation of nodal rotation vectors via the Rodrigues formula. However,as first found by Crisfield and JeleniĆ^[8],these commonly used strategies spoil the objectivity of the strain measures with respect to the rigid rotation. As a remedy,in the same work and later in Ref.^[9],they proposed a new co-rotational interpolation approach by decomposition of the total rotation into a reference and a local rotation. The reference rotation is unique for each element,while the standard interpolation is used to the local rotation. They have shown that the proposed method is invariant and path-independent. Romero and Armero^[10] and Betsch and Steinmann^[11] eliminated the explicit reference to the rotation field and reformulated the geometric exact rod model in terms of two director vectors. Through the direct interpolation of the director fields,the difficulties associated with the interpolation of finite rotations have been avoided. IbrahimbegoviĆ and Taylor^[12] stored the rotation variables and curvature strains at quadrature points so as to avoid the interpolation of rotation fields. Their strategy also leads to a strain objective formulation. Ghosh and Roy^[13] tackled this problem by using the quaternion interpolation method from computer vision.

Given the wide spectrum of interpolation methods available in the literatures,we believe that there is still a necessity to provide a more rational and geometric reasoning that reveals the essential relations between discrete approximations and the symmetries of the underlying continuum model. Moreover,we hope that it can serve as a guideline for future development of finite element algorithms. Through serious investigations on the non-objective interpolations,we abstract a general statement phasing in the mathematical language,i.e.,the objectivity of strain measure is preserved in a particular discrete approximation setting,if and only if the interpolation operator commutes with the special Euclidean group SE(3) = SO(3) ${\times}\mathbb{R}^3$ (the group of rigid body motion),or it is SE(3) -equivariant. As the major contribution of the current work,this theorem can serve as a convenient criterion and an insightful guideline to judge and design invariant formulations.

With this geometrical understanding,we reexamine the strategy of finite element implementations first discovered by Simo et al.^[14],when he tried to propose the energy-momentum conserving algorithms,which are used by IbrahimbegoviĆ and Taylor^[12] in their attempt to address the issue of objective finite element formulations. For the sake of convenience,we refer to such formulation as Simo's formulation. Several researchers^{[8, 13, 15]} hold the opinion that formulations based on the interpolation of the iterative and incremental rotation vector are not objective. Since the approach suggested by IbrahimbegoviĆ and Taylor^[12] falls into this category,they think that it is not objective. Even the authors themselves think that their formulation is objective provided that the applied loads are treated as follower forces and moments. However,the theoretic analyses and numerical experiments presented here show that their opinions are not correct.

To understand the objectivity of Simo's formulation,we should keep in mind that knowledge of the rotation field is required only at the quadrature points,since the exact integration is replaced by the numerical quadrature. Accordingly,one can store the rotation field ${\Lambda}$ and the rotational strain measures K at quadrature points. Moreover,by using the interpolation of the incremental rotation vector to update these stored values,the equivariant condition can be enforced at these points. Therefore,the symmetries are preserved exactly by this approximation. As a result,the objectivity of strain measures is guaranteed. In comparison with the widely accepted method proposed by JeleniĆ and Crisfield^[9],despite the extra storage requirement,the required amount of arithmetic decreases considerably.

The outline of this paper is as follows. In the next section,we give the essential mathematical tools from the Lie group theory. The three-dimensional geometrically exact rod model and the frame-invariance of the strain measures are presented in Section 3. We focus on the understanding of the invariant formulation from the perspective of symmetry in Subsection 4.1. With this geometric interpretation in hand,the invariance of Simo's formulation is demonstrated in Subsection 4.2. Several numerical examples are presented in Section 5 to validate these theoretical results and to assess the performance of the concerned algorithm.

2 Mathematical results of rotation group

It has already been mentioned that the presence of the special orthogonal group of $\mathbb{R}^3$ among the configuration space is an intrinsic characteristic of the geometrically exact rod theory. A sound knowledge of the rotation group is indispensable to understand the geometrically exact rod and its numerical treatment. To construct a frame-invariant formulation,the most critical step is the interpolation of finite rotation. The mathematical tools such as the exponential map and its derivative allow us to reason in terms of geometric quantities,and therefore the ability to deal with this interpolation is enhanced. In this section,we only present the results required for describing and analyzing the finite element algorithms concerned in this article. We refer to Marsden and Ratiu^[16],Iserles et al.^[17] and reference therein for further information.

2.1 Elementary facts on rotation group

The rotation group SO(3) is the compact subgroup of the general linear group GL(3,$\mathbb{R}$),consisting of proper orthogonal matrices,

(1)

Its Lie algebra so(3) is the linear space of skew-symmetric matrices,with Lie bracket $[\cdot,\cdot]$ being the ordinary matrix commutator,and is defined as

(2)

The Lie algebra \{so(3) ,$[\cdot,\cdot]\}$ is identified with $\{\mathbb{R}^3,\times\}$ via the standard Lie algebra isomorphism,which is called the hat map or the skew map $(\hat{\cdot}): {\mathbb{R}}^3\mapsto \mathrm{so}(3) $ defined as

(3)

In Cartesian coordinates,this map has the matrix representation,

(4)

The quaternion parametrization of rotation group is crucial in our numerical implementation. Since a quaternion requires only four parameters,substantial computational and storage savings result from working and storing rotation in this format instead of using the matrix form. In addition,they bypass the coordinate singularity typically associated with the three-dimensional parametrization (the Euler angle,the rotation vector). For this reason,quaternions are used widely in the computer graphics,the computer vision,the computational mechanics,the control theory,etc.

Basically,quaternions are four-vector $({{q}_{0}},q)$ assigned with the non-commutative multiplication rule,

(5)

In particular,the unit quaternions,which,as a set,equal the unit sphere $ {\mathrm{S}}^3 \subset {\mathbb{R}}^4 $,form a group under the quaternion multiplication. The classical Euler-Rodrigues formula gives an explicit expression that associates a rotation matrix $ {Q} $ to a unit quaternion $ q \in {\mathrm{S}}^3 $,

(6)

This parametrization of SO(3) in terms of unit quaternion is of crucial importance in computational mechanics,and plays an essential role in our finite element implementation.

2.2 Exponential map

Numerical algorithms dealing with rotation variables often need to convert infinitesimal rotation vectors to finite rotations. The natural choice is the exponential map expm,defined for a matrix group via the standard infinite series,

(7)

A remarkable fact specific to the rotation group is that the exponential map exp: ${\mathbb{R}}^3 \mapsto$ SO(3) admits a closed-form representation given by the Rodrigues formula,

(8)

where ${\textbf{a}} \in \mathbb{R}^3$,$ a = \|{\textbf{a}}\| $,and $E_3$ is a $3\times3$ identity matrix. The $\exp({\textbf{a}}) $ can be regarded as a finite rotation with the rotation angle $a$ about the axis ${\textbf{n}} :{\textbf{a}}/a$.

The Rodrigues formula (8) can be rewritten in terms of unit quaternion. The quaternion corresponding to a rotation of the angle $a$ about the unit vector ${\textbf{n}}$ is simply

(9)

In practical computations,this quaternion form exponential is always used instead of the formula (8) . Considering the procedure for updating the rotation variable $ q \in \mathrm{S}^3 $,the solution of the linearized problem gives the infinitesimal rotation vector $ {\theta} $,in view of Eq.(9) ,the finite rotation associated with this vector takes the form

where $q$ can be updated by the quaternion multiplication (5) : $q \gets \delta q \circ q$. The proposed procedure is advantageous to the approach using the orthogonal matrix multiplication. Since it provides an enhanced computational efficiency due to the quaternion operation,and reduces the storage requirements by handling only four quaternion parameters.

2.3 Derivative of exponential map

In the geometrically exact rod theory,the left representation of spatial derivative of rotation field is none other than the material rotational strain measure. Since interpolations of rotation field rely on the exponential maps (8) and (9) ,the derivative of exponential map $\mathrm {dexp}$ arises naturally in the approximation of strain measures.

For a curve $t\mapsto {\textbf{a}}(t) $ in $ \mathbb{R}^3$,the tangent vector of the curve $t \mapsto \exp({\textbf{a}}(t))$ at the point $\exp({\textbf{a}}(t))$ is given by the derivative of exp at ${\textbf{a}}(t)$,

(10)

Taking the identification of the tangent bundle $T \mathrm{SO}(3) $ with the product ${\mathrm{SO}}(3) \times \mathbb{R}^3$ via right trivialization,the differential of exp,dexp,is a function

(11)

such that

(12)

Moreover,one can define it as the right trivialised derivative of exp,

(13)

where $T_{g} R_{g^{-1}} $ is the derivative of right translation $ R_{g^{-1}} $ at $ g = \mathrm{exp}({\textbf{a}}(t)) $.

The dexp can also be considered as a map of the ${\mathbb{R}}^3$ into the space of linear operators on ${\mathbb{R}}^3$,

In fact,the correspondence $ {{a}} \mapsto \mathrm {dexp}_{{a}} $ is an analytic function of the skew-symmetric matrix $\widehat{{a}}$,

(14)

Just like the function $\exp({{a}})$,it admits a closed-form expression,

(15)

Moreover,the direct calculation of the formal power series $ \widehat{{a}} (\mathrm {expm}(\widehat{{a}})- E_3) ^{-1} $ yields the inverse matrix of $ \mathrm {dexp}_{{a}}$,i.e.,

(16)

The left trivialized version of $ \mathrm {dexp}$ is also useful in several circumstances. From the definition of adjoint operator ${\mathrm{Ad}}_g$,it can be immediately verified that

where $g =\exp({\textbf{a}})$. From the formal power series of dexp (14) ,and with the help of the fact that $ {\mathrm{Ad}}_{\exp({{a}})} = \mathrm {expm}(\widehat{{a}}) $,one can obtain

Thus,we can arrive at the left trivialized derivative of exp by simply changing the sign: ${\textbf{a}} \mapsto -{\textbf{a}}$. Therefore,in companion with (12) ,we also have the relation

(17)

Remark 1 The matrix-valued function $\mathrm {dexp}_{{a}}$ frequently appears in the disguised form in the literatures concerning the geometrically exact rod theory. The operator ${T}({a})$ defined in Ref. ^[7] is equivalent to $\mathrm{dexp}_{{a}}$. In the work of Simo and Vu-Quoc^[18],the operator ${T}({a})$ refers to $\mathrm{dexp}^{-1}_{{a}}$. Crisfield and JeleniĆ preferred to use ${T}({a})$ as $\mathrm{dexp}_{-{a}}$. It seems ambiguous to choose the same notation ${T}$ to denote the different but very closely related object. By adopting the more clear and versatile concept dexp,numerical algorithms involving the differential of exponential map can be expressed in the more consistent and precise way.

3 Geometrically exact rods

We summarize in this section the geometric exact rod model capable of undergoing finite extension,shear,twist,and bending. The goal is to provide the notation and concepts required for the remainder of this article. Particularly,the frame invariance of strain measures and its implication on the rod configurations are discussed in detail. This result elevates the conceptual level at which we can design our algorithms,and allows us to deal with interpolation from a symmetric-based argument.

Throughout this work,we denote $S$ the arc-length parameter of the undeformed rod,and the prime $( \cdot )' $ is used to mark the derivative with respect to $S$. And let $\{{e}_i\}^{3}_{i=1}$ denote the standard fixed orthonormal basis in $\mathbb{R}^{3}$.

3.1 Kinematics

The configuration of a rod is described by a space curve ${\varphi}(S)$,which corresponds to the centerline,and an orthonormal frame $\{{d}_i\}^{3}_{i=1}$,called the directors,that serves to locate the orientation of a cross section of the rod. The cross section plane is spanned by the pair ${d}_1$ and ${d}_2$,and its normal ${d}_3$ is required to satisfy the condition ${d}_3(S)\cdot {\varphi}'(S) > 0$ which limits the amount of shearing and prevents the degenerate rod configurations.

Out of all possible configuration,a reference one is chosen to coincide with the undeformed state. In this configuration,it is convenient to choose the triad $\{{d}^0_i\}^{3}_{i=1}$ with ${d}^0_1$ and ${d}^0_2$ directed along the principal axes of the cross section and ${d}_3^0$ tangent to the curve ${\varphi}_0$. However,it should be noted that the vector field ${d}_3$ need not be tangent to the deformed centroid,so that such model is capable of representing the shear deformation,while ${d}_1$ and ${d}_2$ are still directed along the principal axes of the cross-section.

The orthonormal frame $\{{d}_i(S)\}^{3}_{i=1}$ along the centerline relative to a fixed basis $\{ {e}_i \}^{3}_{i=1}$ in $\mathbb{R}^3$ is uniquely specified by a map $ S \mapsto {\Lambda}(S) \in $ SO(3) ,such that ${d}_i = {\Lambda} {e}_i$. In contrast with the standard Cosserat rod theory,we regard elements in ${\mathrm{SO}}(3) $ as the basic variables and make no further reference to directors field. Accordingly,a Cosserat rod configuration is completely specified by a pair of curves: $S \mapsto {\varphi}(S) \in \mathbb{R}$ and $S \mapsto {\Lambda}(S) \in \mathrm{SO}(3) $.

We denote by $\partial I$ the parameters of two end points of the rod and assume that the possible configurations of the rod have prescribed values $\tilde{{\Phi}} = (\tilde{{\varphi}},\tilde{{\Lambda}})$ on a subset $\partial_\Phi I\subseteq \partial I$. The foregoing discussion implies that any element in the configuration space

(18)

defines an admissible configuration of the rod. By introducing the linear space $\mathcal{V}$ of test functions

(19)

the admissible variations associated with any ${\Phi} = ({\varphi},{\Lambda}) \in \mathcal{Q}$ span the tangent space $T_{\Phi}\mathcal{Q}$ characterized as

(20)

Geometrically,an element in the tangent space $T_{\Phi}\mathcal{Q}$ can be thought as the tangent vector of a curve in $\mathcal{Q}$

(21)

for any $({u},{\theta}) \in \mathcal{V}$. By taking derivative with respect to $\varepsilon$ at $\varepsilon=0$ and observing ${\Phi}_{\varepsilon}|_{\varepsilon=0} = {\Phi}$,we obtain

(22)

In the physics and mechanics literature,the pair $({u},{\theta})$ is traditionally referred to as a variation of the configuration. The first argument ${u}$ is interpreted as an infinitesimal displacement of the line of centroid,and the second argument ${\theta}$ is regarded as a superposed infinitesimal rotation onto the body frame.

3.2 Objectivity of strain measures

The strain measures associated with a configuration $({\varphi},{\Lambda}) \in \mathcal{Q}$ can be defined as

(23)

The translation strain measure $ {\Gamma} $ collects the axial and shear strains. The $ \widehat{{K}} $ is the skew matrix associated with the rotational strain measure $ {K} $,which includes the bending strains and the torsion strain. It can verified that both strain measures are unaffected by spatial rigid body motions and vanish at the reference configuration. Consider an arbitrary rigid body motion $ ({Q},{r}) \in {\mathrm {SO}}(3) \times \mathbb{R}^3 $ which transforms the rod configuration $ ({\varphi},{\Lambda}) $ to $ ({\varphi}_{\ast},{\Lambda}_{\ast}) := ({Q}{\varphi} + {r},{Q}{\Lambda}) $. The strain measures for $ ({\varphi}_{\ast},{\Lambda}_{\ast}) $ can be written as

Since ${Q}^t {Q} = E_3 $,we have $ {\Gamma}_{\ast} = {\Gamma} $ and $ {K}_{\ast} ={K} $.

Moreover,the strain measures completely determine the configuration up to a rigid motion. To see this,we first give a proposition concerning relations of two curves in SO(3) . With the help of the left trivialized tangent transformation of exponential map (17) ,it can be proved by directly taking derivative. These results are also useful in our finite element implementation.

Proposition 1 Let $ g_{1}(t)$ and $g_{2}(t)$ be two curves in SO(3) ,and $ h(t) = g_{2}(t) g_{1}(t) $,the following relation holds for their derivative:

(24)

Moreover,$if\; g_2(t) :=\mathrm{exp}(\zeta(t)) $,where $ \zeta(t) $ is a curve in $ \mathbb{R}^3(\equiv {\mathrm{so}}(3) )$,then the above relation can be reduced to

(25)

where $ \xi(t) $ and $ \eta(t) $ are the left trivialized tangent vectors of $ g_{1}(t) $ and $ h(t) $,

(26)

Proposition 2 If two configurations $ ({\varphi},{\Lambda}),({\varphi}_{\ast},{\Lambda}_{\ast}) \in \mathcal{Q} $ have the same strain measures $ ({\Gamma},{K}) $,then they are related to each other by a rigid body motion,

where $ ({Q},{r}) \in$ SO(3) $\times \mathbb{R}^3 $.

Proof First,let us assume that the two rotation fields $ {\Lambda}(S),{\Lambda}_{\ast}(S) \in {\mathrm {SO}}(3) $ are related by ${\Lambda}_{\ast} ={\Omega}{\Lambda}$,where $ \Omega$ is a curve in ${\mathrm{SO}}(3) $ needed to be determined here. Substituting them into Eq.(24) ,in view of the definition of $ {K} $ (see (23) ),yields

(27)

The condition $\widehat{{K}}_{\ast} = \widehat{{K}} $ implies that $ {\Omega}' = {0} $. Therefore,there is a constant rotation $ {Q} \in {\mathrm{SO}}(3) $ such that $ {\Lambda}_{\ast}(S) = {Q} {\Lambda}(S)$.

With this result in hand,we consider the relation between ${\varphi}$ and ${\varphi}_{\ast} $. From the definition of $ {\Gamma} $ (see (24) ),the condition ${\Gamma}_{\ast} = {\Gamma} $ gives $ {\varphi}'_{\ast}(S) = {Q}{\varphi}'(S) $,which implies that $ {\varphi}_{\ast}(S) = {Q}{\varphi}(S) + {r} $ for some $ {r} \in {\mathbb{R}}^3 $. Therefore,the result is proved.

In the finite element analysis,the objectivity of strain measures does not extend automatically to its numerical approximation. The proposition presented above requires an invariant finite element formulation,which should satisfy the condition that,for an arbitrary configuration $ {\Phi} \in \mathcal{Q} $,the approximation $ {\Phi}^h_{\ast} $ of the transformed configuration $ {\Phi}_{\ast} $ can be obtained by applying the rigid body motion to the approximation $ {\Phi}^h $ of the original configuration. The mathematical idea behind this situation is the equivariant property of approximation operator $ \mathcal{I} $,and will be elaborated next section.

3.3 Governing equations and symmetry

We let ${n}(S)$ and ${m}(S)$ denote the contact resultant force and contact resultant couple,respectively,acting on the cross section at $ S $. On the Neumann boundary $\mathrm{\partial_R}I$,$ {n} $ and ${m}$ are prescribed with the values $ \tilde{{n}} $ and $\tilde{{m}}$. Let $\overline{{n}}(S)$ and $\overline{{m}}(S) $ denote the prescribed body force and couples per unit of the reference length. The well-known local forms of equilibrium equations are given by

(28)

supplemented with the boundary conditions

(29)

The standard conditions $\partial_\Phi I\cup \partial_RI = \partial I$ and $\partial_\Phi I \cap \partial_RI = \varnothing$ are assumed to hold.

For the linear elasticity material,the constitutive law is given by the linear relation between the material contact force ${N} := {\Lambda}^t {n} $ and the material contact couple ${M} := {\Lambda}^t {m} $,and the strain measures $ {\Gamma} $ and ${K}$ are

(30)

where

(31)

The constants $E$ and $G$ are the Young's modulus and the shear modulus,respectively,$A$ is the cross-sectional areas of the rod,and $A_1$ and $A_2$ are the effective areas when taking account of the cross-section distortion. $EA$ is regarded as the elastic axial stiffness of the rod cross-section,and $GA_1$ and $GA_2$ denote the shear stiffnesses along the axes ${{d}}_1$ and ${{d}}_2$. $EI_1$ and $EI_2$ can be interpreted as the elastic bending stiffnesses of the rod cross section relative to principal axes ${{d}}_1$ and ${{d}}_2$. The term $GJ$ is referred to the elastic torsional stiffness.

The finite element formulation is based on the Galerkin weak form of the governing equation (28) . Formally,the construction of the weak formulation of equilibrium equations is obtained by taking the dot product of (28) with an arbitrary test function $ ({{w}},{{\nu}}) $ in $ \mathcal{V} $,integrating over the domain $ [0,L]$ and using integration by parts. The result is

(32)

where $G_{\mathrm{int}} $ is a part of the weak form representing the elastic deform effect,and expressed as

(33)

while $G_{\mathrm{ext}}$ is the virtual work of the applied load and boundary stress resultants,and given by

(34)

The variational equation (32) together with the local constitutive law (30) comprises the weak formulation of the boundary value problem for the geometric exact rod.

We turn to the symmetries possessed by the equilibrium state. From the invariance of strain measures and the constitutive relation,we notice that,under the transition,

the stress resultants ${n},{m} $ will transform to $ {Q}{n},{Q}{m} $,respectively. An obvious,but nonetheless important consequence of this transformation rule is that if $ ({\varphi},{\Lambda}) $ is a solution of (28) with the boundary condition (29) ,then $ ({\varphi}_{\ast},{\Lambda}_{\ast}) $ is also an equilibrium state,with the transformed Dirichlet boundary condition $ ({{Q}}\tilde{{\varphi}}+ {r},{{Q}}\tilde{{\Lambda}})$,and the transformed applied external forces and couples $({Q}{{\overline n}},{Q}{{\overline m}})$,$({Q}{{\tilde n}},{Q}{{\tilde m}})$.

This observation gives us a direct way to confirm the invariance of numerical approximations. In order to verify that a given formulation is invariant under the superposed rigid body motion,we need to check that the solution $({\varphi}^h_{\ast},{\Lambda}^h_{\ast}) \in \mathcal{Q}$ accompanied with the transformed displacement boundary condition and the transformed load condition,can be obtained by superposing the rigid body motion $ ({Q},{r}) $ on the original approximation $ ({\varphi}^h,{\Lambda}^h) $,that is

We will make use of this fact in Subsection 5.3 to demonstrate the objectivity of the formulation concerned in this article.

4 Geometrical understanding of invariant finite element formulations

Let us now turn to understand the invariant finite element formulation from the geometrical point of view. The seminal work of Crisfield and JeleniĆ^[8] discussed it from the invariance of approximated strain measures ${\Gamma}$ and ${K}$ under the superposed rigid motion. Instead,we examine it by studying the symmetry of this nonlinear rod and its implication on the finite element interpolation. The central result can be phrased concisely in the language of group theory,i.e.,a particular finite element formulation preserves the objectivity of strain measures if and only if the interpolation operator commutes with the special Euclidean group SE(3) =SO(3) $\times \mathbb{R}^3$ (the group of rigid body motion) or is SE(3) -equivariant. This geometric reasoning allows us to understand invariant formulations in a conceptual uniform and clear way.

Armed with this result,we examine the invariance of Simo's formulation. In contrary to some researcher's opinions that this formulation is non-objective,we can prove that it is invariant under an arbitrary rigid motion,therefore,the objectivity of the continuum model is preserved. The crux of this proof lies in the fact that since the numerical quadrature is used to compute the involved integrals,the equivariant condition can be relaxed to hold only at the Gauss point. We will cover them in detail later.

4.1 Geometric criterion

A conforming finite element approximation to the weak form involves the interpolation of the arbitrary test functions in $\mathcal{V}$,leading to a finite dimensional approximating subspace $\mathcal{V}^h \subset \mathcal{V}$,

(35)

Here,$N_{\mathrm{pt}}$ is the number of nodal points,and $N_A$ denotes the global shape function corresponding to the node $A$.

At each configuration $({\varphi}_n,{\Lambda}_n)$,the space of admissible variations $T_{{\varphi}_n}\mathcal{Q} \equiv \mathcal{V} $ (tangent space at $({\varphi}_n,{\Lambda}_n)$,see Subsection 3.3) is approximated by $ \mathcal{V}^h $. Accordingly,the incremental displacements and rotation fields $ ({u},{\theta}) $ are interpolated in the standard manner,

(36)

In addition to the standard interpolation for the test function and the incremental field,the deformed geometry,determined by the current configuration $ ({\varphi},{\Lambda}) \in \mathcal{Q} $,is also needed to be interpolated. Since the change of shape and the stress state cannot change when the rod system rotates or translates as a whole,it is natural to require the transformation of approximated configurations to conform with this symmetry.

To be more precise,we abstract a geometric criterion from this physical thought. Let $ \mathcal{Q}^{h} \subset \mathcal{Q}$ be a set of approximating function in the configuration manifold $ \mathcal{Q} $,we may regard an interpolation method as a map $ \mathcal{I}: \mathcal{Q} \mapsto \mathcal{Q}^h $ such that

(37)

where $ \mathcal{I}_1 $ and $ \mathcal{I}_2 $ are interpolations for the position and rotation field. To ensure the frame invariance of an approximation requires that the interpol ant $ \mathcal{I} $ is equivariant respect to the rigid body motion (see Marsden and Ratiu^[16]),that is to say,for all $({\varphi},{\Lambda})$ in $ \mathcal{Q} $,$ \mathcal{I} $ commutates with the action of arbitrary $ ({Q},{r}) \in \mathrm{SO(3) } \times \mathbb {R}^3$,

(38)

(39)

These relations can be clearly depicted in the following commutative diagram:

(40)

The mnemonic diagram shows that approximation of configuration $ ({\varphi}^{h}_{\ast},{\Lambda}^{h}_{\ast}) $ can be obtained through the rigid motion of the approximating configuration $ ({\varphi}^{h},{\Lambda}^{h}) $. Therefore,the invariance of approximated strain measures is guaranteed. And from Proposition 2 in Subsection 3.2,we can conclude that to ensure the invariance of strain measures in a discrete setting,it is necessary and sufficient that the employed interpolation satisfies the equivariance relation.

To appreciate the meaning of conditions (38) ,(39) and understand their constrains on the construction of the numerical approximation,we first consider the standard interpolation of the position field. For an arbitrary $ {\varphi} $,its interpolation takes the form

(41)

where $ {\varphi}^{A} $ is the nodal value of $ {\varphi} $ at the node $A $. Taking account of the interpolation for the incremental displacement (36) ,the approximated position at the step $n+1$ can also be written as

(42)

where ${\varphi}_{n+1}^{A} = {\varphi}_{n}^{A} + {u}^{A}$. Therefore,the map ${\varphi}$ determining the line of centroid is both interpolated and updated in a standard additive fashion. Moreover,from the linearity of $ \mathcal{I}_1 $ and the fact that $ \sum\limits_{A=1}^{N_{\mathrm{pt}}} N_A(S)=1 $,we can directly verify that

(43)

Therefore,the operator $ \mathcal{I}_1 $ satisfies the relation (38) .

When we come to deal with the spatial interpolation of rotation fields,two fundamental difficulties arise here since the SO(3) group is a non-linear manifold. The first difficulty is how to preserve the orthogonality within each element domain. The second is how to fulfill the equivariance condition (39) . To see this,we observe that the standard interpolated matrix field ${\Lambda}^h(S) = \sum\limits_{A=1}^{N_{\mathrm{pt}}} N_A(S) {\Lambda}^{A} $ cannot be still in SO(3) except at nodal points. The most strait forward method to overcome this difficulty is to observe that in the canonical chart of SO(3) (rotation vector parametrization),vectors in the Lie algebra can be safely interpolated in the standard fashion,after taking exponential of the obtained result,an approximation of rotation field can be computed as

(44)

where $ {\chi}^A $ is the rotation vector parameter for the nodal rotation $ {\Lambda}^{A} $: $ {\chi}^A = \exp^{-1}({\Lambda}^{A}) $.

However,this interpolation method does not satisfy the equivariance condition (39) . A concrete example will help us to understand this problem. Given a typical linear element with local shape functions $ N_1(\xi) = \frac{1}{2}(1-\xi) $ and $ N_2 = \frac{1}{2} (1+\xi) $,the nodal rotation is prescribed with the values $ {\chi}^1 = (0,0,0) $ and $ {\chi}^2 = (0,\pi/2,0) $. The interpolation method (44) leads to a approximated rotation field $ {\Lambda}^h(S) $. Given a constant rotation $ {Q}_{\theta} $ with the angle $ \theta $ about the axis ${e}_3 $,let ${\chi}^{a}_{\theta} = {\mathrm{exp}}^{-1}({Q}_{\theta}{\mathrm{exp}}({\chi}^a))$,and $ {\Lambda}^{h}_{\theta}(\xi) $ denote the interpolation field $ \mathrm{{exp}} ( \sum\limits_{a=1}^{N_{\mathrm{en}}} N_a(\xi) {\chi}^{a}_{\theta} ) $. In Fig. 1(a),we give the stereographic projection images of $ {\Lambda}^{h}_{\theta}(\xi) {e}_3 $. The stereographic projection of $ {Q}_{\theta} {\Lambda}^h(\xi) {e}_3$ is a ray with the azimuth $\theta $. When the rotation angle $ \theta$ increases,the deviation of $ {\Lambda}^{h}_{\theta}(\xi) {e}_3$ (solid line) from ${Q}_{\theta} {\Lambda}^h(\xi) {e}_3 $ (dashed line) becomes large. In Fig. 1(b),the rigid rotation angle $ \theta$ increases further,the line $ {\Lambda}^{h}_{\theta}(\xi) {e}_3$ distorts significantly. When $\theta$ reaches $2\pi$,the line ${Q}_{2\pi}{\Lambda}^h(\xi){e}_3$ returns to the original state,however,$ {\Lambda}^{h}_{2\pi}(\xi){e}_3 $ has a bizarre shape. From this numeric experiment,we can conclude that $ {\Lambda}^{h}_{\ast}(\xi)\neq{Q}{\Lambda}^h(\xi)$ for all $\xi \in(-1,1) $. It is this failure of the equivariance condition (39) that spoil the objectivity of approximated strain measure.

4.2 Invariance of Simo's formulation

In order to preserve the objectivity of the continuum theory,Crisfield and JeleniĆ proposed a interpolation strategy based on a standard interpolation (configuration-independent) of the current local rotation with respect to a reference element-wise constant rotation. Its invariance under rigid motion is well understood and has been completely discussed in several papers^{[8, 9, 15]}. For the purpose here,we just focus on Simo's formulation,and examine its invariance from the geometric perspective discussed previously. So far,to our best knowledge,the invariance of this formulation has not been well recognized.

We now discuss how the condition (39) can be guaranteed by the aforementioned method. To do this,it is important to observe that when trying to construct numerical approximations to the variational equation (32) ,knowledge of the rotation field is required only at the quadrature points $\xi^e_g$ inside each element $e$. Therefore,the interpolation could only provide values and derivatives of rotation fields at these quadrature points. Formally,we regard this approximation operator $\mathcal{I}_{2}$ as a map

(45)

Furthermore,if the proper update of these data is performed,it is possible to enforce the equivariance conditions (38) and (39) to hold only at these Gauss points,

(46)

(47)

for all ${Q} \in \mathrm{SO(3) } $. We call the above relation the weak equivariance condition. In view of the definition of $ {\Gamma}$,one can immediately find that this condition ensures the invariance of strain measures at Gauss points. It will become clear that,since rotation fields are updated by incremental rotation vectors ${\theta}$,what matters here is the interpolation of ${\theta}$ fields,there is no necessity for the actual interpolation of the rotation field in this algorithm.

We are now set to detail the update procedure of rotation. Assume that at the step $n$,we have the following values stored at each Gauss point $\xi^e_g$:

Conceptually,for a given interpolated incremental rotation field ${\theta}^h = \sum\limits_{a=1}^{N{\mathrm{en}}}{N_a}{\theta}^a$,a rotation field $ {\Lambda}^{h}_n $ at $\xi^e_g$ is updated via exponential of $ {\theta}^h $. In the practical implementation,we update its quaternion form $ q^{h}_n(\xi^e_g)$ via the following formula：

(48)

where $\Delta q^{h} = (\cos\frac{\theta^h}{2},\; \frac{{{\theta}}^h}{\theta^h} \sin\frac{\theta^h}{2} )$. Taking account of Eq.(42) ,this type of rotation update implies that the rotational strain measures ${K}^{h}_n(\xi^e_g)$ can be updated in the following manner:

(49)

At the end of this update procedure,we store these new values $(q^{h}_{n+1},\: {K}^{h}_{n+1})$ for the next iteration.

With this information of the rotation field and the interpolated centerline $\varphi^h$,the strain measure ${\Gamma}^{h}_{n+1}$ at each Gauss point can be computed directly from its definition (23) ,

(50)

Moreover,the stress resultants $({n}^h,{m}^h)$ at $\xi^e_g$ are immediately computed by direct use of the constitutive equation,

(51)

The procedure described above can ensure the weak equivariance conditions (46) and (47) . First,from (43) ,the standard interpolated centerline ${\varphi}^h $ automatically satisfies (46) . To establish (47) ,we need to prove that values ${\Lambda}^{h}_{n}$ and ${K}^{h}_{n}$ at the current step transform properly under an arbitrary rotation ${Q} \in$ SO(3) . By applying ${Q}$ to the rotation field $ {\Lambda}^{h}_{n}(\xi^e_g)$,it behaves correctly: $ {\Lambda}^{h}_{n}(\xi^e_g) \mapsto {Q} {\Lambda}^{h}_{n}(\xi^e_g) $. Since the rotation vector associated with ${Q}$ is a constant,therefore,${\theta}' ={0}$. By making use of the updated formula (49) ,it is evident that ${K}^{h}_{n}(\xi^e_g)$ remains unchanged.

What remains to be proven is that values ${\Lambda}^{h}_{n+1}$ and ${K}^{h}_{n+1}$ at the next step transform properly. Consider the action of ${Q}$ on the incremental rotation vector ${\theta}^h \mapsto{Q}{\theta}^h$,since ${\mathrm{exp}}({Q}{\theta}^h)$ can be rewritten in the form ${Q} {\mathrm{exp}}({\theta}^h) {Q}^t $,the updated rotation transforms in the following way:

The invariance of ${K}^{h}_{n+1}(\xi^e_g)$ comes from the fact that ${\mathrm{dexp}}_{({Q}{\theta}^h)} = {Q} {\mathrm{dexp}}_{{\theta}^h} {Q}^t$. We observe that the increment in the update formula (49) is invariant under an arbitrary rotation ${Q}$,

Since ${K}^{h}_{n}(\xi^e_g)$ is invariant,${K}^{h}_{n+1}(\xi^e_g)$ is also unchanged under the rotation ${Q}$. Consequently,for an arbitrary rotation ${Q}$,the ${\Lambda}^{h}(\xi^e_g)$ gets multiplied by ${Q}$,but the ${K}^{h}(\xi^e_g)$ remains unchanged.

In conclusion,by storing the rotation field and its derivative (rotational strain measure) at the Gauss points,the information suffices to complete the finite element implementation,and preserves the frame invariance of the continuum model,while bypassing the need for an explicit interpolation of the rotation field.

Remark 2 The authors in the paper^[12] briefly mentioned the frame-invariance of this discrete approximation of the strain measures at the current step "$n$",and did not consider the part concerning the invariance at the next step "$n + 1$". Based on the weak equivariance conditions (46) and (47) ,the mathematical reasoning presented here provides a more solid and complete analysis on the invariance of Simo's formulation.

5 Numerical examples

The algorithms discussed previously have been implemented for the research purpose in the Python programming language. In this section,we use this code to run several illustrative numerical simulations in order to assess the aforementioned formulation and to further confirm its invariance under the rigid body motion.

5.1 Roll-up of beam

Let ${e}_x$,${e}_y$,${e}_z$ be the fixed unit vectors of the right-handed Cartesian coordinate system in ${\mathbb{R}}^3$. A straight rod of length $L$ and bending stiffness ${EI}$ has the initial configuration $({\varphi}_0,{\Lambda}_0) = (S{e}_x,E_3) $,and it is subjected to a concentrated end moment $M$ in the $z$-direction of the Cartesian coordinate system. Since there is no applied external force,the equation ${n}' + \overline{{n}}= {0}$ gives ${n} ={0}$,which implies that ${\varphi}' = {\Lambda} {e}_x$. From these results,the second equation in (28) can be reduced to ${m}' = {0}$. With the Neumann boundary condition ${m}(L) = M {e}_z$,we infer that ${m}(S) = M {e}_z$. Setting ${\Lambda} = {\mathrm{exp}}(\theta {e}_z)$,the last equation implies $\theta(S) = \theta_1S/L $,where $\theta_1 = ML/{(EI)}$ is the rotation angle at the free end. Substituting this result into ${\varphi}' = {\Lambda} {e}_x$,we find

Integrating this equation yields

(52)

It is clear that the centerline is a circular curve with the radius $r = L/\theta_1 = {EI}/M$. An applied end moment $2 \pi{EI}/L $,will force the rod to deform into the full closed circle. With a moment $4 \pi {EI}/ L $,the rod will wind around itself twice.

The numerical simulation of these analytic results with a finite element model consisting of ten linear elements is discussed herein. In Fig. 2,a sequence of deformed configurations with the free end rotation angle $\theta_1 = 0,\frac{\pi}{3},\cdots ,2 \pi$ is obtained by a series of bending moments $ M=EI \theta_1 / L $. It can be checked that these numerical simulation results are consistent with the analytic solution (52) ,and the accuracy of this algorithm is verified.

5.2 Helical beam

In this example,we take a cantilever beam defined similar to the previous example,and subject it to an increasing sequence of moments $M_n$ and the out of plane force $ F_n$ acting at the free end both in the direction ${e}_z$. This example was first studied by IbrahimbegoviĆ^[19] by using a different finite element formulation. The material properties of the cantilever are as follows: the length $L = 10$,the axial stiffness ${EA} = 10^4$,the shear stiffness ${GA}=10^4$,the bending stiffness $EI = 10^2$,and the torsional stiffness $GJ = 10^2$. We choose the maximum $M^*$ of $M_n$ to be $100 \pi $. From the previous discussion,if there is no out of plane force,the applied end moment $M^{*}$ will make the rod wind around itself five times. The simultaneous action of the moment and the force produces a helical deformed shape. The maximum of $ F_n $ is chosen to be $F_{*} = 25$. A sequence of deformation obtained for different $(F_n,M_n)$ is shown in Fig. 3. Note that $ (F_n,M_n) $ are determined by the load parameter $ \lambda_n \in [0,1] $,such that $ (F_n,M_n) = (\lambda_n F_{*},\lambda_n M_{*}) $.

The most surprising phenomenon in this case is that while the rod keeps around itself,the out-of-plane displacement takes oscillatory values as presented in Fig. 4.

The purpose of this example is to show that the quaternion parametrization works in situations with very large rotation angles ($ > 2 \pi $). If we use the rotation vector parametrization,the non-uniqueness problem,which entails the rank-deficiency of the stiffness matrix,will occur when the rotation angle $ \|{\chi}\| = 2 \pi $. As proposed by M\"{a}kinen^[20],this kind of singularity can be tackled through the coordinate transformation,

where ${\chi}= \|{\chi}\|$. However,the use of quaternion parametrization enables us to update rotation fields via quaternion multiplication,without special treatment to tackle singularity typically associated with parametrization using rotation vectors or Euler angles.

5.3 Right-angle cantilever

In the last example,we examine an L-shaped cantilever clamped at one end,subjected to a concentrated vertical force applied at its free end and support rotations about different axes. This example was first studied in Ref. ^[9]. The structure is composed of two equally long mutually perpendicular legs of the length $ L = 10 $,which are parallel to the $x$- and $y$-axes,respectively. It is modeled with five elements in each leg. The material properties are given by these constants,i.e.,$ {EA} = {GA} = 10^6 $,$ {EI} = {GJ} = 10^3 $. We use this example to check objectivity of our formulation.

The simulation consists of two stages. In the first stage,a concentrated force $ {f} = -5{e}_z $ is applied at the free end of the cantilever,and at the end of this phase,the rod will deform to a new state $ ({\varphi}_b,{\Lambda}_b) \in \mathcal{Q} $. Then,a rotation is imposed at the clamped end (the support). As we discuss in Section 3.3,if the applied tip force $ {f} $ follows the support rotations,from the symmetry of the governing equation (20) ,the final configuration will be $ ({Q}{\varphi}_b,{Q}{\Lambda_b}) $.

The results of numerical simulations are reported here. In the first phase,the nodal force applied at the free end nodal point increasing to its final value in five equal load steps produces the configuration $ ({\varphi}^{h}_b,{\Lambda}^{h}_b) $ as plotted in Fig. 5. To confirm the frame invariance of the resulted approximation,we need to verify that the numerical solution corresponding to the support rotation ${Q}$ should be equal to $ ({Q}{\varphi}^{h}_b,{Q}{\Lambda}^{h}_b) $.

When the bending stage finishes,a rotation with the angle $ \alpha $ about the $z$-axis or $x$-axis is imposed at the clamped nodal point,and the concentrated force $ {f}_{\alpha} $ also rotates with the support

In both case,one finds that for an arbitrary node $ A \in [1,N_{\mathrm{pt}}] $,the computed translational and rotational variables satisfy the following relations within the error of the order $O(10^{-14}) $ (see Fig. 5(b),Fig. 5(c)):

These simulations convince us that the rigid rotation of the whole system will not disturb the already established stress state,and discrete nodal variables rotate uniformly without changing the deformed shape. Therefore,the objectivity of the current formulation is verified.

Finally,we consider the case in which the direction of the applied tip force ${f}$ is kept fixed during the support turning around the $ x $-axis. From Fig. 6,it can be observed that when the rotation angle $ \alpha $ gradually increases from 0 to $ 2 \pi $,the cantilever does not rotate as a rigid body,and its shape changes significantly due to the nontrivial work contribution of the fixed force $ {f} $. It is more convenient to consider this process in a frame that follows the support rotation. In this rotating reference frame,this cantilever exhibits a periodic deformation under the action of the cyclic force $ \exp(-\alpha {e}_x) {f}$. It is demonstrated in Fig. 6(b) which plots the deformed shapes during one load cycle in this frame. The difference between the motion in this process and the rigid body rotation is further illustrated in Fig. 7 by comparing their tip positions.

In the numerical simulation,the solution does not repeat itself precisely when the support makes an entire revolution (although the difference which is of the order ${O}(10^{-3}) $ could be in practice quite acceptable). In fact,the results at the end of any whole turn are different from those obtained previously. Figure 8(a) presents the values of the free-end vertical displacement $W$ at the end of each whole revolution. From this figure,we can discern a shift from turn to turn. Some researchers^[13] think that the reason for this non-physical phenomenon is that the discrete approximation cannot preserve the invariance of strain measure under the superposed rigid-body rotation. However,as we have previously shown,the support rotation and the fixed load together have changed the existing state of stress established at the end of the first loading phase,and the whole structure does not rotate as a rigid body. Therefore,a more reasonable explanation of the difference in deformed shapes after each complete revolution should be attributed to the path-dependence feature of the computed response. As noted by Crisfield and JeleniĆ^[8] and IbrahimbegoviĆ and Taylor^[12],this deficiency,found in almost formulations based on the interpolation of incremental and iterative rotations,can also be traced back to the nonlinearity of the rotation manifold. In practical applications,this kind of error can be reduced significantly if the mesh is sufficiently refined. Let ${p} $ and ${p}^{+}$ denote respectively the free-end position before and after a whole turn. The quantity $e = \|{p}^{+} - {p}\|/L$ can be used as a measure of the difference in the computed results at the end of each turn. In Fig. 8(b),it can be observed clearly that $e$ is monotonically decreasing as the number of nodes in each leg increases.

6 Conclusions

The main concern of this article is to examine finite element formulations of the three-dimensional geometrically exact rod,with a particular attention paid to the geometric understanding of the condition that ensures the frame invariance in the resulting discrete approximations. To achieve the main goal,we first provide essential tools from the Lie group theory to facilitate our understanding of the mathematical structure underlying the problem at hand. Next,after carefully studying the properties of the strain measures and their restrictions placed upon rod configurations,we bring forward a novel viewpoint that the concept of equivariance,which has been elaborated upon in this work,is the key to clarifying several subtle points of interpolation process for configurations. With this recognition,the vague statement "approximation inherits invariance" can be rephrased precisely in terms of the mathematical language as "the interpolation is equivariant".

The Simo's formulation concerned here relies crucially on the observation that the knowledge of rotation fields is only required at Gauss points due to the replacement of exact integral by the numerical quadrature. By storing the rotation field and its derivative (the rotational strain measure) at Gauss points,and updating these data properly,the equivariance condition can be enforced to hold at these points. In consequence,the objectivity of strain measures can be guaranteed. Within the framework of this procedure,iterative rotation updates can be carried out in terms of quaternions. From the computational standpoint,the use of quaternion proves to be the optimal choice that avoids singularity typically associated with parametrization employing rotation vectors or Euler angles,provides an enhanced computational efficiency with the quaternion algebra substituted for the matrix multiplication,and reduces the secondary storage requirements by handling only four parameters.

A number of numerical simulations have been presented to demonstrate the performance of this formulation and to further assess our theoretical findings. The last example deserves a special attention,since it can be used to verify numerically the frame invariance of a discrete approximation under superposed rigid body motions.