See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/356793275

The Motion Formalism for Flexible Multibody Systems

Article in Journal of Computational and Nonlinear Dynamics · December 2021

DOI: 10.1115/1.4053148

CITATIONS READS
2 242

2 authors:

Olivier A. Bauchau Valentin Sonneville

University of Maryland, College Park University of Maryland, College Park
251 PUBLICATIONS 5,939 CITATIONS 34 PUBLICATIONS 491 CITATIONS

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Olivier A. Bauchau on 20 April 2023.

The user has requested enhancement of the downloaded file.

 The Motion Formalism for Flexible Multibody Systems∗
Olivier A. Bauchau∗ and Valentin Sonneville†
∗
Department of Aerospace Engineering, University of Maryland
College Park, Maryland 200742
†
Technische Universität München; Chair of Applied Mechanics

Abstract
This paper describes a finite element approach to the analysis of flexible multibody systems.
It is based on the motion formalism that (1) uses configuration and motion to describe the
kinematics of flexible multibody systems, (2) recognizes that these are members of the Special
Euclidean group thereby coupling their displacement and rotation components, and (3) resolves
all tensors components in local frames. The goal of this review paper is not to provide an in-
depth derivation of all the elements found in typical multibody codes but rather to demonstrate
how the motion formalism (1) provides a theoretical framework that unifies the formulation of
all structural elements, (2) leads to governing equations of motion that are objective, intrinsic,
and present a reduced order of nonlinearity, (3) improves the efficiency of the solution process,
and (4) prevents the occurrence of singularities.

1 Introduction
Flexible multibody systems consist of a collection of rigid and flexible bodies, such as flexible
joints, beams, plates, and shells, interconnected by kinematic joints and force elements [1, 2]. The
configuration of all these components is described by position and orientation coordinates. Similarly,
a frame is described by the position of its origin and the orientation of three mutually orthogonal
unit vectors. Clearly, the configuration of a rigid body and that of a frame are concepts that can
be used interchangeably and so are the words “frame” and “configuration” in the sequel.
Describing the configuration of a kinematic joint involves two frames and the joint imposes
constraints on the relative transformation between these two frames. For instance, revolute joints
limit the relative motion of their two frames to a rotation of arbitrary amplitude about a material
axis attached to the joint. In the sequel, the terms “frame transformation” and “relative motion”
are used interchangeably.
Similarly, a flexible joint or bushing can be viewed as two frames, called “handles,” connected by
a three-dimensional elastic medium. The relative motion of the handles deforms the elastic medium,
generating reaction forces and moments at the handles. Sometimes, flexible joints are subjected
to small deformation only; in other cases, the elastic medium is made of an elastomer that can
∗
Journal of Computational and Nonlinear Dynamics, 17, pp 030801 1–9, March 2022.

1
accommodate finite deformations. The relative motion of the handles remains small in the former
case but is moderate in the latter.
The configurations of beams and shells are described by one- and two-dimensional fields of
frames, respectively, both in their reference and deformed configurations. The deformation field of
the structural element is then related to the relative motion fields between the two configurations.
Clearly, frames and their relative configurations play a fundamental role in the description of
multibody systems. Rigid bodies, beams, and shells are examples of Cosserat solids also called
“directed continua.” The brothers Cosserat [3] were the first to develop a mechanics theory of such
media. Later, Naghdi [4] developed a theory of plates and shells based on these concepts.
This paper focuses on the finite element modeling of flexible multibody systems. In this ap-
proach, a library of finite elements is developed that includes rigid body, flexible joint, lower/higher
pair kinematic joint, beam, and shell elements [5, 6]. Each structural component maps into one or
more finite elements. This modular approach allows the description of systems of arbitrary topology.
The equations of motion of the complete system are obtained by assembling the contributions of
the various elements of the model. This procedure leads to a non-minimal set of unknown variables
constrained by nonlinear algebraic equations that are enforced via Lagrange’s multiplier technique.
The high-level description of the problem presented above underlines the fundamental role played
by frames and their relative configurations. The equivalence between a rigid body and a frame
stems from the assumption that a rigid body cannot deform and hence, a single frame defines its
configuration unequivocally. The description of flexible joints requires the use of two frames and
their relative motion yields the joint deformation measure. The description of lower-pair kinematic
joints also requires the use of two frames and the nature of the joint imposes restrictions on the
allowable joint relative motion.
Similarly, because the cross-section of a beam is commonly assumed to remain rigid [7], a
frame can be attached at any point along its axis and the resulting field of frames defines the
beam’s configuration unequivocally. As the beam deforms, two neighboring frames undergo slightly
different motions and their relative configuration, or, more mathematically, the spatial derivative
of the field of frames, yields the curvature field that characterizes the deformation of the beam.
Finally, because the normal material line of a shell is commonly assumed to remain rigid [8, 9, 10],
a frame can be attached at any point over its mid surface and the resulting two-dimensional field
of frames defines the shell’s configuration. Note that a material line does not define a frame
unequivocally because it can rotate freely about this material line; this rotation, often called the
“drilling degree of freedom” has to be treated carefully. Here again, spatial derivatives of the frame
field yield the curvature field that characterizes the deformation of the shell.
The description of all these structural components in their reference and deformed configurations
follows the same approach. Similarly, the configurations of the system at times tn and tn+1 are
represented by two sets of frames denoted {F|n } and {F|n+1 }, respectively. Based on the equations
of motion of the system, time integration schemes must determine frames {F|n+1 } knowing frames
{F|n } and the associated velocities. Clearly, the task of the time integrator is to determine the
relative motions that bring frames {F|n } to frames {F|n+1 }.
The discussion above shows that frames describe the system but also embody the fundamental
assumption on which the analysis is based. Relative frame motion can be categorized by their
magnitude that can be finite, moderate, or infinitesimal. Finite relative frame motions are found
in kinematic joints, for instance the relative rotation, displacement, or both found in revolute,
prismatic, or cylindrical joints, respectively. Finite frame motions also describe the change of

2
configuration of the system from its reference to its deformed configuration.
Typically, moderate relative frame motions are found between the two handles flexible joints,
generating the deformation of the elastic medium. Similarly, because an accurate simulation of
system dynamics requires the use of small time step sizes, the relative motions from frames {F|n }
to frames {F|n+1 } are moderate. The term “moderate” implies that while the relative motion may
be small, it is not infinitesimal and cannot be treated via a linearization operation.
Finally, infinitesimal relative frame motions are found in structural elements such as beams or
shells. The relative motion of two neighboring sections of the beam that are infinitesimally close is
infinitesimal, giving rise to finite curvature. This situation is treated through the tools of calculus:
the curvature field is the spatial derivative of the motion field. Because frames are functions of
time, their temporal derivative give rises to the velocity field. The distinction between motions
of finite, moderate, or infinitesimal magnitude is important because the mathematical handling of
finite frame motions is fraught with singularities [6, 11]. Consequently, different tools are used to
handle finite and moderate or infinitesimal relative motions.
The manipulation of frames involves the treatment of displacement components that form a
linear space but also that of rotation components, sometimes referred to as orientation or attitude
variables. Rotations form a nonlinear manifold called SO(3), the Special Orthogonal group, and
hence, must be treated with care. Indeed, all parameterizations of rotation involve singularities [6,
11] and furthermore, careless manipulation of rotation can lead to the loss of objectivity of the
formulation.
The equations of motion of a system consist of (1) the equilibrium equations expressed in terms
of the components of vectors, including strain and curvature and linear and angular velocity and (2)
the kinematic relationships that relate these quantities to the position and orientation coordinates
of the bodies. The analyst is then faced with a choice: which frame should these tensor components
be resolved in? Because the equations of motion are invariant under a change of frame, physics
does not provide an answer to this question. Rather, this choice is dictated by convenience and
numerical efficiency.
The sectional forces in a beam are most naturally expressed in a frame attached to the beam’s
cross-section. Indeed, when resolved in the sectional frame, the force components can be interpreted
as the axial and two transverse shear force components. If the same force vector is resolved in the
inertial frame, its components cannot be interpreted easily because the beam’s section is at an
arbitrary orientation with respect the the inertial frame. The same observation applies to the
components of deformation measures and velocity vectors or to the relative motion vectors in a
kinematic joint.
These local frame coordinates are obtained naturally through the consistent use of frame kine-
matics: frames are elements of SE(3), the Special Euclidean group. This rigorous framework, some-
times called screw theory [12], has become the lingua franca in the theoretical kinematics [13, 14]
and robotics communities [15, 16, 17].
Within the mathematical structure provided by the special Euclidean group, position and ro-
tation coordinates are coupled rather than being treated as independent entities. Although most
commonly used in structural mechanics, the latter approach can lead to inconsistent formulations
that lack objectivity and to numerical problems, such as the locking phenomenon in finite element
implementations. By placing frame transformations at the heart of the formulation, the govern-
ing equations exhibit a reduced level of nonlinearity because they are invariant under rigid-body
motions, which de facto, are “filtered out.”

