Dynamics of a flexible body: a two-field formulation

A two-field formulation of the nonlinear dynamics of an elastic body is presented in which positions/orientations and the resulting velocity field are treated as independent. Combining a nonclassical description of elastic velocity that includes the convection velocity due to elastic deformation with floating reference axes minimizing the relative kinetic energy due to elastic deformation provides a fully uncoupled expression of kinetic energy. A transformation inspired by the classical Legendre transformation concept is introduced to develop the motion equations in canonical form. Finite element discretization is achieved using the same shape function sets for elastic displacements and velocities. Specific attention is brought to the discretization of the gyroscopic forces induced by elastic deformation. A model reduction strategy to construct superelement models suitable for flexible multibody dynamics applications is proposed, which fulfills the essential condition of orthogonality between a rigid body and elastic motions. The problem of expressing kinematic connections at superelement boundaries is briefly addressed. Two academic examples have been developed to illustrate some of the concepts presented.


Introduction
The floating frame of reference formulation (FFRF) is a well-established method for the dynamic description of an elastic body undergoing arbitrarily large motion. Its principle consists in trying to simplify the problem by splitting the general motion between the rigid body motion of a reference frame and an elastic deformation with respect to that frame. The instantaneous position and orientation of the reference frame are however undeterminate. They can thus be chosen following two main criteria that are hopefully compatible: -The average relative and rotation from the reference configuration should be kept minimal, so that the assumption of linear elasticity can be adopted in the reference frame.
M. Géradin mgeradin@gmail.com 1 TUM-IAS, Lichtenbergstrasse 2a, 85748 Garching, Germany 2 University of Liège, Liège, Belgium -The description of the kinetic energy and of the resulting inertial forces should be made as simple and accurate as possible. This can be achieved through maximal decoupling between kinetic energies associated with rigid and elastic motions, with a minimum of neglected terms.
Following these two criteria leads naturally to the determination of a reference frame undergoing a global motion that does not follow the motion of any point of the material support but follows its own trajectory, thus referred to as a floating reference frame. The floating frame concept was introduced in 1976 by Fraeijs de Veubeke [10] and was soon exploited in the context of spacecraft dynamics e.g. in [5].
The development of the FFRF to describe the elastic behavior of a flexible component in the context of flexible multibody dynamics started about 10 years later [1,17,21], and its description into Shabana's textbook [20] greatly contributed to its popularity. The method is still receiving a lot of interest nowadays, as demonstrated by the constantly growing number of publications on the subject [4,18,22,24,26,27].
Indeed a number of important questions relevant to is efficient implementation remain open, such as: -Is it feasible to choose the floating frame of reference to achieve total decoupling between kinetic energies associated with rigid and elastic motions without any approximation? -If not, which approximations can be made to achieve this goal? -What is the most efficient approach to discretize the relative elastic motion? -Under which circumstances are the gyroscopic forces of elastic origin negligible? -How appropriate are the modal reduction methods classically used in linear structural dynamics in the FFRF context? -What do we define as elastic velocity? -Should elastic displacements and velocities obey the same discretization? -How to select the system unknowns in the most efficient way? -Does the floating center of the reference frame need to be kept as an unknown of the problem?
The development of the FFRF is classically based on the determination of a reference frame such that orthogonality is maintained between rigid body and elastic motions. It results from the assumption that the kinetic energy associated with elastic deformation is minimum. The resulting axes are generally referred to as Tisserand axes. In practice, however, the concept of Tisserand axes is difficult to implement exactly because of the nonlinearity of the orthogonality conditions to be fulfilled.
The FFRF version described in this paper (and already introduced in [12]) aims at getting around this difficulty by adopting an appropriate definition of elastic velocity. The latter takes into account the convection term due to angular rotation around the instantaneous center of mass. Adopting such a representation of the elastic velocity is the key to full decoupling of the system kinetic energy into rigid body and elastic deformation contributions without any approximation. Such a decoupling is, however, submitted to the essential condition of orthogonality between a rigid body and elastic velocity fields.
Developing this new FFRF version imposes treating the elastic displacement and velocity fields as independent, as was already proposed in [12]. It will be shown in the first part of the paper that the derivation of the formulation can be best achieved by introducing a transformation inspired by the Legendre transformation concept [8], from which a canonical form of Hamilton's principle results.
The remainder of the paper proceeds as follows.
The full discrete expression of elastic inertia forces is developed in Sect. 5 devoted to finite element discretization. The matrix form of the gyroscopic forces resulting from elastic deformation is given. Despite their apparently simple form, getting their explicit expression would normally require developing the associated matrices at the finite element creation level. Strategies exist, however, to get them from manipulations on the linear stiffness matrix [12] when the finite element model results from the discretization of a 3D continuum. It is also observed that rotational equilibrium around the instantaneous center of mass is affected by elastic deformation since the angular momentum of the body becomes deformation-dependent.
The development of an appropriate reduction strategy of the elastic finite element model (Sect. 7) is subject to the essential condition of orthogonality. Consequently, the component mode synthesis (CMS) widely used not only in structural dynamics, but also in flexible multibody dynamics [7,11] is not suitable in this case since it generates a modal reduction basis consisting of static and vibration modes obtained from fixed boundary configurations. The model reduction scheme adopted here results thus from the combination of attachment and vibration modes obtained in free-free configuration as was proposed in [13] and [16].
As illustrated in Sect. 8, different ways are possible to interface the superelements resulting from the proposed FFR formulation with other components of a complex multibody system. Most of the connection strategies are based on the use of Lagrange multipliers, but the possibility also exists to make use of the nonlinear kinematic relationship between boundary elastic degrees of freedom and nodal frames. Selecting the most appropriate interfacing method might depend on the architecture of the software environment in which the superelement model is exploited. It is generally advisable to adopt an interconnection method that encapsulates the constraint equations linking nodal elastic displacements to nodal frames inside the superelement model.
The academic examples presented in Sect. 9 consist of simple beam systems undergoing large rotation. In both examples, the numerical model consists of superelement models obtained through assembly of linear beam finite elements and reduced through the methodology described in Sect. 7. The gyroscopic terms resulting from elastic deformation have not been included. However, post-processing the results of both examples shows that the ratio between kinetic energy of elastic origin and total kinetic energy is quite small, so that the influence of the elastic gyroscopic forces is likely to be insignificant in either case. On the other hand, the second example presented highlights the fact that the generally neglected effect of geometric stiffening may become significant at high rotation speeds.

