Two-field formulation of the inertial forces of a geometrically-exact beam element

An independent velocity field is introduced via Legendre transformation of the kinetic energy of a geometrically-exact beam, leading to a first-order system of twice as many governing equations as a one-field formulation. Nevertheless, the new field does not have to be assembled across elements and can be eliminated at the element level, so that the assembled system has the same size as a one-field formulation. Furthermore, because the new field does not have to satisfy the compatibility equations that the original velocity field is subjected to, its finite-element discretization is simpler and leads to simplified inertial forces.


Introduction
Geometrically-exact beam theory is a popular approach for the modeling of slender structural components having one of their dimensions much larger than the other two [1,2].The resulting kinematic description of such component reduces to the position field of a threedimensional reference curve along the beam, on which a rotation field describes the orientation of the beam cross-sections.Rotations form a Lie group, namely the special orthogonal group SO (3), for which the appropriate mathematical framework has been studied extensively.Because the position field and the rotation field of a beam are inherently coupled, the kinematics is best described under the special Euclidean group SE(3) [3][4][5].The governing equations of geometrically-exact beams are typically obtained by variational principles and take the form of six non-linear partial differential equations with second-order derivatives in time and space.
The finite-element method is typically used to discretize spatially geometrically-exact beams and assemble them into complex flexible multibody systems [6][7][8].Nevertheless, the standard finite-element tools are not capable of handling consistently the coupled kinematics and the non-commutative nature of the rotation group.In particular, several issues such as shear locking and objectivity plagued early finite-element discretizations of geometricallyexact beams.Novel developments have been carried out successfully towards their mitigation, see e.g.[5,[9][10][11][12].
Standard finite-element procedures were readily applied early on to the velocity field and seem to lead to satisfactory finite-element discretizations of the inertial forces.Such techniques are however geometrically inconsistent because the resulting discretized velocity field cannot satisfy the Lie bracket with a discretized strain field, namely non-trivial compatibility equations due to the Lie group structure of rotations.In contrast, the interpolation formulas that produce a compatible, discretized velocity field are fairly involved and configuration-dependent.It appears nonetheless that these geometric inconsistencies decrease under mesh refinement and affect the numerical results in a less noticeable manner than those of the strains [13].
Motivated by recent developments by the second author [14,15], we consider here a two-field formulation of the inertial forces of a geometrically-exact beam.It consists of introducing an additional velocity field by Legendre transformation of the kinetic energy.The resulting governing equations are 12 partial non-linear differential equations with the same second-order in space as the one-field formulation, but first-order derivatives in time.Firstorder formulation in time offers some practical advantages in contrast to typical secondorder formulation in time, such that the use of standard first-order time-integration methods, real-valued eigenvalue computation including velocity effects, and the enforcement of kinematic constraints at velocity level.
From the finite-element point of view, several practical advantages result from the use of the additional velocity field.First, because it is an independent vector field, the additional velocity field can be discretized using standard finite element interpolation functions.This leads to a simpler expression of the inertial forces, because the use of the complicated interpolation functions of the velocity field required for the compatibility equations is partly circumvented.Secondly, the additional velocity field can be discontinuous across elements.Consequently, it can be treated as a field internal to the finite elements and eliminated by condensation (Schur complement).Furthermore, a degree of interpolation that is lower than the one used for the configuration can be used, which opens the door for further computational savings.
The paper is structured as follows.The geometrically-exact beam formulation is reviewed in Sect. 2 and the Legendre transformation of the kinetic energy is used to introduce an independent velocity field.A finite-element discretization is discussed in Sect. 3 and a solution strategy for the two-field formation is derived in Sect. 4. Benchmark examples are discussed in Sect. 5.The papers closes with some conclusions and perspectives in Sect.6.

Beam kinematics
A reference, undeformed configuration of a beam is described by x 0 (α) ∈ R 3 , the coordinates of a reference curve of length L parameterized by coordinate α ∈ [0, L], and R 0 (α) ∈ SO(3), a field of rotation matrices that accounts for the orientation of the beamcross sections.In the current, possibly deformed configuration, the beam is described by positions x(α, t) ∈ R 3 and orientations R(α, t) ∈ SO (3), where t ∈ R is the time coordinate.
The kinematics of a geometrically-exact beam can be conveniently described by a frame field H (α, t) : x(α, t), R(α, t) ∈ SE (3), which couples position and orientation [3][4][5]8].Among other possibilities, elements of SE(3) can be represented as 4 × 4 homogeneoustransformation matrices: The deformation field of the beam is measured by comparison of the spatial gradient of the frame field of the current configuration with that of the reference configuration: where κ(α, t) ∈ R 3 is the axial vector field associated with κ(α, t) = R T d α R ∈ so(3) and To ease the notation, differentiation with respect to space is often denoted with a prime: The velocity field of the beam is denoted V (α, t) ∈ R 6 and defined as the time derivative of the current configuration of the frame field: To ease the notation, differentiation with respect to time is often denoted with an upper dot: Similarly, the variation operator is denoted d δ and the associated field To ease the notation, variations are often indicated by a δ in front of the quantity: The commutativity of differentiation operation on the frame field, e.g.
where the right-hand sides of these relationships do not vanish because the composition operation on SE(3) is not commutative, with The following operator also appears in the subsequent development: •T X = X•, where