3
 This paper refers to the “motion formalism” as a formalism that (1) uses configuration and mo-
tion to describe the kinematics of flexible multibody systems, (2) recognizes that these are members
of the Special Euclidean group thereby coupling their displacement and rotation components, and
(3) resolves all tensors components in local frames. The approach has been embraced by a few au-
thors [18, 19, 20] for the modeling of rigid multibody systems. Yet, it remains largely ignored by the
mechanics and structural mechanics communities, although it opens the door to the development
of geometrically consistent and efficient solution procedures.
The goal of this review paper is not to provide a detailed description of the finite element
implementation of all the elements found in typical multibody codes but rather to demonstrate how
the motion formalism provides a theoretical framework that unifies the formulation of the entire
library of elements and improves the efficiency of the solution process.
This paper is structured as follows. The framework of the motion formalism for flexible mutli-
body systems is staged in section 2, where frame kinematics is discussed together with the related
topics of representations and parameterizations. The general form of the equations of motion and
a time integration method are described in section 3. The formulation of typical structural compo-
nents is sketched succinctly in section 4.

2 Motion formalism: A game of frames

The kinematics of flexible multibody systems [5, 6, 17] can be described in terms of body-attached
frames and frame transformations that are defined below.
1) Body-attached frames are attached to rigid bodies and are often called “material frames.”
For instance, a material frame can be attached to the center of mass of a rigid body to define the
position of its center of mass and its orientation. When dealing with flexible bodies, this concept
must be generalized. For instance, a material frame attached at a point of a beam defines the
configuration of the cross-section at that point. Although the section is assumed to remain rigid,
thereby defining a material frame, the actual section deforms under the applied loading. Clearly,
the frame now defines an average configuration of the cross-section, typically defined through an
energy closure condition [21, 22]. Body-attached frames are denoted Fp : (p, Rp ), where p is the
material point at the origin of the frame and Rp is the attached basis of mutually orthonormal
vectors ~rp,1 , ~rp,2 , ~rp,3 .
2) Frame transformations express the relationship between two material frames, i.e., the opera-
tion that brings the first frame to the second. In general, this operation involves both displacement
and rotation components. As discussed earlier, frame transformations are associated with relative
motions that can be of finite, moderate or infinitesimal magnitude. The frame transformation from
material frame Fp to material frame Fq is denoted p F q : (pq,
~ p Rq ) where pq
~ is the vector that joints
p
point p to point q and Rq is the rotation that brings orientation Rp to orientation Rq .
To enable actual computations, coordinates are introduced. An inertial reference frame is set
up that is equipped with the system of Cartesian coordinates E : (e, E), where e is the origin of
the coordinate system and E is the reference basis of mutually orthonormal vectors E ~ 1, E
~ 2, E
~ 3.
Position and orientation coordinates with respect to the inertial reference are assigned to body-
attached frames and denoted E F p : (E xp , E Rp ) where E xp ∈ R3 are the coordinates of point p

4
resolved in E and E Rp ∈ SO(3)1 is a matrix whose i-th column stores the coordinates of base
vector ~rp,i resolved in E. Frame transformations are inherently relative to body-attached frames
and therefore described by local coordinates of displacement and rotation, resolved in one of the
p
related body-attached frames. Frame transformation
  F q is resolved in body-attached
 frame Fp and
is assigned coordinates p F q : p dpq , p Rq = E RTp (E xq − E xp ), E RTp E Rq . To ease the notation, the
prescript of body-attached frames is often dropped when there is no risk of confusion.

Remark. Interpreting a body-attached frame as a frame transformation from the inertial frame
of reference to itself should be avoided. Indeed, such frame transformation does not possess an
intrinsic geometric meaning because its definition depends on the arbitrary choice of the inertial
frame. In contrast, body-attached frames exist independently of the selection of a reference frame
and so do the frame transformations between them. Body-attached frames define the configuration
of bodies with respect to the inertial frame, thereby enabling the application of Newton’s laws.
In contrast, frame transformations describe the relative configuration of components objectively
and the coordinates of frame transformations are inherently local. The behaviors of body-attached
frame and frame transformations differ under a change of reference frame: the coordinates of the
body-attached frames change whereas the local coordinates of the frame transformations do not,
see details in section 2.2.

2.1 Representation of frame operations

Coordinates of body-attached frames and of frame transformations are elements of SE(3) = SO(3)n
R3 , the special Euclidean group, which is defined as SE(3) = {(u, R) | u ∈ R3 , R ∈ SO(3)} with
composition operation (u3 , R3 ) = (u1 , R1 ) ◦ (u2 , R2 ) = (u1 + R1 u2 , R1 R2 ), neutral element (0, I),
and inverse elements (u, R)−1 = (−RT u, RT ). These coordinates are related as
E
F q = E F p ◦ pF q (1)

The right-hand side of this equation is also called the composition operation. As suggested by the
notation, relative transformations results from the composition of two body-attached frames, i.e.,
p
F q = E F −1 E p p r
p ◦ F q , and can be further composed, i.e., F q = F r ◦ F q .
Equation (1) expresses the composition operation in a rather abstract manner through group
operation ◦. To perform actual computations, it is desirable to identify an algebraic structure,
called a representation that implements the group operation. Numerous representations have been
studied in the literature [19, 23, 24, 25] and a few are reviewed below.
Motivated by the action on a point, a popular representation is the 4 × 4 homogeneous trans-
formation matrix  
R u
F := H = T (2)
0 1
1
SO(3) = {R ∈ R3×3 | RT R = I, det R = 1} is the special Orthogonal group. A group (G, ◦) is a set G of
elements q together with a composition operation ◦ whose properties are (i) closure, i.e., q1 ◦ q2 = q3 ∈ G, (ii)
associativity, i.e., (q1 ◦ q2 ) ◦ q3 = q1 ◦ (q2 ◦ q3 ), (iii) the existence of a neutral element e ∈ G such that q ◦ e = e ◦ q = q,
and (iv) the existence of inverse elements q −1 ∈ G such that q ◦ q −1 = q −1 ◦ q = e. For SO(3), the composition
operation is the matrix product, i.e., R1 R2 ∈ SO(3), whose neutral element is identity matrix I, and the inverse
elements are R−1 = RT .

