https://doi.org/10.1007/s11044-017-9606-3
On the use of absolute interface coordinates
in the floating frame of reference formulation for flexible
multibody dynamics
Marcel Ellenbroek1 · Jurnan Schilder1
Received: 3 May 2017 / Accepted: 28 November 2017
© The Author(s) 2017. This article is published with open access at Springerlink.com
Abstract In this work a new formulation for flexible multibody systems is presented based on the floating frame formulation. In this method, the absolute interface coordinates are used as degrees of freedom. To this end, a coordinate transformation is established from the abso-lute floating frame coordinates and the local interface coordinates to the absoabso-lute interface coordinates. This is done by assuming linear theory of elasticity for a body’s local elastic deformation and by using the Craig–Bampton interface modes as local shape functions. In order to put this new method into perspective, relevant relations between inertial frame, coro-tational frame and floating frame formulations are explained. As such, this work provides a clear overview of how these three well-known and apparently different flexible multibody methods are related. An advantage of the method presented in this work is that the result-ing equations of motion are of the differential rather than the differential-algebraic type. At the same time, it is possible to use well-developed model order reduction techniques on the flexible bodies locally. Hence, the method can be employed to construct superelements from arbitrarily shaped three dimensional elastic bodies, which can be used in a flexible multi-body dynamics simulation. The method is validated by simulating the static and dynamic behavior of a number of flexible systems.
Keywords Flexible multibody dynamics· Floating frame formulation · Corotational frame formulation· Inertial frame formulation · Craig–Bampton method · Model order reduction
1 Introduction
Flexible multibody dynamics is concerned with the study of machines and mechanisms that consist of multiple deformable bodies. Although the elastic deformation within a single body
B
J. Schilderj.p.schilder@utwente.nl M. Ellenbroek
m.h.m.ellenbroek@utwente.nl
1 Faculty of Engineering Technology, University of Twente, P.O. Box 217, 7500AE Enschede,
Fig. 1 Graphical representation of a single flexible body using three different formulations. The inertial frame formulation is expressed in the absolute nodal coordinates (left). The corotational formulation uses a corotational frame for each element (center). The floating frame formulation uses a floating frame for each body (right). Figure was made using InkScape
can often be considered as small, the large rigid body rotations between different bodies make the problem to be of a geometrically nonlinear nature.
The kinematics of a body can be described by the motion of a set of coordinate frames, each frame being rigidly attached to a material point on the body. Connections between bodies are introduced by kinematic constraints. These relate the motions of the coordinate frames of the so-called interface points located at different bodies.
Three essentially different commonly used descriptions are available for flexible multi-body systems: the inertial frame formulations, the corotational frame formulations and the floating frame formulations. These formulations are significantly different in the way they describe a body’s elastic behavior. Hence, there are important differences in the choice of the degrees of freedom and consequently in the way kinematic constraints between bodies are enforced. Figure1gives a graphical interpretation of each formulation. A comprehensive overview of these formulations, their background and particularities can be found in [1].
The inertial frame formulation is essentially a nonlinear finite element formulation of which the details can be found in standard textbooks, such as [2]. Elastic deformations are described using the nonlinear Green–Lagrange strain definition, which describes large rigid body rotations correctly. In this formulation, each body is discretized using a set of global interpolation functions. The degrees of freedom are the absolute positions and orientations of the coordinate frames located at the nodes of the finite element mesh. Constraints between bodies can be simply enforced by directly equating the appropriate degrees of freedom of the nodes that both bodies have in common at a certain interface point. For applications to flexible multibody systems, the geometrically exact beam theory is a commonly used example of an inertial frame formulation [3]. However, difficulties in frame invariance may arise due to the interpolation of the large rotations [4,5].
The corotational frame formulation can be interpreted as the geometrically nonlinear ex-tension of a linear finite element formulation [6,7]. A corotational frame describes the large rigid body motion of an element with respect to the inertial frame. Small elastic deforma-tions within the element are superimposed using the linear finite element matrices, based on the Cauchy strain definition. The nonlinear finite element model is obtained from the linear finite element model by simply pre- and post-multiplying the mass and stiffness matrices with the rotation matrices corresponding to the corotational frames. Note that the formula-tion is still able to geometrically describe geometrical nonlinear elastic deformaformula-tions within a single body due to the fact that each element is given its own corotational frame. That is, provided that the small strain assumption within an element holds.
Similar to the inertial frame formulation, the corotational frame formulation uses the ab-solute nodal coordinates as degrees of freedom. Consequently, both methods also enforce constrains in the same way. However, in order to arrive at this formulation, a unique
kine-matic relation must be established that expresses the coordinates of an element’s corota-tional frame in terms of the element’s absolute nodal coordinates. For a variety of standard elements, this is commonly done by locating the corotational frame in one of the element’s nodes or by expressing the corotational frame coordinates as a certain weighted average of the element’s nodal coordinates. Alternatively, it is possible to define the orientation of each corotational frame by demanding zero elastic deformation at its location [8]. A drawback of standard corotational frame formulations is that they do not distinguish between rigid and flexible bodies. That is, rigid bodies are modeled as flexible bodies with a large stiffness, resulting in a less efficient formulation for multibody systems that contain both rigid and flexible bodies.
The floating frame formulation can be interpreted as the flexible extension of a rigid multibody formulation. A floating frame describes the large rigid body motion of a body with respect to the inertial frame. When linear theory of elasticity is used to describe the flexible behavior locally, mass and stiffness matrices are obtained from a body’s linear finite element model. These system matrices can be reduced using a wide variety of comprehen-sive well-developed model order reduction techniques, such as the Craig–Bampton method [9,10].
In the floating frame formulation, the degrees of freedom thus consist of the absolute co-ordinates of the floating frame and a set of local generalized coco-ordinates used to describe the body’s local linear elastic behavior. Because the absolute interface coordinates are not part of the degrees of freedom, the kinematic constraints are typically highly nonlinear equations in terms of many degrees of freedom. As no analytical solution of these equations might be found, the constraints are commonly enforced using Lagrange multipliers, increasing the total number of unknowns. The resulting equations of motion are combined differential-algebraic equations [11,12].
The possibility of using well-developed model order reduction techniques to reduce com-putational cost makes the floating frame formulations the preferred formulation in many situations in which the elastic deformation within a body can be considered small. The dis-advantage of this formulation in satisfying kinematic constraints could be eliminated if it is possible to express the floating frame coordinates and the local elastic deformation directly in terms of the interface coordinates. The search for such a kinematic coordinate transfor-mation has led to the development of superelements in flexible multibody formulations. Methods have been developed in which the floating frame is located at an interface point [13] or expressed as the weighted average of the interface coordinates [14]. The first op-tion introduces an unwanted discriminaop-tion between the interface points, which makes the results dependent on the interface point chosen. Moreover, in general a better accuracy is obtained when the floating frame is close to the body’s center of mass. In the second option the motion of the floating frame cannot be interpreted as the motion of a material point, but only as the body’s average rigid body motion.
In this paper, a new method is presented, with which it is possible to obtain a floating frame formulation in terms of interface coordinates only, that does not suffer from the dis-advantages mentioned above. The method offers the combined advantage of being able to enforce constraints without the use of Lagrange multipliers and still have the possibility to use Craig–Bampton based model order reduction techniques on the bodies’ linear finite element models. In fact, the method demonstrates the interchangeability of the absolute in-terface coordinates and the combination of the floating frame coordinates and local elastic coordinates corresponding to these Craig–Bampton modes. It will be shown that in the dy-namic equations local coordinates can be substituted by global coordinates, and vice versa.
Essential for the method presented here is the fact that the Craig–Bampton modes are able to describe rigid body motions. In the floating frame formulation, these rigid body modes
must be eliminated in order to describe the system’s motion uniquely. However, in this work this property of the Craig–Bampton modes is used to eliminate the floating frame coordi-nates and express the local elastic degrees of freedom solely in terms of the global motion of the interface points. This is done by demanding that the elastic body has no deformation at the location of the floating frame. This requirement is met without the need to locate the floating frame in an interface point.
The subsequent sections are organized as follows: In Sect.2, the kinematic description of a coordinate system attached to a material point of a flexible body is introduced using a position vector and a rotation matrix. In Sect.3, the kinematics of the floating frame for-mulation are introduced. A relation is developed that expresses the local coordinates of a material point on a flexible body as the difference between the absolute coordinates of that material point and the floating frame coordinates. Section 4 introduces the static Craig– Bampton modes to describe a body’s local linear elastic displacement field. It is shown that the generalized coordinates corresponding to these modes can be expressed in terms of the absolute interface coordinates and the floating frame coordinates. In Sect.5, the floating frame coordinates are expressed in terms of the absolute interface coordinates by demand-ing zero deformation at the location of the floatdemand-ing frame. In Sect.6, the obtained coordinate transformations are applied to the standard equation of motion of the floating frame formu-lation. At this point, the kinematic constraint equations are applied without using Lagrange multipliers. The resulting equations of motion can be solved incrementally using numerical time integration. Section7demonstrates the new method on a number of test cases. The most important conclusions finalize the paper.
2 Kinematics of a material point on a flexible body
Consider a flexible body with a material point Pj. A coordinate frame Ej is rigidly attached
to this material point such that the pair{Pj, Ej} defines a Euclidean coordinate system. As
{Pj, Ej} defines both the location and orientation of the frame attached to Pj, the pair will
be referred to as the generalized position, or simply the position of Pj. This position can
be expressed relative to another position{Pi, Ei} by the (3 × 1) position vector ri,ij and the
(3× 3) rotation matrix Ri
j. The position vector r i,i
j defines the position of Pj(lower index j )
relative to Pi(second upper index i) and its components are expressed in the coordinate
sys-tem{Pi, Ei} (first upper index i). The rotation matrix Ri
jdefines the orientation of Ej(lower
index j ) relative to Ei(upper index i) expressed in{Pi, Ei}. Figure2shows a graphical
rep-resentation of the position of a material point using the position vector and rotation matrix. The two notations{Pj, Ej} and {ri,ij ,Rij} will both be used throughout this work, depending
on the context, to identify a position. Fig. 2 The position of Pjwith
respect to Piusing a position
vector and rotation matrix. Figure was made using InkScape
Ri
j also defines a coordinate transformation that can be used to transform a vector that
is expressed in coordinate system j to a vector that is expressed in coordinate system i. For example, the components of the position vector of Pj with respect to Piin frames{Pi, Ei}
and{Pj, Ej} are related as:
ri,ij = Ri jr
j,i
j , (2.1)
where it is seen that in this notation the first upper index has changed. The rotation matrix is an orthogonal matrix of the proper kind, which means that its determinant equals+1 and its transpose equals its inverse, which also represents the inverse coordinate transformation such that Ri j −1 = Rj i, RijR j i = I, (2.2)
with I the (3× 3) unity matrix. Differentiating (2.2), it can be shown that the time derivative of the rotation matrix equals a skew symmetric matrix times the rotation matrix itself. This can be written as ˙Ri j= ω i,i j R i j, (2.3)
in which the tilde operator is introduced such that when applied on a (3× 1) vector a, it yields a skew symmetric (3× 3) matrix ˜a:
˜a = ⎡ ⎣ a03 −a03 −aa21 −a2 a1 0 ⎤ ⎦ . (2.4)
In (2.3), the tilde operator is applied on the vector ωi,ij , which is the angular velocity vector of frame{Pj, Ej} with respect to {Pi, Ei} with its components expressed in {Pi, Ei}. The
linear velocity of{Pj, Ej} with respect to {Pi, Ei} expressed in {Pi, Ei} is simply denoted
by˙ri,ij . The linear and angular velocities can be combined in the (6× 1) velocity vector vi,ij as follows: vi,ij ≡ ˙ri,i j ωi,ij . (2.5)
To make a clear distinction between the inertial reference frame and any other coordinate system,{PO, EO} is used to denote the inertial frame. In this case {rO,Oj ,RO
j} defines the
absolute position of material point Pjand vO,Oj defines its absolute velocity.
3 Kinematics of a material point using the floating frame formulation
Figure3shows a flexible body to which a floating frame is attached in material point Pj.The absolute position of the floating frame with respect to the inertial frame located in PO
is denoted by{rO,Oj ,RO
j}. The position of another material point Pkrelative to the floating
frame is denoted by{rj,jk ,Rjk}. For the rotation from Pkto POholds:
RO k = R O jR j k. (3.1)
The position of Pkin the inertial frame can be expressed in terms of its relative position and
the absolute position of the floating frame as
Fig. 3 The position of material point Pkrelative to POusing
floating frame Pj. Figure was
made using InkScape
The absolute linear velocity of Pkin the inertial frame is found by taking the time derivative
of (3.2): ˙rO,O k = ˙r O,O j + ω O,O j R O jr j,j k + R O j ˙r j,j k . (3.3)
For the absolute angular velocity of Pkwe get
ωO,Ok = ωO,Oj + ROj ωj,jk . (3.4) Equations (3.3) and (3.4) can be rewritten and combined to
˙rO,O k ωO,Ok = RO j 0 0 RO j I −˜rj,jk 0 I RjO 0 0 RjO ˙rO,O j ωO,Oj + RO j 0 0 RO j ˙rj,j k ωj,jk . (3.5) For this, several relevant properties of skew symmetric matrices are used. A convenient overview of these properties is given in [11]. Equation (3.5) can be written in compact form using the definition of the velocity vector (2.5):
vO,Ok = RO j −˜rj,j k RjO vO,Oj + RO j vj,jk . (3.6)
Here the notations[RO
j] and [−˜r j,j
k ] are introduced to simplify and shorten the notation.
They represent the (6× 6) compound rotation matrix and (6 × 6) compound position ma-trices in Eq. (3.5), respectively. The velocity vj,jk is the local velocity of Pkcaused by the
elastic deformation of the body expressed in the floating frame. Equation (3.6) can be refor-mulated as
vj,jk = RjO vO,Ok − −˜rj,jk RjO vO,Oj . (3.7) Equation (3.7) shows that the relative velocity of the elastic deformation in the floating frame is defined by the difference of the absolute velocities in Pkand Pj.
4 Relation between the local elastic velocities and absolute velocities
It is assumed that the kinematics of a flexible body in a floating frame can be described with linear finite element models. The material behavior satisfies Hooke’s law and the strains are approximated with the linear Cauchy strain tensor [15,16]. Based on this model, the dynamic behavior of a flexible body is fully described in the floating frame{Pj, Ej} by abe applied by assuming that the deformation of a body is a linear combination of selected modes. Here the static Craig–Bampton modes (also known as interface modes or boundary modes) will be used as a reduction basis. The generalized coordinates qj,jk corresponding with these modes are the small elastic displacements uj,jk and rotations θj,jk of the interface points Pkon the body:
qj,jk = θj,jk uj,jk . (4.1)
The position vector rj,jk is the relative location of interface point Pkin the body reference
frame{Pj, Ej} and can be expressed as
rj,jk = xj,jk + uj,jk xj,jk , (4.2) where the position vector xj,jk is the location of Pkon the undeformed body, which is
con-stant. Using (2.2) in rewriting (4.2), the elastic displacement uj,jk can be expressed in terms of the absolute position of Pkand the absolute position of the floating frame Pj:
uj,jk = rj,jk − xj,jk = RjOrO,Ok − rO,Oj − xj,jk . (4.3) For small deformations, the rotation matrix Rjk can be related to the nodal rotations θj,jk directly, since within the framework of linear theory all parameterizations of Rjk yield the same results. To this end, consider the Taylor expansion of Rjk at t= t0+ t around the
undeformed configuration and use the exponential map [17] to approximate the rotation matrix as:
Rjkθj,jk + tωj,jk = expθkj,j+ tωj,jk ≈ I +θj,jk + tωj,jk . (4.4) Since the vector ωj,jk defines the local angular velocity at t= t0, Eq. (4.4) holds for small t .
With this relation, it follows that the nodal rotations θj,jk at t= t0can be derived from the
rotation matrix as θj,j k ≈ 1 2 Rjk− Rj Tk . (4.5) Moreover, the time derivative of θj,jk at t= t0is approximately the angular velocity
vec-tor ωj,jk . This gives a relation between the relative velocity and the generalized coordinates corresponding to the Craig–Bampton modes:
˙qj,j k = ˙θj,j k ˙uj,j k ≈ ωj,jk ˙uj,j k = vj,j k . (4.6)
This relation, in combination with (3.7) shows that it is possible to relate the relative velocity of a material point˙qj,jk to the absolute velocities vO,Ok and vO,Oj :
˙qj,j
k =
RjO vO,Ok − −˜rkj,j RjO vO,Oj . (4.7) Using a floating frame formulation, the dynamic equations contain the absolute coordinates of the floating frame, as well as the relative coordinates of the interface points. With result
(4.7), it is possible to eliminate the local elastic velocities from the dynamic equations and replace them by the absolute velocities of the floating frame and interface points [18]. At this point, the bodies can be coupled directly and Lagrange multipliers are no longer required to enforce constraints.
Because in this work the floating frame is not located at an interface point, the Craig– Bampton modes can still describe six rigid body motions. As the rigid body motion is al-ready described by the floating frame, including both makes the system singular. Six addi-tional constraints should be imposed on the Craig–Bampton modes in order to remove this singularity. These additional constraints enable the possibility to express the motion of the floating frame in terms of the motion of the interface points Pkonly. This yields a multibody
formulation solely in terms of the absolute coordinates of the interface points.
5 Elimination of the floating frame from the kinematic description
Figure4shows a flexible body with floating frame Pj and two interface points Pkand Pl.
The relative positions of the interface points with respect to the floating frame depend on the Craig–Bampton degrees of freedom qj,jk and qj,jl . In the floating frame formulation, the equations of motion are expressed in terms of the absolute floating frame coordinates {rO,O
j ,ROj } and the relative coordinates q j,j k and q
j,j
l (dashed arrows). The idea is to
refor-mulate the equations of motion by replacing these coordinates by the absolute motion of the interface points, defined by{rO,Ok ,RO
k} and {r O,O
l ,R
O
l } (solid arrows).
Suppose that the body has N interface points and that Pj is not an interface point. The
deformation of Pj due to the Craig–Bampton modes related to interface point Pkare
col-lected in the (6× 6) modes matrix Φjk. The elastic deformation qj,jj in point Pj is then a
superposition of Φjkand the generalized coordinates qj,jk in the interface points:
qj,jj =
N
k=1
Φjkqj,jk . (5.1)
In order to remove the 6 rigid body modes from the set of Craig–Bampton modes, the re-quirement is made that the elastic body has no elastic deformation in Pj. Unfortunately,
since the relation between the relative position of the interface points and the absolute posi-tions of the interface points and the floating frame is nonlinear, these additional constraints cannot be solved explicitly on the position level. However, this can be done on the velocity Fig. 4 Local position of the
interface coordinates in the floating frame (dashed arrows) and absolute position of the interface coordinates in the inertial frame (solid arrows). Figure was made using InkScape
level. To this end, Eq. (4.7) is substituted in the time derivative of (5.1). It follows from the requirement of zero elastic deformation that this must also be zero:
˙qj,j j = N k=1 Φjk RjO vO,Ok − −˜rj,jk RjO vO,Oj = 0. (5.2) Rewriting yields a relation between the velocity of the floating frame and the absolute ve-locities of the interface points:
Qj RjO vO,Oj = N k=1 Φjk RjO vO,Ok , Qj ≡ N k=1 Φjk −˜rj,jk . (5.3) The (6× 6) matrix [Qj] consists of 6 elastic displacement vectors defined by the sum of the
Craig–Bampton modes multiplied with the local positions[−˜rj,jk ] of the interface points. The matrices[−˜rj,jk ] define the 6 rigid body modes of the deformed body in the interface points calculated relative to Pj. When the body is undeformed, the position vectors rj,jk equal
the position vectors xj,jk and the resulting displacements in[Qj] become unit displacements
and so it equals the unity matrix I. Since rj,jk equals xj,jk on leading order,[Qj] can be
assumed regular in all cases and can be inverted. Hence, (5.3) can be solved for the absolute velocity of the floating frame:
vO,Oj = ROj Qj −1
N
k=1
Φjk RjO vO,Ok . (5.4) This can be written in compact matrix–vector notation as follows:
vO,Oj = RO j [Zj] RjO vO,O. (5.5)
Here the (6N× 1) vector vO,Ocontains the absolute velocities of all interface points and the
(6N× 6N) rotation matrix [RjO] is assembled from the (6 × 6) rotation matrix [RjO]:
vO,O≡ ⎡ ⎢ ⎢ ⎣ vO,O1 .. . vO,ON ⎤ ⎥ ⎥ ⎦ , RjO ≡ ⎡ ⎢ ⎢ ⎣ [Rj O] . .. [Rj O] ⎤ ⎥ ⎥ ⎦ . (5.6)
Moreover, the (6× 6N) matrix [Zj] is defined using a (6 × 6N) matrix [Φj
CB] as follows: Zj ≡ Qj −1 Φj CB , ΦjCB ≡ Φj1 . . . ΦjN . (5.7) Upon substitution of (5.5) into (4.7), it is possible to express the local velocities of the deformed body in terms of the absolute velocities of the interface points, which can be written compactly as ˙qj,j=[Φj rig] I −[Zj] I RjO vO,O. (5.8)
Here the (6N× 1) vector ˙qj,j contains the relative velocities of all interface points due to
elastic deformations of the body expressed in the floating frame and the (6N× 6) matrix [Φj
rig], which is the matrix of rigid body modes, is defined as
˙qj,j≡ ⎡ ⎢ ⎢ ⎣ ˙qj,j 1 .. . ˙qj,j N ⎤ ⎥ ⎥ ⎦ , Φjrig ≡ ⎡ ⎢ ⎢ ⎣ [−˜rj,j 1 ] .. . [−˜rj,j N ] ⎤ ⎥ ⎥ ⎦ . (5.9)
In short, (5.8) can be written as ˙qj,j= Tj Rj O vO,O, Tj ≡[Φj rig] I −[Zj] I = I − Φjrig Zj . (5.10)
The combination of Eqs. (5.5) and (5.10) makes it possible to rewrite the dynamic equa-tions of motion in floating frame formulation to a formulation in terms of the inertial frame formulation. This is demonstrated in the next section.
6 The equations of motion in absolute interface coordinates
The derivation of the equations of motion in terms of the absolute interface coordinates are derived based on the equations of motion in the floating frame formulation. For each flexible body in a multibody system, the standard equations of motion can be written as
M¨q + Qv+ Kq = Qe+ Qc. (6.1)
Here M and K are the mass and stiffness matrices, respectively. The vector of generalized coordinates q is the set of the floating frame coordinates and the local elastic coordinates corresponding to the Craig–Bampton modes. Qvis the vector of quadratic velocity inertia
forces, Qeis the vector of externally applied forces and Qcis the vector of constraint forces.
A detailed derivation of Eq. (6.1), based on the principle of virtual work, can be found in standard textbooks on multibody dynamics, e.g., [12]. Equation (6.1) can be partitioned as follows: Mrr Mrf Mf r Mff aO,Oj ¨qj,j + Qv,r Qv,f + 0 0 0 Kff 0 qj,j = Qe,r Qe,f + Qc,r Qc,f . (6.2)
Here the subscripts r and f refer to the rigid and flexible partitions of the matrices and vectors, respectively. aO,Oj is the absolute acceleration of the floating frame. The matrices Mff and Kff are directly obtained from the linear finite element model of the body on
which the Craig–Bampton reduction is applied.
In order to rewrite Eq. (6.2) in terms of the absolute interface coordinates, a coordinate transformation is used. On the velocity level, this coordinate transformation is obtained from combining Eqs. (5.5) and (5.10):
vO,Oj ˙qj,j = ⎡ ⎣[ROj][Zj][R j O] [Tj][Rj O] ⎤ ⎦ vO,O= AvO,O, (6.3)
where A shortly denotes the coordinate transformation matrix. Differentiating (6.3) with respect to time yields the coordinate transformation on the acceleration level. By means of the product rule, this transformation consists of a transformation of the accelerations themselves and a squared velocity term containing the time derivative of the transformation matrix:
aO,Oj
¨qj,j
= AaO,O+ ˙AvO,O. (6.4)
Substituting this coordinate transformation into (6.2) yields the equation of motion in terms of the absolute interface coordinates. By direct computation of the matrix multiplications concerned with this coordinate transformation, it is found that the resulting equation of motion can be written in the following form:
ROj [Mff]
RjO aO,O+ QOv + ROj Tj T[Kff]qj,j= QOj. (6.5)
Here the main inertia term can be recognized as the body’s local Craig–Bampton mass ma-trix rotated to the global frame multiplied with the absolute accelerations of the interface points. Other inertia forces that result from the coordinate transformation (6.4) are quadratic in the velocity and are therefore combined with the quadratic velocity inertia terms in the original equation of motion and represented by QO
v. This term is not taken into account
in standard corotational formulations. The procedure as explained here could be used to in-clude these additional inertia effects in the corotational formulation as well. However, taking this term into account is only significant in case of systems operating at high velocities. In the elastic term, the Craig–Bampton stiffness matrix is multiplied with the local elastic co-ordinates. The pre-multiplication with[Tj]T can be interpreted as an operation that extracts
the elastic deformations from the absolute interface coordinates by eliminating the rigid body component, see also the definition of[Tj] in Eq. (5.10). The resulting elastic forces
are transformed to the global frame by the rotation matrix. QO
j contain the forces applied on
the interface points expressed in the global frame. The constraint forces are also included simply in QO
j , because the constraints are enforced at the location of the interface points.
As explained in Sect.5, it is not possible to construct an explicit coordinate transforma-tion between the floating frame coordinates, local interface coordinates and absolute inter-face coordinates on the position level. Hence, due to the elastic term the equation of motion (6.5) is not entirely formulated in terms of the absolute interface coordinates. However, when numerically integrating in time, the equation of motion does not need to be solved for the large absolute position of the interface points. Instead, it is only solved for the small incre-ment in the interface coordinates that occurs during the time increincre-ment. The time-discretized equations are linear in this position increment and tangent to the current configuration space. For that reason, the same coordinate transformation matrix as in (6.3) can be applied on this position increment. In this way, the problem that needs to be solved at every time step is formulated completely in terms of the absolute interface coordinates.
The location of the floating frame at the next time step can be determined from integrating (5.5). Theoretically, however, the numerical error may cause the floating frame to drift. For that reason, the location of the floating frame is determined from the absolute interface coordinates on the position level using a Newton–Raphson procedure. As an initial estimate, the floating frame position of the current time step is used. In practice, only few Newton– Raphson iterations are required as the time steps are sufficiently small. With the absolute position of the interface coordinates and the floating frame determined, the local elastic deformation can be determined directly.
The presented method can also be used to solve large deformation, static problems. By simplifying (6.5) for this case, the equilibrium equations for one body are obtained:
ROj [Tj]T[Kj]qj,j= QOj. (6.6)
These equations are nonlinear in the elastic deformations qj,jas both[RO
j ] and [Tj] depend
on them. However, these equations can be solved incrementally using a standard updated Lagrangian approach. In this way, equations are obtained in terms of the position increment, which can be rewritten in terms of the absolute interface coordinates, similarly as discussed for the dynamic case.
7 Validation
In order to validate the method presented in this work, its solution to four problems is com-pared with that of other standard software. These validation problems consist of a static cantilever beam subjected to a large vertical tip force, a 2D and 3D slider–crank mechanism with a flexible connector and a 3D spinning beam on a spherical joint. The solution of the static problem is compared with Spacar and Ansys. The solutions of the dynamic problems are compared with Spacar and MSC/Adams. Spacar is a finite element based multibody soft-ware that uses the corotational formulation [19]. Ansys uses an inertial frame formulation for their nonlinear finite element analyses. MSC/Adams uses a floating frame of reference formulation.
For the static problem, a cantilever beam with circular cross section was modeled with 10 bodies. The total length of the beam is 1 m. The outer radius of the cross-section is 0.01 m has a wall thickness of 0.001 m. The Young’s modulus is 70.0E9 Pa. The beam was incrementally loaded starting at 100 N and increasing to 10000 N. Results have been obtained with the new method, Spacar and Ansys. The computed deformed beam shapes are shown in Fig.5for applied loads of 100, 500, 2000, and 10000 N. The figures show good agreement between the new code and both Spacar and Ansys.
The dynamic 2D slider–crank problem is adopted from [20] and shown in Fig.6. The rigid crank with length of 0.15 m is rotating with a constant angular velocity of 150 rad/s. The flexible connector with length of 0.3 m has a uniform circular cross-section with a diameter of 0.006 m and is made of steel. In the simulation a Young’s modulus of 0.2E12 Pa and a mass density of 7.87E3 kg/m3is used. The end of the connector is connected to
a slider with a mass half the mass of the connector. The slider is able to translate without friction on its base.
The angular velocity of the crank introduces an initial linear velocity and an angular velocity of the connector, assuming no deformation. The connector is modeled with two bodies. Hence, this model has three nodes: the two interface nodes and a half-way node at the location of the floating frame.
As output the displacement of the midpoint perpendicular to the undeformed beam was determined and plotted against the crank angle. The results are shown in Fig.7. This figure also shows the results obtained with Spacar in which longitudinal deformations due to nor-mal forces are suppressed. Moreover, results obtained with MSC/Adams are included. As Fig.7shows, the new method agrees very well with the results obtained with Spacar. The results obtained with MSC/Adams show small differences.
The dynamic 3D slider–crank mechanism is adopted from [21] and shown in Fig.8. The physical properties of the mechanism are the same as in the 2D case described above. The
Fig. 5 Deflection of a cantilever beam subjected to a tip force compared with Spacar and Ansys. Figure was made using Matlab for plotting and Adobe Illustrator for labels
Fig. 6 2D Slider–crank mechanism with flexible connector. Figure was made using InkScape
Fig. 7 Mid-point deformation of the connector compared with Spacar and MSC/Adams. Figure was made using Matlab for plotting and Adobe Illustrator for labels
horizontal position of the rotation axis d is 0.1 m. In the initial configuration, the crank is oriented vertically upward.
Fig. 8 3D Slider–crank mechanism with flexible connector. Figure was made using InkScape
Fig. 9 Mid-point deformation of the connector compared with Spacar and MSC/Adams. Figure was made using Matlab for plotting and Adobe Illustrator for labels
Fig. 10 Flexible beam on a spherical joint. Figure was made using InkScape
As output the displacement of the midpoint in its local y-direction was determined and plotted against the crank angle. The results are shown in Fig.9. It can be seen that also in this case the results obtained with the new method are very close to the results obtained with Spacar. The results obtained with MSC/Adams again show a small difference in comparison with the other two methods.
The 3D spinning beam on a spherical joint is adopted from [14] and shown in Fig.10. The physical properties, prescribed loads and simulation settings are the same as described in this reference: The beam has length 141.42 mm, cross-section 9.0 mm2and moment of
inertia 6.75 mm4. The material has a mass density 7.8E-3 kg/mm−3 and Young’s Modulus
2.1E6 N/mm2. The beam is modeled with two bodies. A torque of 200 N mm is applied
about the vertical axis during the first 10.2 seconds. After 15 seconds, an impulsive vertical tip force of 100 N is applied.
Fig. 11 Angular velocity about the vertical axis at the base of the beam compared with Spacar and MSC/Adams. Figure was made using Matlab for plotting and Adobe Illustrator for labels
As output the absolute angular velocity about the vertical axis at the base of the beam was determined and plotted as a function of time. The results are shown in Fig.11. It was observed in [14] that different methods show different results only after the impulsive verti-cal force is applied. In Fig.11it can be seen that all methods used here produce very similar results even after this moment. Comparing all results with those published in [14], it can be seen that the new method presented in this work is most comparable to results of the nonlinear finite element formulation used in [14].
8 Conclusion
Describing the kinematics of a flexible multibody system comes down to the kinematic for-mulation of the motion of the interface points. In the inertial frame and corotational frame formulations, the absolute interface coordinates are part of the degrees of freedom, allowing for a direct application of the constraints. This is in contrast to the floating frame formula-tion, which requires the use of Lagrange multipliers. In this work, it has been demonstrated that the absolute floating frame coordinates and the local elastic coordinates that appear in the equation of motion of a floating frame formulation can be replaced by the absolute in-terface coordinates. Consequently, the new method does not require Lagrange multipliers to enforce the kinematic constraints. In this way a floating frame formulation in terms of a minimal set of coordinates is successfully obtained. Validation of the method with static and dynamic problem found in literature has shown to yield reliable results in all cases.
The use of Craig–Bampton modes as local shape functions is crucial in the presented procedure. The rigid body motions contained in these Craig–Bampton modes are employed to eliminate the floating frame coordinates from the system. The problem of establishing the kinematic coordinate transformation presented in this work for the floating frame formula-tion is in fact comparable to the problem of expressing the corotaformula-tional frame coordinates in terms of the element’s absolute nodal coordinates for a corotational frame formulation. Only here this problem is encountered at the level of an entire body instead of on the level of a single element. An advantage of the new method is that it can be applied to systems that consist of arbitrarily shaped three-dimensional bodies that have an arbitrary number of in-terface points. Also in this sense, it can be seen as a generalization of the corotational frame formulation, which is only developed for a limited number of parameterized finite elements, such as beams, plates and shells.
As such, this paper not only provides valuable insight in the relations between differ-ent multibody formulations, but it also offers the possibility to reduce geometric nonlinear systems by applying important modal order reduction techniques in a body’s local frame. This is found a convenient way of creating so-called superelements in a flexible multibody formulation. The use of the term superelement emphasizes the striking similarity between the floating frame and corotational frame formulations at this point. For each flexible body, the tangent mass and stiffness matrices, reduced to the interface points can be obtained from detailed linear finite element models. These system matrices can directly be applied in the flexible multibody analysis of the entire system.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Inter-national License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
References
1. Wasfy, T.M., Noor, A.K.: Computational strategies for flexible multibody systems. Appl. Mech. Rev. 56(6), 553–613 (2003)
2. de Borst, R., et al.: Non-linear Finite Element Analysis of Solids and Structures, 2nd edn. Wiley, New York (2012)
3. Géradin, M., Cardona, A.: Flexible Multibody Dynamics, a Finite Element Approach. Wiley, New York (2001)
4. Jelenic, G., Crisfield, M.A.: Interpolation of rotational variables in nonlinear dynamics of 3D beams. Int. J. Numer. Methods Eng. 43, 1193–1222 (1998)
5. Romero, I., Armero, F.: An objective finite element approximation of the kinematics of geometrically exact rods and its use in the formulation of an energy-momentum conserving scheme in dynamics. Int. J. Numer. Methods Eng. 54, 1683–1716 (2002)
6. Rankin, C.C., Brogan, F.A.: An element independent corotational procedure for the treatment of large rotations. J. Press. Vessel Technol. 108, 165–174 (1986)
7. Crisfield, M.A.: A consistent co-rotational formulation for non-linear, three-dimensional, beam elements. Comput. Methods Appl. Mech. Eng. 81, 131–150 (1990)
8. Moita, G.F., Crisfield, M.A.: A finite element formulation for 3-D continua using the co-rotational tech-nique. Int. J. Numer. Methods Eng. 39, 3775–3792 (1996)
9. Craig, R.R., Bampton, M.C.C.: Coupling of substructures for dynamic analysis. AIAA J. 6(7), 1313– 1319 (1968)
10. Craig, R.R.: Coupling of substructures for dynamic analyses: an overview. In: Structures, Structural Dynamics and Material Conference, 41st AIAA/ASME/ASCE/AHS/ASC, Atlanta (2000). AIAA-2000-1573
11. Haug, E.J.: Computer Aided Kinematics and Dynamics of Mechanical Systems: Basic Methods. Prentice Hall, New York (1989)
12. Shabana, A.A.: Dynamics of Multibody Systems. Cambridge University Press, Cambridge (1998) 13. Cardona, A., Géradin, M.: Modeling of superelements in mechanism analysis. Int. J. Numer. Methods
Eng. 32, 1565–1593 (1991)
14. Cardona, A.: Superelements modelling in flexible multibody dynamics. Multibody Syst. Dyn. 4, 245– 266 (2000)
15. Sokolnikoff, I.S.: Mathematical Theory of Elasticity. McGraw-Hill, New York (1956) 16. Timoshenko, S., Goodier, J.N.: Theory of Elasticity. McGraw-Hill, New York (1934)
17. Cardona, A.: An integrated approach to mechanism analysis. Ph.D. Thesis, Université de Liege (1989) 18. Ellenbroek, M.H.M.: On the fast simulation of the multibody dynamics of flexible space structures. Ph.D.
Thesis, University of Twente (1994)
19. Jonker, J.B., Meijaard, J.P.: A geometrically non-linear formulation of a three-dimensional beam element for solving large deflection multibody system problems. Int. J. Non-Linear Mech. 53, 63–74 (2013) 20. Jonker, J.B.: Dynamics of Spatial Mechanisms with Flexible Links. Ph.D. Thesis, Technical University
of Delft (1988)
21. IFToMM Technical Committee for Multibody Dynamics: Library of Computational Benchmark Prob-lems.www.iftomm-multibody.org/benchmark, consulted on 14-09-2017