Strain energy
The strain energy V of the beam is defined as where K(α) is the 6 × 6 symmetric, positive-definite stiffness matrix.

One-field formulation
The kinetic energy K of the beam is defined as where M(α) is the 6 × 6 symmetric, positive-definite mass matrix.

Two-field formulation
In order to treat the velocity field and the frame field as independent fields, a Legendre transformation of the kinetic energy is now introduced: where W (α, t) ∈ R 6 is an independent velocity field introduced by relaxation of the kinematic relationship in time

Equations of motion
The equations of motion of the beam can be derived from Hamilton's principle, which states that the actual trajectory between two time instants t i and t f is such that the action functional is stationary for vanishing variations at t i and t f : where L = K − V is the Lagrangian of the beam.
With the help of relationship ( 6), the variations of the strain energy defined in Eq. (10) reads

One-field formulation
With the help of relationship (7), the variations of the kinetic energy defined in Eq. ( 11) reads Eventually, the six equilibrium equations that govern the frame field along the beam are the following six second-order partial differential equations

Two-field formulation
With the help of relationship (7), taking the variations of the two-field kinetic energy in Eq. (12) yields The first two terms on the right hand side define the two-field inertial forces (compare to Eq. ( 16)) whereas the last term recovers the compatibility conditions (13) between the two velocity fields.Eventually, the governing equations for the two-field formulation are These twelve differential equations involve second-order derivatives in space, namely E , and first order derivatives in time, namely V and Ẇ .

Finite-element discretization
The beam is discretized in space according to the finite element method presented in [16].
The frame field in an element is approximated by an interpolation of a set of uniformly distributed N > 1 nodal frames H k (t) ∈ SE(3) according to the implicit formula where f k (η) are Lagrange's polynomials of degree N − 1, : SE(3) → R 6 is a minimal parameterization of SE(3) [6], H (η, t) is the interpolated frame of the current configuration and η ∈ [−1, 1] are the finite element coordinates.The same interpolation formula is used in the reference configuration, with 0 For convenience, we use bold fonts to refer to the whole set of a node-related quantity.Accordingly, H denotes the set of all the nodal frames H k 's.For vectorial quantities, the bold notation refers to a 6N-dimensional finite-element array that stacks the nodal quantities, which can then be used in standard matrix-array operations.This includes δU and V, which stack all the nodal δU k 's and V k 's, respectively.With a slight abuse of notation, we also use when it is convenient to emphasize that the dependency on H is strictly on the relative nodal transformation k 's.
Spatial differentiation of interpolation formula (21) in the undeformed and current configuration leads to the discretized gradients where , with T the tangent operator associated with parameterization .The discretized deformation measures are then obtained by evaluation of Eq. ( 2): Their variation leads to d δ E = D( )δU and the discretized internal forces are given by where the 6N -dimensional array of discretized internal forces reads The linearization of the discretized internal forces around some * reads where the tangent stiffness matrix is given by Here, for simplicity, the contribution of the derivatives of the interpolation functions, i.e. the second term on the right-hand side, is neglected.