5
together with the rules of matrix algebra: the matrix product is the composition operation. The
neutral and inverse elements are the identity and inverse matrices, respectively.
The motion tensor [6, 25] stems from the spatial transformations of lines and is the adjoint
representation of the homogeneous transformation matrices
 
R ũR
F := C = (3)
0 R

together with the rules of matrix algebra. The two above representations involve redundant co-
ordinates because they use 12 coordinates, namely 3 for the position/displacement and 9 for the
orientation/rotation, but the latter must satisfy the 6 orthonormality conditions RT R = I.
In computer applications, a reduced number of coordinates is appealing as it may decrease the
computational cost of the group operation and reduce storage requirements [26, 27]. A popular and
numerically efficient choice of coordinates is the pose, a seven-dimensional array
 
q
F := s = q 
 (4)
u

where q and q are the scalar part and the vectorial part of a unit quaternion [6, 28], i.e., q 2 +q T q = 1.
Unit quaternions, also called Euler parameters, are related to the rotation invariants, angle φ and
the coordinates of unit rotation vector n: q = cos (φ/2) and q = sin (φ/2)n. The composition
operation, s3 = s1 ◦ s2 , becomes q3 = q1 q2 − q T1 q 2 , q 3 = q1 q 2 + q2 q 1 + q̃1 q 2 , and u3 = u1 + R1 (q1 , q 1 )u2 ,
where R1 (q1 , q 1 ) = I + 2q1 q̃1 + 2q̃1 q̃1 . The neutral and inverse elements are {1, 0T , 0T } and s−1 =
{q, −q T , −(RT (q, q)u)T }, respectively. The pose can also be represented in matrix form such that
the composition operation reduces to the standard matrix product
 
A(q, q) 0 0
q −q T
 
s :=  0 R(q, q) u , where A(q, q) =
 (5)
T T q qI + q̃
0 0 1

Note that, due to the unit quaternion condition, A(q, q) is an orthogonal matrix, i.e. AT A = I.
Dual algebra and its connection to the transformation of lines reveal other representations [24,
29, 30, 31]. It is characterized by the introduction of “bookkeeping” parameter ε which satisfies the
rule εn = 0, n ≥ 2. A dual number is denoted a + ε b, where a and b are referred to as the primal
and dual parts, respectively. Dual orthogonal matrices provide an alternative representation of the
Euclidean group
F := R = R + ε ũR (6)
together with the rules of matrix algebra. This equation provides the same information as the
motion tensor (3); orientation/rotation matrix R and matrix ũR now become the primal and dual
parts of an orthogonal dual matrix.
Dual algebra also leads to the dual unit quaternion representation, in which dual quaternions
are made of 8 coordinates and are of the form

F := Q = (q̂, q̂) (7)

6
where q̂ = q + ε q ◦ is a dual scalar and q̂ = q + ε q ◦ a dual array of coordinates that satisfy the two
unit dual quaternion conditions: q 2 + q T q = 1 and qq ◦ + q T q ◦ = 0. Unit dual quaternions are related
to frame invariants: dual scalar φ̂ = φ+ε nT u and dual vector n̂ = n+ε ũn. Indeed, q̂ = cos(φ̂/2) =
cos(φ/2) − ε (nT u) sin(φ/2)/2 and q̂ = sin(φ̂/2)n̂ = sin(φ/2)n + ε (sin(φ/2)ũn + (nT u) cos(φ/2)n/2).
The composition operation, Q3 = Q1 ◦ Q2 , becomes q̂ 3 = q̂ 1 q̂ 2 − q̂ T1 q̂ 2 and q̂ 3 = q̂ 1 q̂ 2 + q̂ 2 q̂ 1 + q̂ 1 × q̂ 2 .
The neutral and inverse elements as (1, 0̂) and Q−1 = (q̂, −q̂), respectively. Dual quaternions can
be cast into a matrix representation such that the composition operation reduces to the standard
matrix product: Q := A(q, q) + A(q ◦ , q ◦ ) with matrix A defined in eq. (5).

2.2 Superimposed Euclidean transformations

A change of reference frame consists in resolving the coordinates of body-attached frames in a
system of Cartesian coordinates S : (s, S) that differ from system E : (e, E). Under such change,
the coordinates of body-attached frames transform as
S
F p = SF E ◦ EF p (8)

where S F E : (S e, S RE ) is a frame transformation made of the coordinates of e and E in S. A

rigid-body motion is an identical operation but E is regarded as fixed and S F E is interpreted as the
motion of a given body-attached frame from one configuration to another.
Both operations are so-called “superimposed Euclidean transformations.” They transform
the coordinates of body-attached frames while leaving those of frame transformations unchanged
because they are resolved in the body-attached rather than the reference frame. Introducing
transformation (8) in eq. (1) shows that the coordinates of frame transformations remain un-
affected, or invariant, under superimposed Euclidean transformations: p F q = S F −1 S
p ◦ Fq =
E −1 S −1 S
F p ◦ F E ◦ F E E F q = E F −1 E
p ◦ F q.

2.3 Derivatives of frame operations

Consider a parametric trajectory F (t), where t ∈ R denotes time. To simplify the notation, the
derivations presented in this section are carried out in the dual orthogonal matrix representation (6),
but similar results can be obtained with any representation. Because of the orthogonality, RT R =
I = I + ε 0, the time derivative of a dual orthogonal matrix takes the following form

dt R = RṼ (9)

where frame velocity Ṽ (t) = ω̃+ε ṽ is a dual skew-symmetric matrix made of the local coordinates of
the angular velocity ω̃(t) = RT dt R and of the linear velocity v(t) = RT dt u. In light of the discussion
in section 2.2, this frame velocity remains invariant under superimposed Euclidean transformations
and is thus referred to as the local frame velocity. Alternatively, the frame velocity can be resolved
in the inertial frame, dt R RT = RṼ RT , and it often referred to as the fixed-pole velocity [25].
It is not invariant under superimposed Euclidean transformations and is not used in the current
framework.

7
 Similarly, spatial derivatives and variations of frame fields over flexible bodies lead to local frame
derivatives,
dα R = RK̃, K̃ = κ̃ + ε k̃ (10)
˜
dδ R = RδΥ, ˜ = δθ
δΥ ˜ + ε δυ ˜ (11)
To ease their manipulation in the context of mechanics, coordinates of derivatives are associated
with six-dimensional arrays that stack the dual and primal parts, e.g. V = {v T , ω T }.
Because frame fields are continuous functions of space and time, derivative operators commute,
dt dα R = dα dt R, leading to the following compatibility relationships dt K̃ − dα Ṽ = K̃ Ṽ − Ṽ K̃ that
can be recast in terms of arrays
 
κ̃ k̃
dt K − dα V = K̂V = −V̂ K, K̂ = (12)
0 κ̃
Similar relationships hold for the other combinations of derivatives
ˆ V
dt δΥ − dδ V = δΥ (13a)
dδ K − dα Υ = K̂ δΥ (13b)

2.4 Subgroups
Subgroups2 of SE(3) provide a convenient framework for the description of frame transformations
in components with restricted kinematics such as lower-pair joints. For frames that belong to a
k-dimensional subgroup, frame velocities become
V = A V̄ (14)
where A is a 6 × k constant, full-rank matrix whose columns span the subalgebra of allowable
infinitesimal motions and V̄ is a k-dimensional array of velocity coordinates. Similarly, variations
take the form δΥ = A δΥ.¯
Subgroups for which V̂ A vanishes identically for any V̄ are called Abelian subgroups. Among
other interesting properties, their elements commute, F 1 ◦ F 2 = F 2 ◦ F 1 , and the commutativity
¯ = dδ V̄ , because A is full-ranked.
of derivative operators yields dt δΥ = dδ V or, equivalently, dt δΥ
Frame transformations for all lower-pair joints are elements of an Abelian subgroup of SE(3), with
the exception of spherical joints.

2.5 Parameterization of frame coordinates

The representations presented in section 2.1 involve additional or redundant parameters. In con-
trast, “non-redundant” or “minimal-set” representation involve six parameters only and are called
“parameterizations.” The frame coordinates become
F := P (15)
2
(H, ◦) is a subgroup of (G, ◦) if H is a subset of elements h of G such that (i) H is closed under ◦, i.e.,
h1 ◦ h2 = h3 ∈ H, (ii) the neutral element of H is that of G, i.e., h ◦ e = e ◦ h = h, and (iii) the inverse elements of
H are also in H, i.e., h−1 ∈ H.

8
Parameters P form six-dimensional arrays {pT , pTu }, where p and pu are three-dimensional arrays
of orientation/rotation parameters and position/displacement, respectively. The most common pa-
rameterizations fall under the category of “vectorial parameterizations” [6, 31, 32, 33], for which
p = f (φ)n and pu = T −T (p)u where T (p) is the tangent operator associated with the parameteri-
zation of orientation/rotation.
The coordinates of frames can also be expressed as dual array P = p̂ = f (φ̂)n̂ = p + ε pu where
f (φ̂) is a dual function of dual scalar φ̂ = φ + ε nT u and unit dual vector n̂ represents Chasles’
line [6, 29]. The Euler-Rodrigues and Cayley-Gibbs-Rodrigues parameterizations are used widely
and correspond to f (φ̂) = 2 sin(φ̂/2) and f (φ̂) = 2 tan(φ̂/2), respectively. The Euler-Rodrigues
parameterization is particularly simple to manipulate and will be used in the sequel.
Alternatively, three parameters can be obtained from a geometric description of rotation, such
as Euler angles, but while these are used widely in flight mechanics and astronomy, they are not
recommended for the analysis of flexible multibody systems.

Remark. In the kinematic description of classical multibody formulations, position/displacement

and orientation/rotation are viewed as two separate entities. Parameterizations consist of (1) orien-
tation/rotation parameters, p = f (φ)n, and (2) the position/displacement itself, u. In the motion
formalism, the parameters associated with the position/displacement are coupled to those of the
orientation/rotation because no parameterization exists for which T (p) = I, ∀p. For infinitesimal
values of the orientation/rotation parameters, T (p ≈ 0) ≈ I and the two approaches coincide:
pu ≈ u.
Unfortunately, minimal-set frame coordinates cannot be both global and singularity-free [11]:
configurations always arise where the coordinates grow without bound or are not defined uniquely.
While clever rescaling algorithms have been devised to overcome these problems [6], none is fully
satisfactory. The manipulation of minimal-set coordinates involves nonlinear operations and the
level of nonlinearity increases with the magnitude of the parameters. On the other hand, parame-
terizations are desirable because they form a minimal set, afford a clear geometric interpretation,
and can be related to physical measurements directly.
Clearly, representations and parameterizations both play an important role in the analysis of
multibody systems. When dealing with frame transformations involving relative motions of finite
magnitude, representations must be used. In contrast, when dealing with frame transformations
involving relative motions of moderate of infinitesimal magnitude, both representations and param-
eterizations can be used but typically, the latter are advantageous.
The mapping between a representation and a parameterization expressed as
p
F q = F(p P q ) ⇔ p P q = P(p F q ) (16)
p p
Operator F is the map that builds representation F q from parameterization P q and operator P is
the map that extracts parameterization p P q from representation p F q . Note that all parameterization
satisfy p P q = − q P p . The mapping expressed by eq. (16) should be performed only when the
parameterization is singularity-free.

2.5.1 Derivatives of parameterization

The time derivative of the parameters are related to the array of velocities as
dt P = T −1 (P )V (17)

9
where T (P ) is the tangent operator associated with the parameterization. For frame transforma-
tions, the derivatives of the map introduced in eq. (16) leads to dt p P q = P (p P q ) {V p ; V q } where

P (p P q ) = −T −1 (− p P q ) T −1 (p P q )
 
(18)

Similar relationships hold for other derivatives, e.g. dδ p P q = P (p P q ){δΥp ; δΥq }. Note that identity
p
T −1 (p P q ) = P̂ q + T −1 (− p P q ) holds for any parameterization.

2.5.2 Parameterization of subgroups

In practice, it is useful to relate frame transformations with a minimal set of coordinates. The
parameterization of a frame that belongs to a subgroup can be introduced as parameters α = A ᾱ
and the mapping from/to another representation in eq. (16) becomes
p
F q = F(A p ᾱq ) ⇔ p ᾱq = P̄(p F q ) (19)

where P̄ = (AT A)−1 AT P. The time derivative of these parameters takes the form

dt ᾱ = T̄ −1 (ᾱ)V̄ (20)

¯ Note that eq. (20) is only

where T̄ −1 (ᾱ) = (AT A)−1 AT T −1 (A ᾱ)A. Similarly, dδ ᾱ = T̄ −1 (ᾱ)δΥ.
usable within the range of validity of a selected parameterization but the choice f (φ̂) = φ̂ is always
well-defined for Abelian subgroups because T̄ −1 (ᾱ) = I, ∀ᾱ and gives rises to the simple governing
equations for the parameters: dt ᾱ = V̄ .

2.6 Spatial interpolation

The classical interpolation schemes [34, 35] used in the finite element method were developed to
interpolate displacement fields, which form a linear space. These schemes cannot be used directly
for interpolating nonlinear manifolds.
Figure 1 illustrates the typical interpolation problem for
ξ
four-node beam elements: the components of the frames at- E S
=k
E
F 0 F
= j
E
F E=FN
= k
tached at the nodes of the beam and the components of the =
E
interpolated field are denoted E F k , k = 0, . . . , N and E F ξ , re- ξ0 = -1 F
= ξ ξN = +1
spectively. Non-dimensional coordinate ξ ∈ [−1, 1] is defined S
= jk

along the axis of the beam.

Figure 1: Spatial interpolation
Numerous authors [36, 37, 38] have developed schemes for
along a beam
the interpolation of motion that can be written in the following
form
XN
hk (ξ) ξ S k = 0 (21)
k=0

where shape functions hk (ξ) are Lagrangian polynomials and ξ S k are a parameterization of frame
transformation ξ S k = E F −1 E ξ
ξ ◦ F k . Because the components of S k are resolved in the unknown
interpolated frame, eq. (21) is an implicit, nonlinear interpolation scheme, in contrast with its

10
explicit, linear counterpart used in classical finite element formulations. With the help of eq. (18),
the differentiation of (21) leads to the interpolated curvatures

N
!−1 N
X X
−1 ξ
K(ξ) = hi (ξ)T ( S i) dξ hk (ξ) ξ S k (22)
i=0 k=0

In general, interpolation scheme (21) must be solved numerically but a closed-form solution can
be found [39] for the Euler-Rodrigues parameterization [6] as outlined hereafter. First, given nodal
configurations F k , the relative transformations at the nodes are evaluated as j S k = E F −1 E
j ◦ F i , see
fig 1. Next, the corresponding Euler-Rodrigues parameters, j rk = P E.R. (j S k ) are computed together
p
with dual scalars j η k = 1 − j rTk j rk /4. Relative motion parameters ξ S k are then
N
ξ 1X
Sk = hj (ξ) j rk (23)
λ j=0

where dual scalar λ(ξ) is defined as

N
X N
X
λ2 (ξ) = 1 − 2 hj (ξ)hk (ξ)(1 − j η k ) (24)
j=0 k=j+1

Finally, the interpolated motion is obtained as E F ξ = E F k ◦ F T (k S ξ ). Any node k can be used for
the interpolation process.
Expression (24) involves a double summation but a single summation extending over the Nr =
N (N + 1)/2 independent relative motions at the nodes suffices
X
λ2 = 1 − 2 gβ (ξ)(1 − ηβ ) (25)
β

where gβ (ξ) = hj (ξ)hk (ξ). The following notation was introduced β • = N

P P PN
j=0 k=j+1 •, where the
summation extends for 0 ≤ β < Nr . Indices (j, k) = (0, 1) correspond to β = 0 and the subsequent
pairs are (0, 2), · · · (0, N ), (1, 2) · · · (1, N ), · · · , defining a unique correspondence between index β
and the associated values of indices j and k.
The above closed-form solution relies on the parameterization of the relative motions at the
nodes; while the motions at each node might be large, the relative motions at the nodes remain
small if the interpolation operation makes sense. In finite element applications involving structures
undergoing small deformations, the relative motions at the nodes remain very small and the present
approach provides a closed-form solution of the problem that is singularity free. Because it is
expressed in terms relative motions only, the interpolation scheme is objective [40]: the interpolated
motion remains unchanged under superimposed Euclidean transformations.
Finite element applications also require the interpolation of the curvature field. For the Euler-
ξ
Rodrigues parameterization, T −1 (ξ S i ) = ξ η i − S̃ i /2 and the first factor on the right-hand side of
ξ
eq. (22) becomes N ξ
= N ξ
P P
i=0 hi (ξ)( η i − S̃ i /2) q i=0 hi (ξ) η i = λ, where the second equality follows

from interpolation scheme (21) and ξ η i = 1 − ξ S Ti ξ S i /4.

11
 PN ξ
The interpolated curvature field reduces to K(ξ) = (1/λ) k=0 dξ hk (ξ) S k and eq. (23) leads
to K(ξ) = (1/λ2 ) β qβ (ξ)rβ , where qβ (ξ) = dξ hk (ξ)hj (ξ) − dξ hj (ξ)hk (ξ). Equation (25) show that
P
λ approaches unity for relative motions of small amplitude and furthermore, λ = 1 at the nodes.
Hence, under the small deformation approximations, the curvature field reduces to
X
K(ξ) = qβ (ξ)rβ (26)
β

Shape functions qβ (ξ) describe the curvature distribution for unit relative motions.
Next, with the help of eq. (18), increments in the relative motions at the
 nodes are related to
−1 −1
increments in the nodal motions as ∆rβ = −T (rβ )∆Υj + r̃β + T (rβ ) ∆Υk . Recasting this
relationship is matrix form yields
∆ř = P̌ (ř)∆Υ̂ (27)
where array řT = {rT0 , rT1 , · · · , rTN } stores the Nr incremental relative motion parameters at the
nodes and array ∆Υ̂T = {∆ΥT0 , ∆ΥT1 , · · · , ∆ΥTN } stores the incremental nodal motions. Line β
of matrix P contains two non-vanishing entries only: if rβ = rjk , these entries are −T −1 (rβ ) and
r̃β + T −1 (rβ ) in columns j and k, respectively. Matrix P depends on the relative motions at the
nodes only.

3 Multibody systems
Within the framework presented above, multibody systems are described as a collection of inter-
acting body-attached frames. The frames that can undergo transformations of arbitrary amplitude
are kept as explicit unknown variables. They include the body-attached frames as well as the frame
transformations in kinematic joints that belong to subgroups of SE(3). Frame transformations
of moderate amplitudes, like those used for the deformation measures of flexible bodies, are not
explicit unknown variables but are evaluated as parameterizations of the relative motions between
body-attached frames.
The current framework falls in the category of mixed-variable formulations because it involves
relative variables (frame transformations) and absolute variables (body-attached frames). This
description results in a non-minimal set of variables because body-attached frames and frames
transformations are related to each other. Each of the M body-attached frames has 6 degrees
of freedom (dofs) and each of the m frame transformations i adds ki < 6 dofs. For each frame
transformation, six Lagrange multipliers are added to enforce the constraints imposed on the cor-
responding body-attached frames. The system is thus described by set {F } of M + m frames
involving N = 6M + m
P
i=1 i dofs and a set of 6m Lagrange multipliers {λ}. The system velocities
k
{V } and accelerations {V̇ } are N -dimensional arrays of coordinates made of six coordinates for
each body-attached frames and ki coordinates for frame transformation i.
In the finite element approach, each element e of the Ne elements in the system involves subset
{•}e of the system variables.

3.1 Equations of motion

Hamilton’s principle requires the variation of the action integral to be stationary for the actual
trajectory of the system between time ti and tf . Lagrange’s multipliers method is used to enforce

12
holonomic constraints and hence, the augmented action integral is used. In the finite element
approach, the action integral is the sum of the contributions of each element
Ne
X Z tf
K e − Ve − {λ}Te {C}e dt = 0
 
δ (28)
e=1 ti

where, for each element e, K e ({F }e , {V }e ) is the kinetic energy, Ve ({F }e ) is the potential en-
ergy associated with internal and external forces, and λe are the Lagrange multipliers that enforce
kinematic constraints {C}e ({F }e ). Taking the variation leads to
Z tf
{δΥ}Te {f }e + {B}Te {λ}e + δ{λ}Te {C}e dt = 0
 
(29)
ti

where {f }e ({F }e , {V }e , {V̇ }e ) are the element forces and {B}e ({F }e ) is the gradient of con-
straint {C}e such that δ{C}e = {B}e {δΥ}e . The resulting equations of motion are second-order,
differential-algebraic equations

{Ḟ } = {F }{Ṽ } (30a)

{f } + {B}T {λ} = 0 (30b)
{C} = 0 (30c)

3.2 Time integration method

The equations of motion can be integrated in time with the generalized-α scheme [41]. Other
schemes [42, 43, 44] are also applicable. Denoting the time step by n, the discretized equations of
motion are

{F |n+1 } = {F |n } ◦ {F(w)} (31a)

T
{r|n+1 } = {f |n+1 } + {B|n+1 } {λ|n+1 } = 0 (31b)
{C|n+1 } = 0 (31c)

together with integration formulæ3 that relate {w}, {V |n } and {V̇ |n } to {V |n+1 } and {V̇ |n+1 }.
Kinematic update (31a) states that each frame at step n + 1 results from an increment from its
value at step n by a parametrized frame transformation, where {w} is a N-dimensional array of
frame parameters. The construction of these incremental frame transformations from the parameters
follows eq. (16) for the body-attached frames and eq. (19) for the frame transformations.
3
The integration formulas of the generalized-α [41] are

{w} = h{V |n } + (0.5 − β)h2 {a|n } + βh2 {a|n+1 }

{V |n+1 } = {V |n } + (1 − γ)h{a|n } + γh{a|n+1 }
 
{a|n+1 } = (1 − αf ){V̇ |n+1 } + αf {V̇ |n } − αm {a|n } /(1 − αm )

where a is an array of algorithmic accelerations and h is the time step size. The algorithm parameters are selected
to achieve a desired spectral radius ρ ∈ [0, 1): αm = (2ρ − 1)/(ρ + 1), αf = ρ/(ρ + 1), γ = 0.5(3 − ρ)/(ρ + 1) and
β = 1/(1 + ρ)2 .

13
 Equations (31) are nonlinear and solved using a Newton iteration procedure. Lineariza-