Kinematics of the flexible body
The current position of a given point X of an elastic body undergoing overall rigid body motion ( Fig. 1) can be expressed in the form where x 0 is the current position of the reference point 0, -R 0 is the current rotation about it, -X is the material (body axes) position of any point on the body, and e = e(X, t) is its elastic displacement from the rigid configuration. The velocity vector may then be expressed in the material frame as 1 being the translation velocity of the reference point 0 expressed in the inertial frame, and being the skew-symmetric matrix of angular velocities about 0 expressed in the material frame. Equation (2) can be split in the form where is the velocity field due to rigid motion and is the velocity field due to elastic deformation. The velocity field due to rigid motion (6) may also be written in the form the rigid body modes being described in the material frame by the 6 × 3 matrix

Kinetic energy
The instantaneous kinetic energy of the elastic body is computed as follows: The first term K 1 represents the rigid body kinetic energy. Assuming that the center of mass of the undeformed body is taken as the origin of the material frame, we get V XρdV = 0 (11) so that it can be rewritten as where m and J 0 are respectively the mass of the body and its tensor of inertia in rigid configuration The amplitudes of the rigid body velocities V 0 = R T 0 v 0 and Ω 0 expressed in the body frame are not uniquely determined but depend on the choice of the instantaneous reference axes for the deformed configuration. Following [10], let us thus assume that the choice of the reference axes for the deformed configuration is made so that the kinetic energy associated with elastic motion K 2 is minimum.
For that purpose, the second term of (10) can conveniently be rewritten in the form in which the material velocity V has to be regarded as invariant to a change of reference. Therefore expressing the condition and thus Conditions (17) mean that the translation and rotation momentums due to elastic motion around the rigid reference configuration are null. Owing to definition (9) of the rigid body modes, they also express orthogonality of the elastic velocity field to rigid body motion Reorganizing the last term in (10) yields As a result, and under the assumption that the material axes coincide with the principal axes of the undeformed body, the system kinetic energy K can be expressed in the fully uncoupled form In the next section, we will refer to this expression of kinetic energy in terms of the velocities (v 0 , Ω 0 , V e ) as the dual form of kinetic energy It is indicated by the superscript. The primal form of kinetic energy results from the substitution of kinematic relationships (3), (4), and (5) in (20). It has for expression