One-field formulation
The consistent interpolation formula for the velocity field is obtained by differentiation with respect to time of interpolation formula (21) and reads where Q( ) is a 6 × 6N interpolation matrix.Similarly, the discretized field of variations is given by δU = Q δU .With these interpolation formulas, the variation of the kinetic energy in Eq. ( 16) becomes where the 6N -dimensional array of discretized inertial forces reads The evaluation of the discretized inertial forces requires a differentiation of Q with respect to time.Furthermore, each occurrence of Q implies an explicit dependence on the frame field itself.The linearization of the discretized inertial forces reads where Contributions of derivatives of Q with respect to the frames are neglected, i.e. the contribution of the inertial forces to the tangent stiffness matrix is neglected.The resulting 6N discretized governing equations for an element take the form

Two-field formulation
In the two-field formulation, the compatibility equations between W and V are enforced weakly through the variational derivation.Like in the one-field formulation, V must be discretized according to Eq. (29) because it is the time derivative of the frame field interpolated via Eq.( 21).Velocity field W , however, is an independent field that can be discretized freely.
In particular, it can be discretized using standard interpolation formulas.Furthermore, it is allowed to be discontinuous across element boundaries and treated as a field internal to an element (see Sect. 4).Consequently, the following discretization of W is considered where N(η) is a 6 × 6M interpolation matrix and W k (t) are nodal values internal to the element.Note that 1 M N , i.e. the number of nodal velocities W may be chosen lower than that of V .In contrast to Eq. (29), Eq. ( 37) is configuration-independent; its differentiation is thus trivial: Ẇ = N Ẇ, δW = N δW and W = N W. With these interpolation formulas, the variation of the kinetic energy in Eq. ( 18) becomes where with Note that in comparison to Eq. (31), Eq. (39) features one less Q and d t Q is not needed.Equation (40) are the discretized compatibility equations, in which M NN is constant and invertible.The linearization of the discretized inertial forces reads where whereas the linearization of the discretized compatibility equations takes the form Contributions of derivatives of Q with respect to the frames are neglected, i.e. the contribution of the inertial forces to the tangent stiffness matrix is neglected.The resulting (6N + 6M) discretized governing equations for an element take the form

Solution strategy for the two-field formulation
The time integration of the two-field formulation can be achieved with the Generalized-α method described in appendix B. The integration method is applied simultaneously to the frames, i.e.H and the related velocities V, and to the independent velocity field, i.e.W and the related time derivatives Ẇ: For one element, the iteration problem to solve the (6N + 6M) governing equations (47) for the corrective increments in Eq. (62) reads where the following relationships between the increments due the integrator were used: U = c 8 V and W = c 8 Ẇ.Because the independent velocity field is defined at the element level and not assembled across elements, it can be eliminated from the solution process by constructing the Schur complement (see appendix C) of the frame field: where and At each iteration, (i) the 6N -dimensional contributions (50) of each element are assembled and the assembled system is solved for the velocity increments V and (ii) the 6M increments Ẇ are obtained subsequently by solving Eq. ( 53) in each element independently.
Eventually, the size of the assembled system is the same as in the one-field case.The only additional cost of the procedure is in the construction the iteration matrix and of the residuals at the element level.