tion of eq. (31b) yields ∆{r} = {M }{∆V̇ } + {G}{∆V } + {K}{∆Υ} + {B}T {λ}, where
{M ({F })}, {G}({F }, {V }) and {K}({F }, {V }, {V̇ }) are the N × N tangent mass, damp-
ing/gyroscopic/centrifugal, and stiffness matrices, respectively. In the finite element approach,
these matrices are typically sparse. Corrections to the unknown variables are computed by solving
linear system
{r}∗
   
{∆n}
ST =− (32)
{∆λ} {C}∗
where {r}∗ and {C}∗ denote the evaluation of the left-hand side of eqs. (31b–31c) at the current
iteration of the Newton procedure and S T ({M }, {G}, {K}, {∆n}, {B}) is the iteration matrix of
the generalized-α scheme4 .

4 Typical Elements
This section provides a succinct description of the various elements found in typical multibody
codes. The aim is to demonstrate the ability of the motion formalism to simplify the formulation
of the various elements and to highlight the desirable properties of the resulting equations.

4.1 The rigid body

A rigid body is modeled as a single frame attached to the body and its associated kinetic energy is
Krgb = {V }Trgb {M }rgb {V }rgb /2, where {M }rgb denotes the 6 × 6 constant mass matrix of the body.
To ease the notation of this section, the reference to the element has been dropped. With the help
of eq. (13a), the virtual work done by the inertial forces
is obtained by variation of the kinetic
energy: dδ K = dδ V P = (dt δΥ + Ṽ δΥ) P = dt δΥ P − δΥT (M V̇ − V̂ T P ), where P = M V are
T T T

the local coordinates of the momentum vector. The inertial forces of the rigid body become

{f }rgb (V , V̇ ) = M V̇ − V̂ T P (34)

and their linearization  

∆{f }rgb = M ∆V̇ + −V̂ T M − P̆ ∆V (35)

where P̆ is defined implicitly by identity P̆ ∆V = V̂ T ∆P . These forces exhibit a second-order

nonlinearity term in the velocity, which accounts for gyroscopic effects.
The inertial forces depend on velocity and acceleration components, as expected, but remain
independent of the position and orientation variables that define the configuration of the body with
respect to the inertial frame. Such formulations are referred to as intrinsic [45, 46]. The same
observations carry over to the linearization.
4
With the linearization of the integration formulas, the iteration matrix of the generalized-α scheme takes the
following form  0
β {M } + γ 0 {G} + {K}{T (w)} {B}T

ST = (33)
{B}{T (w)} 0
where β 0 = (1 − αm )/(βh2 (1 − αf ) and γ 0 = γ/(βh). Note that the use of the local parametrization in eq. (31a) leads
to the tangent operator of parameterization, i.e., {∆Υ} = {T (w)}{∆w}.

14
4.2 The flexible joint
The flexible joint, sometimes called “bushing element” or “force element,” is a common element
found in many multibody systems: it consists of an elastic medium connected to two frames of
the system, called “’handles.” The relative motion between these two handles deforms the elastic
medium, generating forces and moments at the handles. The elastic medium can be thought of as a
set of six concentrated springs, three rectilinear and three torsional springs, connecting two bodies,
but in the general case, the linear and torsional deformations couple elastically.
From eq. (1), the relative motion at the joint is given by p F q where p and q are the handle-
attached frames. The joint deformation measures are selected to be the components of a vectorial
parameterization of the relative motion as introduced in eq. (16)

E({F }flx ) = P(F −1

q ◦ F p) (36)

where {F }flx = {F p , F q }. This selection guarantees that the deformation measures are both objec-
tive and tensorial.
The behavior of the elastic medium is modeled by a deformation energy function, V flx ({F }flx ) =
(E − E 0 )T D (E − E 0 )/2, where symmetric stiffness matrix D, of size 6 × 6 is, in general, fully
populated and E 0 characterizes the undeformed configuration of the joint. Variation of the potential
is dδ V flx = dδ E T f , where the forces in the joint are f (E) = D (E − E 0 ). In view of eq. (18), the
variation of the potential is now recast as dδ V flx = {δΥ}Tflx {f }flx , where the generalized forces at the
handles are
{f }flx ({F }flx ) = P T (E) D (E − E 0 ) (37)
Their linearization ∆{f }flx = {K}f lx {∆Υ}flx leads to tangent stiffness matrix

{K}f lx ({F }flx ) = P T (E)D + H T (f , E) P (E)

 
(38)

where matrix H is defined implicitly by identity H T (f , E)∆E = ∆P T (E)f (E). The generalized
forces depend on the configuration of the two handles only through the deformation measures, which
are functions of their relative motions. Consequently, the forces are invariant under superimposed
Euclidean transformations. The same observations carry over to the linearization. As was observed
for the rigid body, the formulation is intrinsic.
The flexible joint illustrates the clear need for both representations and parameterizations. The
description of the configurations of the two handles requires representations because motions of
arbitrary amplitudes are involved. On the other hand, it is reasonable to assume that the defor-
mation of the joint is moderate: indeed, the flexible joint cannot fold onto itself. The use of a
parameterization for the deformation measures then seems justified.
While parameterizations form convenient choice of deformation measures, a particular selection
must be made. Let E ER = P ER (F −1 −1
q ◦F p ) and E CGR = P CGR (F q ◦F p ) be the deformation measures
based on the Euler-Rodrigues and Cayley-Gibbs-Rodrigues parameterizations, respectively. If two
simulations are performed based on these two deformation measures, different predictions will be
obtained. This is not a flaw of the formulation. Rather, it stems from the fact that different
constitutive laws have been used. Indeed, if the same joint stiffness matrix is used in the two
simulations, the forces applied to the handles, f (E ER ) and f (E CGR ) will differ. Clearly, joint
stiffness matrix D must be adapted to the deformation measure used in the analysis. If obtained

15
experimentally, the deformations measures used to reduce the measured data must match those
used in the analysis.
The parameterizations introduced in section 2.5 remain invariant under a change of material
frame [47, 48, 49]. Consider two other body-attached frames to the joint F a = F p ◦ p F a and
F b = F q ◦ q F b where p F a = q F b = F ? characterizes the change of material frame. Let force and
moment components resolved in frames F a and F b be f ? = C T? f where C ? , the representation
in eq. (3) of F ? , includes a change of orientation and a lever-arm effect. On the other hand, the
deformation measures are E ? = P(F −1 a ◦F b ). Because the deformation energy must remain invariant
under a change of material frame E T? f ? = E T f , which implies that the deformation measures
transform as E ? = C −1 ? E. All vectorial parameterizations of motion satisfy this relationship [32],
thus providing consistent deformation measures.
Classical deformation measures, such as those used in commercial finite element codes, are not
invariant under frame transformation and hence, do not result in invariance of the deformation
energy, a fundamental requirement for the soundness of the formulation.
Clearly, deformation measures must be both objective and tensorial. Because the classical de-
formation measures meet these requirements for infinitesimal deformations only, formulations of
flexible joints based on these deformation measures can only describe their behavior under infinites-
imal deformation conditions. In contrast, because motion-based deformation measures always meet
these requirements, formulations of flexible joints based on these deformation measures remain
meaningful under finite deformations.

4.3 Kinematic joints

Kinematic joints enforce a frame transformation between two body-attached frames that is charac-
terized by kJ < 6 relative dofs. For instance, a revolute joint has kJ = 1 relative dof and a cylindrical
joint has kJ = 2 dofs. This situation is mathematically described as p F q = F J ⇔ F −1 p
J ◦ Fq = I
where p F q , defined in eq. (1), is the frame transformation between body-attached frames p and q
and F J is the particular frame transformation that joint J prescribes. In the current framework, F J
is treated as an unknown variable. This kinematic relationship is conveniently expressed through a
parameterization [50], leading to constraint equation

{C}jnt {F }jnt = P F −1 −1



J ◦ Fp ◦ Fq = 0 (39)

where {F }jnt = {F p , F q , F J }. Its linearization evaluated at constraint satisfaction, {C}jnt = 0,

yields constraint gradient
{B}joint {F }jnt = −C −1

 
J I −A J (40)
where C J is the representation in eq. (3) of F J and AJ is the matrix that characterizes the subgroup
to which joint J restricts the motion, see section 2.4. The constraint equations and their gradients
depend on local frame transformations but not the current position and orientation of the joint
with respect to an inertial frame; hence, it is an intrinsic formulation.
The use of relative variables streamlines the formulation of joint features such as lumped spring-
dampers, controller and sensors within joints. The virtual work contribution from such components,
¯ T f jnt , is evaluated easily as
δΥ

˙ = T̄ −T (ᾱ) K(ᾱ − ᾱ0 ) + D ᾱ˙ + τ̄

 
f jnt (ᾱ, ᾱ) (41)

16
where K and D are kJ × kJ stiffness and damping matrices, ᾱ0 is an kJ -dimensional array of the
natural lengths associated with the elastic behavior, and τ̄ is an kJ -dimensional array of forces and
torques exerted by an actuator.

4.4 Flexible beams

A beam is described as frame field F (α, t), where α is the spatial parameter along the reference line
of the beam. The time derivative of the frame field yields the velocity field, see eq. (9), whereas
the spatial derivative introduced in eq. (10) yields the array of curvatures K and K 0 in the current
and undeformed configurations, respectively.
The strain energy density of the beam is a quadratic function of the curvatures, V (α, t) =
(K − K 0 )T DB (K − K 0 )/2 and its kinetic energy density is a quadratic function of the velocities,
T
K(α, t) = V M B V/2, where D B (α) and M B (α) are the 6 × 6 sectional stiffness and mass matrices,
respectively, computed with respect to the reference axis of the beam.
The beam inertial forces are obtained by following the procedure described in sec-
tion 4.1 for the rigid body to find the nodal inertial forces as {f }beam ({V }beam , {V̇ }beam ) =
{M }beam {V̇ }Rbeam − {Ṽ T }beam {M }beam {V }beam . The constant mass matrix of the element is
{M }beam = LT (α)M B (α)L(α) dα, where matrix L(α) stores the shape functions for the velocity
field. Linearization of the inertial forces follows the procedure described in section 4.1. For simplic-
ity, interpolation matrix L(α) is evaluated in the reference configuration of the beam, resulting in
the constant mass matrix.
Using curvature interpolation scheme (26), the curvature field can be recast as E(α) = Q(α)ř,
where matrix Q(α) stores shape functions qβ (α) and array ř stores the parameterization of the
Nr relative motions at the nodes. Variation of the curvature field now becomes dδ E = Qdδ ř =
QP (ř)dδ Υ̂, where the second equality results from eq. (27).
Variation of the strain energy becomes δ V = dδ E T DB E and introducing the above discretization
leads to δ V = dδ Υ̂T P T QT DB Q ř. Integration over the element then yields the nodal elastic forces

fˆebeam (ř) = P T (ř)Ď ř (42)

R
where the constant stiffness matrix of the element is Ď = QT DB Q dα. Linearization of the
elastic forces yields ∆fˆe = K ∆Υ̂, where the effective stiffness matrix of the element is
beam
h i
K(ř) = P T (ř)D̂ + H T (fˇ, ř) P (ř) (43)

where matrix H T is defined implicitly by identity H T (fˇ, ř)∆ř = ∆P T (ř)fˇ and fˇ = Ď ř.
The nodal elastic forces depend on the relative motions at the nodes only, not on the configura-
tion of the beam respect to an inertial frame: the formulation is intrinsic [46, 51]. These remarks
also apply to the linearization.
Because the interpolation scheme couples displacement and rotation variables the shear-looking
phenomenon is circumvented [37, 39, 52]. Because the mass and stiffness matrices of the element
remain constant, they can be evaluated once only at the onset of the simulation, resulting in
improved computational efficiency.

17
 The kinetic energy in the beam is treated in the same manner as that of a rigid body with
multiple nodes. Identical nodal inertial forces are obtained for two elements. Similarly, the strain
energy in the beam is treated in the same manner as that of a flexible joint with multiple handles.
Identical nodal elastic forces are obtained for the two elements, compare eqs. (37) and (42), and
identical effective stiffness matrices result, compare eqs. (38) and (43). Interpolation tools and
the finite element method are used to derive the constant mass and stiffness matrices of the beam
element whereas the mass matrix of the rigid body and the stiffness matrix of the flexible joint are
assumed to be given or obtained experimentally.
This remarkable similitude between the various elements stems from the unified kinematics
description presented in this paper: configurations are described by body-attached frames, and
velocities and deformations by derivatives of frame transformations, all resolved in material frames.
Modal [53] and shell elements [39] can be formulated in a similar manner, leading to the same
desirable features.

5 Conclusions
A finite element approach for the analysis of flexible multibody systems was presented. It is based on
the motion formalism that (1) uses configuration and motion to describe the kinematics of flexible
multibody systems, (2) recognizes that these are members of the Special Euclidean group thereby
coupling their displacement and rotation components, and (3) resolves all tensors components in
local frames.
The following advantages result. (1) A clear distinction was drawn between body-attached
frames, whose coordinates are inertial, and frame transformations, whose coordinates are local. (2)
Large amplitude motions are handled properly, eliminating the potential occurrence of singularities.
(3) The formulation is intrinsic: position and rotation variables describing the configuration of struc-
tural components with respect to an inertial system are eliminated from the governing equations.
(4) Deformation measures resulting from the systematic use of the formalism are objective and
tensorial. (5) The mass and stiffness matrices of the flexible structural elements remain constant
throughout the simulation, resulting in improved computational efficiency. (6) Parameterizations
are used for local transformations only. (7) Equilibrium equations are invariant under superim-
posed Euclidean transformations because the large amplitude motions are filtered out, leading to a
reduced order of nonlinearity.

References
[1] W. O. Schiehlen. Multibody system dynamics: Roots and perspectives. Multibody System
Dynamics, 1:149–188, 1997.

[2] A. A. Shabana. Dynamics of Multibody Systems. Cambridge University Press, Cambridge,

third edition, 2005.

[3] E. Cosserat and F. Cosserat. Théorie des corps déformables. A. Hermann et Fils, Paris, second
edition, 1909.

18
 [4] P. M. Naghdi. The theory of plates and shells. In S. Flügge, editor, Handbuch der Physik,
volume 2, pages 425–640. Springer-Verlag, Berlin, 1972.

[5] M. Géradin and A. Cardona. Flexible Multibody System: A Finite Element Approach. John
Wiley & Sons, New York, 2001.

[6] O. A. Bauchau. Flexible Multibody Dynamics, volume 176 of Solid Mechanics and Its Applica-
tions. Springer Netherlands, Dordrecht, 2011.

[7] O. A. Bauchau and J. I. Craig. Structural Analysis with Applications to Aerospace Structures.
Springer, Dordrecht, 2009.

[8] E. Reissner. The effect of transverse shear deformation on the bending of elastic plates.
Zeitschrift für angewandte Mathematik und Physik, 12:A.69–A.77, 1945.

[9] E. Reissner. On bending of elastic plates. Quarterly of Applied Mathematics, 5:55–68, 1947.

[10] R. D. Mindlin. Influence of rotatory inertia and shear on flexural motions of isotropic elastic
plates. Journal of Applied Mechanics, 18:31–38, 1951.

[11] J. Stuelpnagel. On the parameterization of the three-dimensional rotation group. SIAM Review,
6(4):422–430, 1964.

[12] R. S. Ball. A Treatise on the Theory of Screws. Cambridge University Press, Cambridge, 1998.

[13] O. Bottema and B. Roth. Theoretical Kinematics. Dover Publications, Inc., New York, 1979.

[14] J. M. McCarthy. An Introduction to Theoretical Kinematics. The MIT Press, Cambridge, MA,
1990.

[15] R. M. Murray, Z. Li, and S. S. Sastry. A Mathematical Introduction to Robotic Manipulation.

CRC Press, Boca Raton, FL, USA, 1994.

[16] J. M. Selig. Geometric Fundamentals of Robotics. Monographs in computer science. Springer,

New York, 2005.

[17] J. Angeles. Fundamentals of Robotic Mechanical Systems. Theory, Methods, and Algorithms,
volume 124 of Mechanical Engineering Series. Springer International Publishing, New York,
2014.

[18] A. Müller and Z. Terze. The significance of the configuration space Lie group for the constraint
satisfaction in numerical time integration of multibody systems. Mechanism and Machine
Theory, 82:173–202, 2014.

[19] A. Müller. A note on the motion representation and configuration update in time stepping
schemes for the constrained rigid body. BIT Numerical Mathematics, 56(3):995–1015, 2016.

[20] A. Müller. Screw and Lie group theory in multibody kinematics. Multibody System Dynamics,
43:37–70, 2018.

19
[21] V. Giavotto, M. Borri, P. Mantegazza, G. Ghiringhelli, V. Carmaschi, G. C. Maffioli, and
F. Mussi. Anisotropic beam theory and applications. Computers & Structures, 16(1-4):403–
413, 1983.

[22] O. A. Bauchau and S. L. Han. Three-dimensional beam theory for flexible multibody dynamics.
Journal of Computational and Nonlinear Dynamics, 9(4):041011 (12 pages), 2014.

[23] J. Rooney. A comparison of representations of general spatial screw displacement. Environment

and Planning B: Planning and Design, 5(1):45–88, 1978.

[24] D. P. Chevallier. Lie algebra, modules, dual quaternions and algebraic methods in kinematics.
Mechanism and Machine Theory, 26(6):613–627, 1991.

[25] M. Borri, L. Trainelli, and C. L. Bottasso. On representations and parameterizations of motion.

Multibody Systems Dynamics, 4:129–193, 2000.

[26] J. Funda, R.H. Taylor, and R.P. Paul. On homogeneous transforms, quaternions, and compu-
tational efficiency. IEEE Transactions on Robotics and Automation, 6(3):382–388, 1990.

[27] J. Funda and R.P. Paul. A computational analysis of screw transformations in robotics. IEEE
Transactions on Robotics and Automation, 6(3):348–356, 1990.

[28] R. A. Wehage. Quaternions and Euler parameters. A brief exposition. In E.J. Haug, editor,
Computer Aided Analysis and Optimization of Mechanical Systems Dynamics, pages 147–180.
Springer-Verlag, Berlin, Heidelberg, 1984.

[29] J. Angeles. The application of dual algebra to kinematic analysis. In J. Angeles and E. Za-
khariev, editors, Computational Methods in Mechanical Systems, volume 161, pages 3–31.
Springer-Verlag, Heidelberg, 1998.

[30] E. Pennestrı̀ and R. Stefanelli. Linear algebra and numerical algorithms using dual numbers.
Multibody System Dynamics, 18:323–344, 2007.

[31] S. L. Han and O. A. Bauchau. Manipulation of motion via dual entities. Nonlinear Dynamics,
85(1):509–524, July 2016.

[32] O. A. Bauchau and J. Y. Choi. The vector parameterization of motion. Nonlinear Dynamics,
33(2):165–188, 2003.

[33] O. A. Bauchau and L. Trainelli. The vectorial parameterization of rotation. Nonlinear Dy-
namics, 32(1):71–92, 2003.

[34] T. J. R. Hughes. The Finite Element Method. Prentice Hall, Inc., Englewood Cliffs, New
Jersey, 1987.

[35] K. J. Bathe. Finite Element Procedures. Prentice Hall, Inc., Englewood Cliffs, New Jersey,
1996.

20
[36] T. Merlini and M. Morandini. The helicoidal modeling in computational finite elasticity. Part I:
Variational formulation. International Journal of Solids and Structures, 41(18-19):5351–5381,
2004.

[37] V. Sonneville, O. Brüls, and O. A. Bauchau. Interpolation schemes for geometrically ex-
act beams: a motion approach. International Journal of Numerical Methods in Engineering,
112(9):1129–1153, 2017.

[38] S. L. Han and O. A. Bauchau. On the global interpolation of motion. Computer Methods in
Applied Mechanics and Engineering, 337(10):352–386, 2018.

[39] O. A. Bauchau and V. Sonneville. Formulation of shell elements based on the motion formalism.
Applied Mechanics, 2(4):1009–1036, December 2021.

[40] L. E. Malvern. Introduction to the Mechanics of a Continuous Medium. Prentice Hall, Inc.,
Englewood Cliffs, New Jersey, 1969.

[41] O. Brüls, A. Cardona, and M. Arnold. Lie group generalized-α time integration of constrained
flexible multibody systems. Mechanism and Machine Theory, 48:121–137, February 2012.

[42] Z. Terze, A. Müller, and D. Zlatar. Lie-group integration method for constrained multibody
systems in state space. Multibody System Dynamics, 34(3):275–305, 2015.

[43] S. Hante and M. Arnold. RATTLie: A variational Lie group integration scheme for constrained
mechanical systems. Journal of Computational and Applied Mathematics, 387:112492, May
2021.

[44] V. Wieloch and M. Arnold. BDF integrators for constrained mechanical systems on Lie groups.
Journal of Computational and Applied Mathematics, 387:112517, May 2021.

[45] A. E. H. Love. A Treatise on The Mathematical Theory of Elasticity. Dover, New York, fourth
edition, 1944.

[46] D. H. Hodges. Geometrically exact, intrinsic theory for dynamics of curved and twisted
anisotropic beams. AIAA Journal, 41(6):1131–1137, June 2003.

[47] O. A. Bauchau, L. H. Li, P. Masarati, and M. Morandini. Tensorial deformation measures for
flexible joints. Journal of Computational and Nonlinear Dynamics, 6(3):0310021–8, 2011.

[48] O. A. Bauchau and S. L. Han. Flexible joints in structural and multibody dynamics. Mechanical
Sciences, 4(1):65–77, 2013.

[49] P. Masarati and M. Morandini. Intrinsic deformable joints. Multibody System Dynamics,
23(4):361–386, 2010.

[50] V. Sonneville and O. Brüls. A formulation on the special Euclidean group for dynamic analysis
of multibody systems. Journal of Computational and Nonlinear Dynamics, 9(4):041002, 2014.

21
 [51] M. Borri and C. L. Bottasso. An intrinsic beam model based on a helicoidal approximation.
Part I: Formulation. International Journal for Numerical Methods in Engineering, 37:2267–
2289, 1994.

[52] V. Sonneville, A. Cardona, and O. Brüls. Geometrically exact beam finite element formu-
lated on the special Euclidean group SE(3). Computer Methods in Applied Mechanics and
Engineering, 268(1):451–474, 2014.

[53] V. Sonneville, O. Brüls, and O. A. Bauchau. Modal reduction procedures for flexible multibody
dynamics. Multibody System Dynamics, 51(4):377–418, April 2021.

22

View publication stats