Canonical form of the variational problem
The motion equations will be derived from the application of Hamilton's principle stated in the form where V int represents the internal energy associated with elastic deformation and the righthand side describes the virtual work of external forces. Let us first consider the contribution to (23) of the kinetic energy. It is obtained by substituting into (23) the primal form (22) of kinetic energy The variation to be performed on (ẋ 0 , R T 0Ṙ 0 ,ė) involves thus the knowledge of the angular variation vector δΘ 0 = vect(R T 0 δR 0 ) and its time derivative. The material expression of angular variation δΘ 0 is obtained by considering the associated matrix relationship Time-differentiation of the angular variation matrix (25) yields Likewise performing the variation the angular velocity matrix (4) yields so that we get the matrix relationship and the associated vector expression Developing variational expression (24) can be achieved by substituting into it variation expressions (25), (28), and (29) and integrating by parts in time whenever appropriate. It leads however to a rather complex second-order expression of elastic inertia forces due to the occurrence of the gyroscopic terms elastic origin.
Second-order derivation in time of the displacement/rotation fields can however be avoided by introducing a transformation between the primal and dual expressions of kinetic energy inspired by the classical concept of Legendre transformation [8, p. 114]. It exploits the fact that the following forms of kinetic energy are quantitatively equivalent: The transformation reads thus Substituting (31) into variational expression (24) and making use of (21) provides the kinetic energy contribution to the canonical form of the variational problem It is a mixed variational expression in which the kinematic variables (ẋ 0 , R T 0Ṙ 0 ,ė) and the associated velocities (v 0 , Ω 0 , V e ) play equal roles and can be varied independently. Expressing the variation of (32) provides two sets of first-order equations which have the meaning of momentum equations on the one hand and inertia forces expressions on the other hand. Their explicit expression will be given later on in the discrete case.
The other contributions to variational expression (23) are now developed. It is assumed that the elastic deformation of the system is governed by an internal potential energy V int resulting from the integration over the body of the strain energy density W where the strain tensor ij is derived from the elastic displacements by a strain-displacement relationship of type D(.) being a strain operator. Assuming linearity between strains and elastic displacements and linearity material behavior, the internal potential energy is described by the quadratic expression where D is the linear strain operator and E is a symmetric, positive-definite matrix of elastic coefficients. External actions result from two contributions, namely external body forces in the volume and surface traction forces on a part Sσ of the boundary on the one hand and imposed displacements on a complementary part S B on the boundary on the other hand.
The virtual work of external loads resulting from the application of body forces f in V and surface traction forces t on Sσ reads By making use of we get with the external force and torque acting on the rigid body Equation (40) shows that the torque acting about the instantaneous center of mass is deformation-dependent. It can be split in the form Finally, prescribed positionx B on S B will be expressed through a 3 × 1 vector of kinematic constraint functions imposed on the boundary through the Lagrange multipliers method To summarize, the canonical variational expression governing the dynamic behavior of the flexible body is obtained in the form and where we have for convenience introduced the notation to denote the angular velocities computed from time differentiation of the rotation operator.

Discretization
The finite element discretization of the elastic continuum can conveniently be achieved using the same set of shape functions N(X) for elastic displacements e and velocities V e : The discretization of the volume and surface integrals in functional (44) requires the evaluation of the following matrices: -The linear, symmetric, and positive semi-definite stiffness matrix -The symmetric and positive definite mass matrix and -The N × 6 discrete matrix of rigid body modes U such that constraints (18) obey the discretized form -The skew-symmetric coupling matrix of gyroscopic origin and the 3 × n matrix S(q) -The vector external forces on the elastic body expressed in the body frame -The external torque on the instantaneous center of mass being deformation-dependent is expressed as The discretization of the constraint function vector (43) to describe imposed position on the boundary is made in all generality by defining at the N B boundary nodes a set of nodal frames It gives rise to a set of 6 × N B kinematic constraints which we write in the generic form where L B is an N B × N extracting from the global set the boundary degrees of freedom the N B degrees of freedom involved in the connection and where w 0 stands for the nodal frame at the instantaneous center of mass Finally, the discretized form of the canonical variational expression (44) reads

Equations of motion
In order to get the system motion equations from the variation of canonical form (63), we first note the following: -Owing to (29) the variation of the angular velocity W 0 has for expression -By making use of equality (55) we get -Assuming for the time being that the boundary nodal frames (59) are themselves variables of the problem, the variation of (60) can be formally expressed as where B w 0 , B q B , and B w B are Jacobian matrices resulting from the differentiation of (60). The first term in the right-hand side of (66) can be further split into Performing the variation of (63), integrating by parts in time whenever appropriate, taking account of results (64), (65), and (66), and reordering the terms yields We thus get the canonical motion equations of the elastic body in the form under the constraint of orthogonality between elastic velocity field and rigid body motion The following observations can be made. (69h) is an essential condition to be fulfilled to obtain the solution to the set of equations (69a)-(69g). In all generality, it should be taken into account using a set of Lagrange multipliers in order to remove the indeterminacy over the position and orientation of the floating reference frame. An alternative consists in adopting a modal reduction basis already orthogonal to rigid body modes. In Sect. 7 a choice of modal reduction basis is proposed which is the property of orthogonality with respect to rigid body modes while keeping local elastic displacements as generalized degrees of freedom.

Remark 2
The rotational equilibrium equation (69c) at the center of mass is expressed in terms of a deformation-dependent expression of the angular momentum The total torque C 0 (q) applied to the instantaneous center of mass is also deformationdependent.

Remark 3
The momentum equation (69f) is nothing else than a discretized form of the kinematic compatibility relationship V e =ė + R T 0Ṙ 0 e. A possible simplification to it would be to replace it with its collocated varianṫ

Remark 4
The equilibrium equation (69e) of the discretized continuum is modified by the addition of the gyroscopic term G( W 0 )q.

Remark 5
The effective computation of the N × N gyroscopic matrix G( W 0 ) and the 3 × N matrix S(q)r is discussed in [12]. In the specific case of a 3D continuum model, the same shape functions are generally adopted for all three directions so that the mass matrix can be split in the form with an identical mass kernel M in all three space directions. It is easily verified that the matrices G( W 0 ) and S(q)r may then be easily computed from the knowledge of this kernel.

Remark 6
Most of the time, the motion equations governing the flexible behavior of the moving body are simplified by assuming linearity with respect to relative elastic displacements and velocities. Equations (69c), (69e), and (69g) are then simplified as Remark 7 When using the linearized approximation of inertia forces as proposed in Remark 6, the resulting system of equations is easily recast in the classical second-order form. There is however a benefit in keeping the canonical form of the motion equations, in the sense that treating velocities as independent authorizes velocity discontinuity on the boundary of the superelement. Avoiding the time-differentiation of kinematic constraints (69g) required to express velocity continuity at interfaces simplifies significantly the computation of the iteration matrix, but at the cost of an increase in the number of degrees of freedom.

Model reduction
Model reduction of the system of equations (69a)-(69f) can be achieved by adopting a modal reduction basis Y which applies equally to displacements and velocities q = Yq r= Yr, q andr being the reduced sets of displacements and velocities. The choice of the modal basis Y is not arbitrary but is submitted to the following requirements: 1. Condition (69g) of orthogonality of the elastic velocity field r to rigid body modes is an essential condition which ought to be fulfilled by the modal reduction basis The classically used component mode synthesis (CMS) method [9] generates a modal reduction basis composed of static boundary modes and vibration modes resulting from fixation conditions on the boundary. They are by construction not orthogonal to rigid body modes. Therefore the CMS method does not fulfill the essential condition (74) and has to be discarded in the present context. 2. The model reduction method adopted should be able to provide a statically correct solution on the boundary B of the reduced model to allow for a correct transmission of reaction forces resulting from interaction between bodies of the global system. 3. Ideally, the boundary displacements q B should remain part of the reduced setq to facilitate kinematic connections on the boundary shared by adjacent bodies.
As demonstrated in [12], conditions 1 and 2 are fulfilled by adopting a reduction scheme of type where F B is a set of free attachment modes defined on the boundary B, p B is the associated nodal intensities, Ξ is a set of free-free vibration modes, and η is the associated modal intensities. The drawback of reduction scheme (75) initially proposed in [15] is its hybrid character, the p B being unknowns of force type.
In [12], use is made of the fact that, based on a transformation proposed in [16] and [13], Equation (75) can be recast in terms of the conjugated boundary nodal displacements q B in the form The matrices F BB , F IB , Ξ B , and Ξ I involved in the transformation result from a splitting of the matrices F B and Ξ into boundary and internal contributions. All discrete canonical equations (69a)-(69g) are not affected by modal reduction. The modified equations are: -Equation (69c) governing rotational equilibrium of the center of mass: -Equation (69e) governing dynamic equilibrium of the elastic continuum: -The elastic momentum equation (69f): -Since the kinematic constraints apply only to boundary displacements, Equation (69g) is simplified as The matricesK = Y T KY andM = Y T MY are the usual reduced stiffness and mass matrices. Getting the explicit expression of the reduced gyroscopic matricesḠ( W 0 ) and S(q)r without returning to the generation step of the finite element model is not a simple task. Developing a simplified formulation of the gyroscopic terms from the knowledge of the linear mass matrix M and integrating them in system models in which their effect might be significant is still an open subject of research [24]. In practice, the nonlinear gyroscopic contributions of elastic origin to inertia forces can be considered as negligible in most applications of multibody dynamics. However, they still play an important role in specific domains of application such as rotor dynamics. In the examples developed in Sect. 9 linear elastic deformation has been assumed, and gyroscopic effects of elastic origin have been considered as negligible. The splitting of the kinetic energy in the form (20) that the present formulation allows provides a tool to measure the relative contribution of elastic deformation to the total kinetic energy of the system.
On the one hand, looking a posteriori to energy distribution as done in both examples presented allows verifying the correctness of this approximation to inertia forces. On the other hand, it will also be observed in the second example that the nonlinear elastic forces induced by large centrifugal forces may play a significant role in the system response.

Boundary connections
Discretizing the constraint function (43) on the boundary is more easily performed through collocation, thus in an inconsistent way. For that purpose, let us assume that the set of degrees of freedom of the finite element model consists of triplets of nodal infinitesimal displacements d P and possibly (for models including beams, plates, or shells) infinitesimal rotations α P , so that At each node P of boundary B the following set of six constraints has to be satisfied: The variations of the nodal frames w k can be expressed as In [12] it has been shown that, in the vicinity of Φ P = 0, we get for the variation of the constraints (85) In order to describe the connectivity between bodies, let us consider the general case of two bodies sharing the same boundary B. Interconnection can be expressed by defining two sets of kinematic constraints sharing the same nodal frames w B In practice, it can be achieved in different ways.

Mixed interface method
In the mixed interface method, the interface constraints result from collecting information from both interfacing domains. Equations (86)-(87) are combined to eliminate between them the common nodal frame w B to generate the single interface constraint The principle of the mixed interface method is illustrated on Fig. 2.
With the mixed assembly method, the geometrically nonlinear treatment of the interface is done outside the superelement models through collection of information on their current position, rotation, and current local deformation state. The global kinematic constraint (88) is generally imposed through a single set of Lagrange multipliers λ B . There is however no direct access to the global state w B of the common boundary.

Local interface method
The principle of the local interface method (Fig. 3) consists in introducing inside body k the spatial frame w k B through a set of Lagrange multipliers λ k B . In this manner, the superelement Fig. 2 The mixed interface assembly method: a single set of constraints Φ B is imposed to interface the bodies through a common set of multipliers λ B . Each body has thus to communicate its current nodal frame w 0 and elastic deformation q B at the interface can be treated Boolean-wise, and it requires a minimum of complexity for its development.

Spatial coordinates method
The spatial coordinates method (Fig. 4) consists in inverting kinematic relationships (82) to get an explicit expression of the local elastic degrees of freedom or, symbolically, Making use of Equation (90) allows to adopt the nodal frames w B as unknowns on the boundary. The hidden complexity of the method results from the observation that taking the variation of the elastic unknowns q B provides the expression Fig. 4 The spatial coordinates assembly method: the elastic boundary degrees of freedom q k B are eliminated to provide an explicit expression of the nodal frames w k B . The latter are external connectors which can be treated Boolean-wise where w B = (x 0 , R 0 ) is the nodal frame at the center of mass. As a result of the existence of the second term in Equation (91) which corresponds to motion variation at the instantaneous center of mass, a coupling is introduced between equilibrium equations of the local elastic model (Equation (78)) and equilibrium equations at the instantaneous mass center (Equations (69a) and (77)).
The spatial coordinates assembly method is described in detail in [12]. Its main advantages are: -encapsulation into the superelement model of the inherent geometric nonlinearity of the connection between local and global degrees of freedom, -absence of Lagrange multipliers to express boundary connections that can be formulated Boolean-wise, -minimization of the number of degrees of freedom.
The drawback of the spatial coordinates assembly method is the increased complexity in the development of the superelement model.

Applications
Both examples presented in this section are simple beam systems undergoing large rotation. In either case, the discrete model consists of superelements obtained from the assembly of linear beam elements. The boundary connectors of the superelements are in all cases the positions and rotations at end nodes. The gyroscopic inertia forces have not been included because the corresponding termsḠ( W 0 ) and S(q)r are not directly deductible from the reduced mass matrixM for a beam-like model. The superelement models have been built using either method 2 or 3 as described in Sect. 8. Both methods provide identical results. Method 2 turns out to be easy to implement but generates a larger set of dofs, while Method 3 provides a minimum set of dofs at the cost of more complex implementation. Validation is obtained in both cases through comparison with a nonlinear beam finite element solution.

Rotating beam with a spherical joint submitted to driving torque and tip load
The rotating beam example sketched on Fig. 5 was first proposed in Reference [6]. It consists of a rotating beam with the following physical characteristics: length 141:42, mass density 7.8E-03, cross Sect. 9.0, moment of inertia 6.75, Young modulus 2.1 E+06, and Poisson ratio 0.3. The beam is linked to the foundation through a spherical joint. The loading is such that the system is undergoing 3D motion, allowing thus to validate the 3D character of the formalism and assess its correct implementation. The beam is first put into rotation in the OXY plane through the application of a torque of 200 Nm around OZ, and a tip load of 100 N is applied next in order to induce flapping motion in the OZ direction.
The beam has been modeled with two superelements described each by their reference node and two boundary nodes, and six vibration modes per superelement have been included in the model. The model has been assembled following method 3 (spatial coordinates).
The time integration has been achieved using the α-generalized method for first-order systems [14] in the form proposed in [3]. The simulation has been performed on the time interval [0., 50.s] with a time step of 0.01 s and a spectral radius of 0.9. Figure 6 displays the 3-dimensional trajectory of the nodes at mid-length and at the beam tip. The motion remains plane during the first 15 seconds of the simulation. The driving torque and its period of application are such that the rotation angle around OZ is still rather small (0.3 rad) when the tip load is applied. The impulse intensity is such that the angular velocities about OY and OZ after application of the tip load have the same order of magnitude. The system undergoes large rotation since the hub angular orientation at t = 50 s is described by the rotation vector (0.22, −1.54, 1.27) rad. The local vibrations induced by the driving motion after application of the tip load impulse are already noticeable on the trajectory plot, but will be further investigated from the next figures displaying velocities.  Figure 7 compares Ω z the angular velocities at the hub (blue curve) and at mid-length and tip nodes (red and green curves). A rather low oscillation of angular velocities is observed during the first 10 seconds, which results from the rather low intensity of the driving torque. A vibration of much higher intensity occurs after application of the tip load, indication that coupling occurs between the responses to both components of the excitation. The slight linear reduction of the average angular velocity which occurs after application of the transverse tip load results very likely from the change of orientation of the initial rotation axis. Zooming has been achieved on the angular responses through two time windows of 1 second as shown on the figure. It shows that the responses remain very similar in frequency content, and that the responses at hub (in blue) and tip (in red) have always opposite phases.
It is also worthwhile comparing the time evolution of the angular velocities about OZ and OY axes as done on Fig. 8. Due to the impulsive character of the transverse excita- tion, an oscillation of much higher amplitude is observed around OY (green curve). On the one hand, the initial torque induces a linear increase of the Ω z angular velocity up to 0.0295 rad/s. On the other hand, the tip load induces immediately an average angular velocity of −0.0434 rad/s about OY. One can verify that the fundamental oscillation period of Ω y is 1.74 Hz, which corresponds exactly to the fundamental vibration frequency of the hinged-free beam.
In order to asses the quality of the results provided by the superelement model, the same problem has been solved using the ODIN multibody software currently in development at Uliege (BE) [19]. Figure 9 displays the comparison of the angular velocity curves obtained both at the beam hub and beam tip. There is a remarkable agreement between solutions during the torque application phase. Some discrepancy occurs after application of the tip load: both responses remain perfectly in phase, but a linear increase of the velocity oscillation amplitude is observed for the superelement solution.
The next results of interest are the beam deflections at the beam tip. They can be computed from the relationship and are displayed on Fig. 10. In the same manner as observed for angular velocities, to the impulse character of the excitation in the Oz direction corresponds a much higher vibration amplitude but with 0 average value. In the Oy direction, a small quasi-static deflection occurs due to the driving acceleration. After extinction of the driving torque, the average value of the response in the OY direction is also 0. It is easy to measure the apparent frequency of the observed vibration. We measure again 1.74 Hz, which matches exactly the fundamental vibration frequency of the hinged-free beam. Figure 11 displays the distribution versus time of the different forms of energy present in the system, the sum of which corresponds to the external work produced.
Three different periods can be distinguished. During the torque driven period the external work produced (black curve) increases quadratically and stabilizes around 30 Nm after One observes that the bulk of the kinetic energy corresponds to rigid body motion (blue curve), which means that almost all the kinetic energy is carried out by the centers of mass of the superelements. The part associated with vibration around the center of mass (green curve) remains very small. Strain energy (red curve) remains very small during the accelera- tion phase but can reach as much as 28% of the total energy after application of the impulse load. Its variation is compensated by equal change in the amplitude of rigid kinetic energy.
Finally, Fig. 12 displays the relative energy balance over the simulation period. It suddenly increases after application of the impulse loading but nevertheless remains always quite small, indicating adequate energy conservation of the time integration scheme with the set of parameters adopted.

Displacement-driven rotating beam
In 1988, Wu and Haug [25] proposed as a benchmark the rotation beam model described hereafter. It differs from the previous one mainly by its driving mode since it is displacementdriven. Due to the occurrence of geometric nonlinear effects resulting from high rotation acceleration, this example has been widely used by several authors [2,7,23] to assess different modal synthesis techniques suitable for flexible multibody dynamics. We will refer to this example as Haug's rotating beam.
The system consists of a cantilevered beam of length L = 8 m connected to the ground via a revolute joint. The root rotation angle φ(t) at the joint is prescribed as follows: where T = 15 s, ω = 4 rad/s, and τ = t/T . For τ < 1, the system undergoes a sharp angular acceleration which is responsible for nonlinear geometric behavior of the system. For τ ≥ 1, the root rotation is driven at a constant angular velocity and the motion is characterized by small amplitude vibrations. The beam is again described through repetition of the same superelement model having one reference node and two boundary nodes, plus six free-free vibration modes as internal dofs. The beam performs about eight complete turns over the simulation period since the final rotation computed from (93) is 50 rad. The rotation vector at the nodes of the model being mapped on the [0, 2π] interval, the computed rotation does not follow the imposed one as described by Equation (93). Figure 13 compares the time evolution of the rotation computed at hub level (red curve) to the imposed one (blue curve). Figure 14 compares the angular rotation imposed on the hub (in blue) to the angular rotation undergone by the center of mass (red curve). Very low angular deviation is observed both at the center of mass and tip nodes. The response appears to be nearly vibration free, meaning thus from a structural standpoint that the response is quasi-static as seen from an observer in the rotating frame. When zooming on both curves, one notes a slight phase delay at both CM and tip nodes, which results from the bending occurring during the acceleration phase. Figure 15 compares likewise the angular velocity imposed on the hub (in blue) to the angular rotation undergone by the center of mass (red curve). Again, some delay is observed between the hub motion and the center of mass velocity response, but still without any oscillation, confirming the quasi-steady character of the response.
The results of main interest are the relative elastic displacements at the beam tip computed as described before by Equation (92). Figure 16 shows the elastic deflection in the transverse direction (red curve, left) and the resulting beam shortening in the axial direction (blue curve, right). The maximum deflection remains moderate since it corresponds to 7.5% of the beam length. A beam shortening of 0.28% is observed in the axial direction. It has geometric origin and results from the fact that, the axial stiffness being quite high, the beam is nearly inextensible. It turns out however that the geometric stiffening due the centrifugal force has a major effect on the effective bending stiffness of the beam which cannot be properly captured by a single superelement model. It can be expected, as observed by other authors [2,6,23,25], that increasing the number of superelements allows proper transfer to the elastic model of the centrifugal force undergone by the centers of masses.  Figure 17 shows how the bending deflection of the beam changes when increasing the number of superelements. One observes a significant decrease of the maximum deflection with the number of superelements. Its value stabilizes with four superelements in the model, which is consistent with the results presented by other authors. A comparison of our results with those reported in the literature is made in Table 1.  Worthwhile noticing is the difference in convergence behavior between the nonlinear finite element and superelement analyses. On the one hand, increasing the number of superelements has a stiffening effect on the response since splitting the beam allows transmitting an axial force through the connections. On the other hand, geometric stiffening is inherent to the nonlinear formulation. Therefore the flexibility of the model-and thus the maximum tip deflection-increases with the finite element discretization.
Finally, it is observed on Fig. 17 that an oscillation of small amplitude and low frequency remains after stabilization of the driving velocity to its final value of 4 rad/s. The apparent frequency measured from the red curve (4SE) is 0.446 Hz to be compared to the frequency of 0.464 Hz for the cantilevered beam at rest. Figure 18 displays the evolution in time of the different types of energies present in the system. The rigid body kinetic energy (blue curve) largely dominates the behavior of the system. Its maximum theoretical value can be computed as 274,43 Nm. The strain energy (red curve) follows the same time evolution as the elastic displacements at the tip (Fig. 17), but its maximum value ( 0.6 Nm) corresponds to a small fraction of the rigid body kinetic energy. The elastic kinetic energy (green curve) remains negligible compared to the other forms of energy present in the system, confirming thus the quasi-static character of the response.

Conclusion
A two-field variant of the FFR formulation has been presented, which should facilitate an exact description of the inertia forces resulting from arbitrary motion and deformation of an elastic body. The formulation leads to a canonical form of the discretized motion equations, implying thus the use of a first-order time integrator for dynamic simulation. The examples presented demonstrate the effectiveness of the formulation to develop simple superelement models for integration in flexible multibody simulations. They also show that the decoupling of kinetic energy inherent to the formulation provides insight on the relative magnitude of the different contributions to the global system energy. Further work is still needed to develop and implement the discrete expression of gyroscopic forces of elastic origin from the only knowledge of the linear mass matrix of the body. The formulation needs also to be further developed to take into account the second-order geometric stiffening effects that may occur at large rotation speeds.
The two-field FFRF has been implemented in a prototype code developed in the MATLAB ® environment, in which elastic gyroscopic effects are currently neglected. A new version of the code is currently developed using the JULIA © language, which will be used to implement all further developments of the formulation.