Numerical examples
Two typical benchmark examples are reproduced here to validate the proposed formulation.We use two-noded beam elements and Euler-Rodrigues parameterization [6].Two versions of the independent velocity field are considered: two nodes (linear interpolation functions) and a single node (constant interpolation function).The results are compared to those obtained with the one-field formulation integrated in time with the generalized-α.

Rotating beam
Wu and Haug [17] described a dynamic problem consisting of a cantilevered beam of length L = 8 m connected to the ground via a revolute joint.The root rotation angle at the joint, denoted φ, is prescribed with the following schedule, where T = 15 s, ω = 4 rad/s and τ = t/T .The motion is planar, characterized by geometric nonlinearities and dominated by the inertial forces (centrifugal effect).For τ < 1, the root rotation undergoes a sharp angular acceleration.For τ 1, the root rotation is driven at a constant angular velocity and the motion is characterized by small amplitude vibrations around the nominal rotational motion.The beam's mechanical properties are as follows: The beam is discretized into 10 two-noded elements.Figure 1 shows the transverse and axial components of the beam's tip displacement resolved in a rotating frame attached at the root of the beam.Clearly, both versions of the proposed two-field formulation match the reference solution.

Lateral buckling of a thin beam
We consider the lateral buckling of a thin beam presented by Bauchau et al. [18], see Fig. 2. A straight beam is clamped at one end and a transverse load is applied through a crank and Fig. 2 Illustration of the benchmark of lateral buckling of a thin beam [18] link mechanism at the other end.As the crank rotates, the beam tip is pushed up.When the buckling load is reached, the beam snaps laterally and undergoes violent oscillations.The crank and link are modeled by flexible beams connected by revolute joints.To simulate an initial imperfection of the system, the tip of the beam is connected to the spherical joint via a rigid-body connection of length d = 0.1 mm.The plane of the crank and link mechanism is offset from the plane of the beam by the same distance d.The main beam has constant rectangular cross-section of width 10 mm and height 100 mm.The crank and link have circular cross-sections of radius 12 mm and 24 mm, respectively.All components are made of aluminum, whose mechanical characteristics are: mass density 2680 kg/m 3 , Young's modulus 73 GPa and Poisson's ratio 0.3.The rotation of the crank is prescribed as φ(t) = π(1 − cos πt/T )/2, for t T and φ(t) = π for t > T , where T = 0.4 s.The problem is simulated for 0.5 s with a constant time step size of 0.5 ms.The spectral radius of the generalized-α scheme is set to 0.
The main beam is discretized into 20 two-noded beam elements.Figure 3 shows the vertical displacement and the first component of the angular velocity, both at the beam's mid-span.The results of the two-field formulation are in good agreement with the reference one-field formulation.

Conclusion and perspectives
The two-field formulation of the inertial forces of a geometrically-exact beam presented in this paper yields results that are in good agreement with a similar one-field formulation of reference.While no significant improvement in computational efficiency was observed in our implementation, several elements indicate that the formulation could open the door to an improved framework thanks to the simplified expression of the inertial forces and the discontinuous nature of the independent velocity field (whose elimination can be parallelized over the elements).The use of standard, highly-optimized first-order integrators could bring additional computational improvements.The formulation could also be beneficial for plate and shell finite-element formulations.
GA, GA, GJ, EI, EI ) with axial stiffness EA = 5.03 MN, shear stiffness GA = 1.94 MN, bending stiffness EI = 566 N•m2, and M = diag(m, m, m, m 11 , m 22 , m 33 ) with mass per unit length m = 0.201 kg/m, and moment of inertia m 22 = m 33 = 22.7 mg • m 2 /m (note that GJ and m 11 are unimportant here because the motion is planar.The problem is simulated for 20 s with a constant time step size of 2 ms.The spectral radius of the generalized-α scheme is set to 0.

Fig. 1
Fig. 1 Axial and transverse displacement of the rotating beam.Reference: , two-node velocity: and single-node velocity: