SSRN 4830498
SSRN 4830498
SSRN 4830498
d
Rotordynamics Local Frame Formulations for 3D Con-
tinuum Finite Elements
we
Andreas Zwölfer1 , Jacob Østerby Holst Rasmussen2 , Ilmar Santos2
1
Technical University of Munich, TUM School of Engineering and Design, Department of Mechanical Engineering,
Chair of Applied Mechanics, Boltzmannstr. 15, 85748 Garching b. München, Germany; Corresponding author:
[email protected]
ie
2
Technical University of Denmark, Department of Civil and Mechanical Engineering Solid Mechanics, Koppels
Allé, 404, 2800 Kgs. Lyngby, Denmark
ev
Abstract
This paper presents the derivation of the nodal-based floating frame of reference formulation equations of motion
stemming from flexible multibody dynamics tailored for rotordynamics applications characterized by large rota-
r
tions around a single spatially fixed axis. In addition to this fundamental assumption, further constraints, such
as those concerning the rigid body translation of the floating frame, are directly integrated into the equations of
motion. By gradually incorporating additional assumptions – reducing complexity and increasing efficiency – the
er
well-established rotating frame rotordynamics formulation emerges from the floating frame of reference formula-
tion. The presented formulations are well-suited either for transient analysis, without necessitating 3D rotations
and differential algebraic equations, or for linearized dynamic analyses, including modal, harmonic, and stability
pe
assessments. They are compatible with any displacement-based finite element and are presented both with and
without component mode synthesis reduction. Moreover, all necessary terms within the equations of motion, such
as spin softening and Coriolis effects, can be directly computed from the assembled finite element mass and
stiffness matrices, as well as nodal reference coordinates.
Keywords: Rotordynamics, Floating Frame of Reference Formulation, Continuum/Solid Finite Elements, Com-
ponent Mode Synthesis, Modal Reduction
ot
1 Introduction
tn
Rotordynamics – the study of the dynamics of rotating systems – plays a crucial role in various engineering
applications, including aerospace, automotive, and energy industries. The accurate prediction and analysis of the
behavior of rotating machinery are essential for ensuring their reliability, efficiency, and safety.
In pursuit of enhanced reliability and efficiency in product optimization, industries are increasingly turning to vir-
tual prototyping over traditional physical testing methods. Advanced engineering simulation tools have, therefore,
rin
beam theory. This approach typically assumes the structure to be axisymmetric, and special treatment is needed
to model structures of asymmetric geometry. In cases where shear effects are significant or geometric features do
not allow for utilizing shaft elements, 3D solid elements, which are surprisingly rarely used in the rotordynamics
literature, are needed [8, 9, 10]. Clearly, the utilization of three-dimensional solid finite element methods would
Pr
offer enhanced accuracy and simplified integration of detailed computer-aided design (CAD) models. Certain 3D
solid rotordynamics modeling options are also available in commercial software such as Abaqus and Ansys for
more than a decade [11], and computational efficiency concerns of these more detailed models may be alleviated
with model reduction techniques [12, 13, 14, 15].
Typically, the objective of rotordynamics formulations is to extend the standard linear elastodynamics finite
element (FE) mass and stiffness matrices with the dynamic effects of interest, e.g. gyroscopic, Coriolis, centrifugal
1
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
softening/stiffening, and perform linearized dynamic analyses, such as harmonic, modal, and stability. These
linear equations of motion (EOMs) employed for typical rotordynamics analyses may be written in the stationary or
rotating frame [4], where the former is in any case limited to axisymmetric rotors, while the latter is conventionally
d
limited to a single rotating structure in an analysis [11]. The rotordynamics EOMs are usually derived from a choice
of FE type and, therefore, restricted to modelling the rotor using the same element type across the whole rotor
we
domain.
The current contribution attempts to solve the aforementioned issues from a flexible multibody (MB) dynamics
perspective. A flexible MB system may be defined as a group of interconnected rigid and deformable components,
which are subjected to large rigid body (RB) translations and rotations, i.e., “bulk motion”, during their intended
operation [16, 17]. Flexible MB dynamics usually refers to the computational strategies employed to determine
(long) time histories of motion, i.e., transient analyses / time integration, but the formulations per se may be used in
ie
rotordynamics applications also characterized by a large RB rotation with superimposed vibrations/deformations.
The EOMs of flexible MB bodies are usually derived via a straightforward extension of RB dynamics, i.e., the
kinematic description involves a RB translation and rotation part, as in RB dynamics, and a deformation term
ev
is simply added. Conventionally the degrees of freedom (DOFs) of a flexible body are then the translation and
rotation of the underlying RB motion / moving frame as well as local deformation measures. The formulation
obtained with these local flexible DOFs is referred to as floating frame of reference formulation (FFRF). The FFRF
has become standard for modeling flexible MBs systems and is without doubt the most widely-used method and
implemented in most commercial software packages, such as RecurDyn, Adams, Simcenter 3D, MotionSolve,
r
SIMPACK. The FFRF can be used with any element type and to set-up the EOMs only the linear elastodynamic
FE mass and stiffness matrix as well as nodal reference coordinates are needed with the recently introduced
er
nodal-based approach [18, 19, 20] implemented in the open-source research code Exudyn [21]. As an additional
advantage, transient phenomena and support motion excitation [22] may be studied with the FFRF.
The literature on rotor dynamics analyzed with the FFRF is sparse [23, 24, 25, 26], which may be due to
the fact that differential algebraic equation (DAE) solvers and 3D rotation treatment are conventionally necessary.
pe
However, this may be circumvented by embedding the main rotordynamics assumption, i.e., large rotation around
one axis fixed in space only, directly within the FFRF EOMs – this is the main novelty of the current contribution
and significantly reduces the complexity and EOMs terms needed in the FFRF leading to a significant increase in
efficiency.
The remainder of this paper is structured as follows: The paper starts with a summary derivation of the nodal-
based FFRF [18, 19] in Sect. 2.1. After that a formulation with large rotation around one spatially fixed axis but
ot
rigid body translation DOFs is derived in Sect. 2.2. Then, it is assumed that rigid body translation only takes place
in a plane orthogonal to the rotation axis – the corresponding formulation is presented in Sect. 2.3. Subsequently,
a formulation without rigid body rotation is presented in Sect. 2.4. And finally, a formulation with prescribed rigid
tn
body motion is derived in Sect. 2.5 – this formulation is equivalent to conventional rotating frame rotordynamics
formulations if it is assumed that the rigid body translation and angular acceleration are zero. Conclusive remarks
are provided in Sect. 3. Appx. A shows how the gravitational nodal forces may be easily calculated with the
assembled FE mass matrix, and Appx. B depicts some of the to the proposed formulations’ specific matrices and a
Campbell plot for a simple example to foster understanding and help those interested to check their implementation.
rin
Throughout the contribution, emphasis is placed on considerations valid for generic geometries modeled with 3D
slid finite elements with displacement degrees of freedom.
respect to (w.r.t.) the origin of the global inertial frame F and their orientations are related by the rotation matrix
A = A(θ ) ∈ R3×3 , where θ ∈ Rnr ×1 denotes a proper rotational parametrization with nr rotational generalized
coordinates.
FFRF splits the global nodal displacements c ∈ R3nn ×1 into translational, rotational, and flexible parts, i.e.,
2
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
[18, 27]
c = Φt τ + Abd − Ibd x + Abd c f , (1)
d
with
we
T
∈ R3nn ×3 ,
Φt = I ... I (2)
3nn ×3nn
Abd = diag(A, . . . , A) ∈ R , (3)
Ibd = diag(I, . . . , I) ∈ R3nn ×3nn , (4)
where I ∈ R3×3 denotes the identity matrix and x ∈ R3nn ×1 as well as c f ∈ R3nn ×1 denote the reference nodal
ie
coordinates and the flexible nodal displacements w.r.t. the floating frame1 , respectively. Differentiating (•)
˙ Eq. (1)
w.r.t. time t and exploiting some inherent properties of the involved matrices as well as some basic kinematic
relationships, see [18, 27], yields
ev
ċ = L(q )q̇ , (5)
r
c>
q = τ> θ > f
, (6)
and er
L = Φt −Abd e
r fG Abd , (7)
where G ∈ R3×nr is implicitly defined by the local angular velocity ω ∈ R3×1 and rotational parametrization as
pe
(i)
ω = G θ̇ , and e r f ∈ R3×3 of all FE nodes associated with the
r f comprises the nn skew-symmetric matrices2 e
nodal position vectors
(i) (i)
r f = x (i) + c f ∈ R3×1 ⇔ r f = x + c f ∈ R3nn ×1 (8)
(1)
rf
e
.
3nn ×3
rf =
e .
. ∈R . (9)
tn
(nn )
r
e
f
− + = Qrem , (10)
dt ∂ q̇ ∂q ∂q |{z}
| {z } | {z } remaining
inertia elastic forces
forces forces
1 >
T (ċ(q , q̇ )) =ċ M ċ, (11)
2
1
V (q ) = c > K cf, (12)
2 f
Pr
with the kinetic energy T , the strain energy V , and the remaining forces not stemming from T or V . Note that
M and K denote the constant FE mass and stiffness matrix from linear elastodynamics, respectively. It has been
1
Note that overlined quantities • are expressed in the floating frame in contrast to quantities expressed in the inertial frame •.
2
Note that the tilde operator (e
•) converts any R3×1 vector in its associated skew-symmetric R3×3 matrix.
3
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
shown by [18] that combining Eq. (1), Eq. (5) and Eqs. (10) to (12) and carrying out differentiation in an elegant
way yields
d
(13)
we
Φ>
t MΦt −AΦ>
t M r fG
e AΦ>
t M τ̈ 0 0
> > > >
r f M θ̈ = − 0 − 0
G e
rf M e
r fG −G e
sym. M c̈ f K cf D ċ f
ie
e ˙
−AΦ> t M ωbd ωbd r f + 2ωbd ċ f − r f G θ̇
f
e e e
> > res
+
˙ + >
G r M ωe ω e r + 2ω e ċ − e r G θ̇ G m res ,
e
f bd bd f bd f f (14)
e ˙
ev
−M ω e ω f
bd bd r f + 2ωbd ċ f − r f G θ̇
e e
where
ω
e = diag(ω, e ∈ R3nn ×3nn ,
e . . . , ω) (15)
bd
r
and with the resultant force vector fres ∈ R3×1 , resultant moment vector m res ∈ R3×1 , and applied nodal forces
f ∈ R3nn ×1 . In addition, linear damping forces with damping matrix D ∈ R3nn ×3nn where introduced in the flexible
part of the EOMs.
er
It is usually required to reduce the number of flexible DOFs within FFRF for efficiency reasons. The well-
established component mode synthesis (CMS) is a widely used approach to do so, see [13] for more information
pe
on the topic. Within CMS the flexible deformation is approximated by a linear combination of nm component
modes, such as vibration eigenmodes, static modes, etc., i.e.3
where Ψ ∈ R3nn ×nm contains column-wise the modes included in the reduction basis and ζ ∈ Rnm ×1 are the
ot
associated modal coordinates. It is important to note that if rigid body DOFs of the floating frame are present in
the formulation, the corresponding modes need to be removed from the reduction basis Ψ .
Substituting Eq. (16) into Eq. (14) and left-multiplying the resulting flexible part with Ψ > yields
tn
ΦTt MΦt −AΦTt M e
r fG AΦTt M Ψ τ̈ 0 0
T T T T
−G er M Ψ θ̈ = − 0 − 0
G e
r Me
f r G
f f
> > >
sym. Ψ MΨ ζ̈ Ψ K Ψζ Ψ D Ψ ζ̇
e ˙
−AΦTt M ω e ω
bd bd r f + 2ωbd ċ f − r f G θ̇
rin
fres
e e
T T
˙
>
+
G r
e
f M ω
e ω
bd
e r + 2ω
bd f
e ċ − e
bd f r f G θ̇ + G m
res
, (17)
>
> ˙ Ψ f
−Ψ M ω e ω
bd bd r f + 2ωbd ċ f − r f G θ̇
e e e
ep
c b
ff vf
| {z } | {z } | {z } | {z }
FFRF mass elastic & damping velocity-dep. resultant
matrix forces forces
pseudo forces
3
The equality sign is used in the following equations for the sake of simplicity even though the CMS introduces an approximation if
dim(ζ) < dim(c f ).
4
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
with the reduced FE stiffness and damping matrices and nodal forces known from linear elastodynamics, i.e.,
>
Kψψ = Ψ K Ψ ∈ Rnm ×nm , (19)
d
>
Dψψ = Ψ D Ψ ∈ Rnm ×nm , (20)
>
f ψ = Ψ f ∈ Rnm ×1 .
we
(21)
The velocity-dependent pseudo forces can be written in a more efficient form to avoid multiplications with “large”
matrices depending on the DOFs, i.e., [19]
˙ ˙
A ωΦe TM e
t r f ω + 2Φ T
t M c
e
f ω + Φ T
t M r
e
f G θ̇
ie
b =
Q T
T T ˙ T ˙ ,
−G ω r f M r f ω + 2 r f M c f ω + r f M r f G θ̇
v e e e e e e e (22)
T e T T ˙ T ˙
Iζ ⊗ ω Ψ Me r f ω + 2Ψ M e cfω + Ψ M e r f G θ̇
ev
since [18, 27, 19]
e > = −ω
ω e , (23)
bd bd
r
ω
e M =Mω
bd
e ,
bd (24)
Φ> e e >
t ωbd = ωΦt , er (25)
ω
e r = −e
bd f r f ω, (26)
>
rf M ω
e ω e e> e
bd bd r f = −ω r f M r f ω, (27)
e e
ω
e Ψ = −Ψ e I ⊗ ω,
pe
bd ζ (28)
ecf = Ψ ζ ⊗ I ,
e (29)
where
3n ×3n
Ψ
e= e
ψ1 . . . ψnm ∈ R n m ,
e (30)
ot
and ⊗ denotes Kronecker’s product. Note that Eq. (26) is valid for any R3nn ×1 block-nodal vector, such as c , x , c f ,
. . . ∈ R3nn ×1 , not only for r f .
Defining the following reduced and constant “inertia-like” matrices, [19]
tn
mI = Φ>
t MΦt ∈ R3×3 , (31)
e = 1 Φ> M xe ∈ R3×3 ,
χ (32)
u
m t
rin
>
Θ = xe M xe ∈ R3×3 ,
u (33)
>
Mψψ = Ψ M Ψ ∈ Rnm ×nm , (34)
>
e =Ψ MΨ
Mψψ e ∈ R3nm ×nm , (35)
ep
e>M Ψ
Mψeψe = Ψ e ∈ R3nm ×3nm , (36)
Mφψ = Φ>
t MΨ ∈ R3×nm , (37)
Mφ ψe = Φ>
t MΨ
e ∈ R3×3nm , (38)
Pr
>
M ex ψ = xe M Ψ ∈ R3×nm , (39)
>
M ex ψe = xe M Ψ
e ∈ R3×3nm , (40)
5
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
avoids large “online” matrix multiplications, since these so-called FFRF invariants may be precomputed once and
for all. Hence, [19]
c = mI,
d
M tt (41)
c = Mψψ ,
M (42)
ff
we
c = AMφψ ,
M (43)
tf
Mtr = −A mχ
c e + M e ζ ⊗ I G,
u φψ (44)
>
c = −G M ex ψ + ζ ⊗ I > M e ,
M rf ψψ (45)
>
c = G Θ + M e ζ ⊗ I + ζ ⊗ I > M >e + ζ ⊗ I > M e e ζ ⊗ I G,
M (46)
ie
rr u xψ
e xψ
e ψψ
and
b =Aω
Q e mχ e + M e ζ ⊗ I ω + 2AM e ζ̇ ⊗ I ω
ev
vt u φψ φψ
e + M e ζ ⊗ I G˙ θ̇ ,
+ A mχ u φψ (47)
| {z }
=0, if Euler parameters are used
r
Q vr u xψ
e xψ
e ψψ
> >
− 2G M ex ψe ζ̇ ⊗ I + ζ ⊗ I Mψeψe ζ̇ ⊗ I ω er (48)
>
> > > ˙
− G Θu + M ex ψe ζ ⊗ I + ζ ⊗ I M ex ψe + ζ ⊗ I Mψeψe ζ ⊗ I G θ̇ ,
| {z }
=0, if Euler parameters are used
pe
b = Iζ ⊗ ω > M >e + M e e ζ ⊗ I ω + 2M >e ζ̇ ⊗ I ω
Q vf xψ
e ψψ ψψ
> > ˙
+ M ex ψ + Mψψ
e ζ ⊗ I G θ̇ , (49)
| {z }
=0, if Euler parameters are used
˙ e ζ̇ ⊗ I ,
cf = Ψ
e (50)
˙
since Ψ, I = const., have been used. Note that G θ̇ = 0 if Euler parameters are used for the rotational
e
tn
−AΦTt M ω e ω
bd bd r f + 2ωbd ċ f
−AΦTt M e
r f AΦTt M τ̈
mI 0 fres
e e
T T T
r f M ω̇ = 0 + rf M ω e ω
bd bd r f + 2ωbd ċ f + m res , (51)
rf Me
e rf −e e e e
−K c f − D ċ f f
c̈ f −M ω e ω
bd bd r f + 2ωbd ċ f
sym. M e e
ep
since
˙
ω = G θ̈ + G θ̇ . (52)
Pr
We now assume large rotations around a fixed axis only, which is mathematically expressed as
ω = ϑ̇e (53)
e = const. (54)
6
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
and the scalar rotation angle/parameter ϑ. In this case Eq. (53),
d
and therefore
we
e = e, (56)
−AΦTt M e AΦTt M τ̈ −AΦTt M e e bd r f ϑ̇2 + 2e
e bde e bd ċ f ϑ̇
ie
mI r fe 0 fres
> eT T > eT e ∗
e r Me r e
f f r f M ϑ̈ =
−e > e 0 +
2e r f M e bd ċ f ϑ̇
+ mres ,
sym. M c̈ f −K c f − D ċ f −M e e r ϑ̇ + 2e ċ ϑ̇
e e 2 e f
bd bd f bd f
ev
(57)
or, to enable the pre-computation of the products of constant matrices when the non-reduced version is utilized,
AΦTt M e AΦTt M −AΦ T
M e e r ϑ̇ 2
+ 2 e ċ ϑ̇
r
τ̈
mI e bd r f 0 fres
e e e
t bd bd f bd f
>
e> ϑ̈ = + T e> ∗
r >e r> 0 + mres ,
e Me e r f e bd M −2r f e bd M e bd ċ f ϑ̇
e
f bd bd f
sym. M c̈ f −K c f − D ċ f
er −M e e e e r ϑ̇2 + 2e e ċ ϑ̇ f
bd bd f bd f
(58)
since there are no gyroscopic (pseudo-)forces for in plane rotations as is known from rigid body dynamics, i.e,
pe
Eq. (27) and Eq. (53),
T T
−e > e
r f Me e bd r f = e >e
e bde eerf M e
r fe (59)
= e · [e × (Θd e)] , (60)
T
Θd = Θd (c f ) = e r f.
rf M e (61)
Hence,
tn
r f e = −e
e e bd r f (63)
and
which represents the scalar moment around e , have been used to arrive at Eq. (58).
Introducing modal reduction as before, Eq. (16) and Eq. (29), and projecting the flexible part of Eq. (57)/Eq. (58)
with Ψ > yields the FFRF EOMs in terms of invariants for large rotations around one axis only, Eq. (53), as outlined
below. Note that Mctt = mI , M c> = AMφψ , and M
ctf = M cff = Mψψ remain unchanged as compared to the generic
Pr
ft
FFRF with modal reduction shown in Sect. 2.1.
The case-specific translation parts of the generalized mass matrix and velocity-dependent pseudo forces are
given by
c = A mχ ∗ + M ∗ ζ ,
M tr u φψ (65)
7
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
and
∗∗ ∗
¦ ©
b =−A
Q v mχ ∗∗
u + Mφψ ζ ϑ̇ 2
+ 2Mφψ ζ̇ ϑ̇ , (66)
t
d
where
χ ∗u = e
eχ u , (67)
we
χ ∗∗
u = e eχ u ,
ee (68)
∗
Mφψ =e
e Mφψ , (69)
∗∗
Mφψ =eeMφψ ,
ee (70)
ie
since, Eq. (24),
e bd = e
Me e bd M (71)
ev
and, Eq. (24),
Φt e
e bd = e
eΦt , (72)
due to the fact that M and Φt are composed out of blocks of the identity matrix and, Eq. (15),
r
e bd = diag e
e e, . . . , e
e . (73)
Furthermore, the case-specific rotation and flexible parts of the generalized mass matrix and velocity-dependent
er
pseudo forces are given by
c = − e T M ex ψ + ζ> Iζ ⊗ e T M e ,
M rf ψψ (74)
and
b = − 2 e T M e Iζ ⊗ e ζ̇ + ζ> Iζ ⊗ e > M e e Iζ ⊗ e ζ̇ ϑ̇,
Q vr xψ
e ψψ (76)
T T 2 T
b = Iζ ⊗ e
Q M ex ψe e + Mψeψe Iζ ⊗ e ζ ϑ̇ + 2Mψψ Iζ ⊗ e ζ̇ϑ̇, (77)
vf
e
ot
since
ζ ⊗ I e = Iζ ⊗ e ζ ⇒ ζ̇ ⊗ I e = Iζ ⊗ e ζ̇. (78)
tn
Introducing the constant and projected invariants, which again may be pre-computed, i.e.,
∗
Θu = e T Θu e ∈ R1×1 , (79)
∗ >
e = Iζ ⊗ e
Mψψ Mψψ
e ∈ Rnm ×nm , (80)
∗ >
Mψeψe Iζ ⊗ e ∈ Rnm ×nm ,
= Iζ ⊗ e
rin
Mψeψe (81)
∗
M ex ψ = e > M ex ψ ∈ R1×nm , (82)
∗
= e > M ex ψe Iζ ⊗ e ∈ R1×nm ,
M ex ψe (83)
finally yields
ep
c = −M ∗ − ζ> M ∗e ,
M (84)
rf xψ
e ψψ
∗ ∗ ∗
c =
M rr Θu + 2M ex ψe ζ + ζ> Mψeψe ζ, (85)
and
Pr
b = − 2 M ∗ e ζ̇ + ζ> M ∗e e ζ̇ ϑ̇,
Q vr xψ
e ψψ (86)
h > i >
b = M e + M e e ζ ϑ̇2 + 2M ∗e ζ̇ϑ̇,
Q
∗ ∗
(87)
vf xψ
e ψψ ψψ
∗
where Θ u is the scalar inertia around the rotation axis e .
8
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
2.3 Nodal-based FFRF with large rotation around a fixed axis and in plane rigid body translation
Again starting with the nodal-based FFRF in terms of the angular velocity, i.e., Eq. (51) and assuming ω = ϑ̇e , as
d
in Sect. 2.2, but now also
τ = Nτ∗ , (88)
we
where τ∗ ∈ R2 contains the preserved rigid body DOFs and
N = n1 n2 = null e >
(89)
and therefore contains two unit vectors, n1 and n2 , that are perpendicular to each other and also to the axis of
rotation e spanning the plane of the preserved rigid body DOFs. Note that null (•) denotes the null space of •.
ie
>
> e.g. 1 0 0
e.g.
For example, if e = 0 0 1 ⇒N = such that e and the columns of N form a right-handed
0 1 0
coordinate system.
ev
Substituting Eq. (53) and Eq. (88) into Eq. (51) and projecting the resulting translation part with N > yields, see
also Eq. (57),
r
f
sym. M c̈ f −K c f − D ċ f
−N > AΦTt M e e bd r f ϑ̇2 + 2e e bd ċ f ϑ̇
e bde
>
N fres
er + e> e
T
rf M e
e bdee bd r f ϑ̇2 + 2e
>
e bd ċ f ϑ̇ + e m res ,
−M e e bdee bd r f ϑ̇2 + 2e e bd ċ f ϑ̇ f
(90)
pe
In this case we define, see also Eq. (55),
A∗ = N > A = N > I cos(ϑ) + ẽ sin(ϑ) , (91)
since
e > N = 0>
⇔ N > e = 0.
ot
(92)
∗ cos ϑ − sin ϑ 0 >
For example, A = if e = 0 0 1 .
sin ϑ cos ϑ 0
Also, and in general,
tn
> 1 0 ∗
N N=I = . (93)
0 1
Hence, see also Eq. (58),
rin
∗ T e e
mI ∗ A∗ ΦTt M e A∗ ΦTt M −A Φt M e bd e bd r f ϑ̇2 + 2e e bd ċ f ϑ̇
∗
τ̈ 0∗ ∗
e bd r f fres
> e> e> ϑ̈ = + T e> ∗
r e Me
f e r
bd bd f r>
f e bd M
0 −2r f e bd M e
e
bd ċ f ϑ̇ + mres ,
sym. M c̈ f −K c f − D ċ f −M e e r ϑ̇ + 2e
e e 2
e ċ ϑ̇ f
bd bd f bd f
ep
(94)
c =M
where M c> , M
c ,Mc ,Qb ,Q
b in terms of invariants are equal as those in Sect. 2.2.
rf fr rr ff v v r f
Again, introducing modal reduction as before, Eq. (16), and projecting the flexible part of Eq. (94) with Ψ > yields
the specific FFRF EOMs terms in terms of invariants for large rotations around one axis, Eq. (53), and in-plane
Pr
9
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
and
∗∗ ∗
¦ ©
b = − A∗
Q v mχ ∗∗
u + Mφψ ζ ϑ̇ 2
+ 2Mφψ ζ̇ ϑ̇ . (98)
t
d
2.4 Nodal-based FFRF with large rotation around a fixed axis and no rigid body translational
we
acceleration
Again starting with the nodal-based FFRF in terms of the angular velocity, i.e., Eq. (51) and assuming ω = ϑ̇e , as
in Sect. 2.2, but now also
τ̈ = 0 (99)
ie
yields, see also Eq. (58),
> > T e>
∗
ϑ̈
r Tf e
e bd M e
e bd r f r Tf e
e bd M 0 −2r f e bd M e bd ċ f ϑ̇ mres
ev
e
= + + , (100)
sym. M c̈ f −K c f − D ċ f −M e bd e bd r f ϑ̇ + 2e bd ċ f ϑ̇
e e 2 e f
Again, introducing modal reduction as before, Eq. (16), and projecting the flexible part of Eq. (100) with Ψ >
c = M
yields the same terms for the flexible and rotation parts, M c> , Mc , M c, Q b , Qb , as in Sect. 2.2 and
r
rf fr rr ff v v r f
Sect. 2.3, but no rigid body translation contributions, i.e., no M
c ,M c =Mc> , M
c =Mc> , and Q
b .
tt tr rt tf ft v t
τ̈ = τ̈(t), (101)
and
ω = ω(t), (102)
ot
we have
e˙ + ω e˙
h i
>
M c̈ f + 2M ω
e + D ċ + M ω e ωe + K c = f − M ωe ωe +ω
bd x − Mφ A τ̈, (103)
tn
bd f bd bd bd f bd bd
or if ω = e ϑ̇ and e = const.
M c̈ f + 2M e
e bd ϑ̇ + D ċ f + M ee bd ϑ̈ + e e bd ϑ̇2 + K c f = f − M e
e bde e bd ϑ̈ + e e bd ϑ̇2 x − Mφ A> τ̈, (104)
e bde
rin
where
ep
Mφ = MΦt , (106)
and the multiplications of the constant matrices may be precomputed. These EOMs may be used for harmonic
analyses and their homogeneous versions for Campbell analyses [29].
Pr
Again, introducing modal reduction as before, Eq. (16), and projecting Eqs. (103) to (105) with Ψ > yields, see
also the definitions of the invariants in Sect. 2.1,
T T T
M ψψ ζ̈ + D ψψ − 2Mψψ
e Iζ ⊗ ω ζ̇ + K ψψ − Mψψ e Iζ ⊗ ω̇ − Iζ ⊗ ω Mψeψe Iζ ⊗ ω ζ = (107)
T T T >
= f ψ + M ex ψ ω̇ + Iζ ⊗ ω M ex ψe ω − M φψ A> τ̈,
10
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
due to Eq. (28) and Eq. (24), or if ω = e ϑ̇ and e = const.,
d
e e e
(108)
we
or if ω = e ϑ̇ and e = const. and ϑ̇ = const., (no rigid body angular acceleration)
∗> ∗>
∗ >
M ψψ ζ̈ + e ϑ̇
D ψψ − 2Mψψ ζ̇ + K ψψ − Mψeψe ϑ̇2 ζ = f ψ + M ex ψe ϑ̇2 − M φψ A> τ̈. (109)
Clearly, if τ̈ = 0 (no rigid body translation acceleration) the last term in the EOMs in this section can be
ie
dropped.
Note that the effects seen in Eq. (103), i.e.,
ev
• centrifugal pseudo (inertia) force: −M ω
e ω
bd bd (x + c f ),
e
e˙ (x + c ) ,
• Euler pseudo (inertia) force: −M ω bd f
r
• body pseudo (inertia) force: −Mφ A> τ̈,
are also available in commercial FE software, e.g., as loads, such as ANSYS and Abaqus, however, according
er
to the authors’ knowledge, the calculation is not based on the assembled mass matrix but rather performed at an
elemental level using shape functions.
pe
3 Conclusions
This paper presented a set of nodal-based floating frame formulation (FFRF) equations of motion stemming from
flexible multibody dynamics tailored for rotordynamics applications characterized by large rotations around a sin-
gle spatially fixed axis. In addition to this fundamental assumption, further constraints, such as those concerning
the rigid body translation of the floating frame, are directly integrated into the equations of motion. By gradually
ot
incorporating additional assumptions – reducing complexity and increasing efficiency – the well-established rotat-
ing frame rotordynamics formulation emerges from the FFRF. In summary, the following constraints are directly
enshrined into the equations of motion:
tn
The presented formulations are well-suited either for transient analysis, without necessitating 3D rotations
and differential algebraic equations, or for linearized dynamic analyses, including modal, harmonic, and stability
assessments. They are compatible with any displacement-based finite element and are presented both with and
without component mode synthesis reduction. Moreover, all necessary terms within the equations of motion, such
Pr
as spin softening and Coriolis effects, can be directly computed from the assembled finite element mass and
stiffness matrices, as well as nodal reference coordinates.
Follow up papers will employ the formulations for the dynamic analysis of a large bore two stroke diesel en-
gine and a highly transient optical application. Furthermore, future work will focus on the inclusion of centrifugal
stiffening in linearized dynamic analyses proposed in Sect. 2.5.
11
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
References
[1] J. S. Rao, History of Rotating Machinery Dynamics. Springer Netherlands, 2011.
d
[2] M. I. Friswell, J. E. T. Penny, S. D. Garvey, and A. W. Lees, Dynamics of Rotating Machines. Cambridge
University Press, 2010.
we
[3] E. Krämer, Dynamics of Rotors and Foundations. Springer Berlin, Heidelberg, 1993.
[4] G. Genta, Dynamics of Rotating Systems. Springer, New York, NY, USA, 2005.
[5] G. Genta and M. Silvagni, “On centrifugal softening in finite element method rotordynamics,” Journal of Ap-
ie
plied Mechanics, vol. 81, no. 1, 2013.
[6] S. Hou, R. Lin, L. Hou, and Y. Chen, “Dynamic characteristics of a dual-rotor system with parallel non-
concentricity caused by inter-shaft bearing positioning deviation,” Mechanism and Machine Theory, vol. 184,
ev
p. 105262, 2023.
[7] H. Guan, Q. Xiong, H. Ma, Y. Yang, J. Zeng, P. Wang, and Y. Bao, “Study on dynamic characteristics of the
gear-dual-rotor system with multi-position rubbing,” Mechanism and Machine Theory, vol. 191, p. 105501,
2024.
r
[8] J. Oh, A. Palazzolo, and L. Hu, “Stability of Non-Axisymmetric Rotor and Bearing Systems Modeled With
Three-Dimensional-Solid Finite Elements,” Journal of Vibration and Acoustics, vol. 142, no. 1, p. 011010,
2019.
er
[9] D. Kumar, “Rotordynamic analysis using 3d elements in fixed and rotating reference frame,” in Proceed-
ings of the ASME 2016 International Mechanical Engineering Congress and Exposition, vol. 1: Advances in
pe
Aerospace Technology, p. V001T03A017, 2016.
[10] D. Kumar, K. Juethner, and Y. Fournier, “Efficient Rotordynamic Analysis Using the Superelement Approach
for an Aircraft Engine,” in Proceedings of the ASME Turbo Expo 2017: Turbomachinery Technical Conference
and Exposition, vol. 7B: Structures and Dynamics, p. V07BT33A009, 2017.
ot
[11] M. S. Kumar, “Rotor dynamic analysis using ANSYS,” in IUTAM Symposium on Emerging Trends in Rotor
Dynamics. IUTAM Bookseries (K. Gupta, ed.), vol. 1011, pp. 153–162, Springer, Dordrecht, 2011.
[12] M. B. Wagner, A. Younan, P. Allaire, and R. Cogill, “Model reduction methods for rotor dynamic analysis: A
tn
survey and review,” International Journal of Rotating Machinery, vol. 2010, pp. 1–17, 2010.
[13] B. Besselink, U. Tabak, A. Lutowska, N. van de Wouw, H. Nijmeijer, D. Rixen, M. Hochstenbach, and
W. Schilders, “A comparison of model reduction techniques from structural dynamics, numerical mathematics
and systems and control,” Journal of Sound and Vibration, vol. 332, no. 19, pp. 4403–4422, 2013.
rin
[14] D. F. Li and E. J. Gunter, “Component Mode Synthesis of Large Rotor Systems,” Journal of Engineering for
Power, vol. 104, no. 3, pp. 552–560, 1982.
[15] L. B. Saint Martin, R. U. Mendes, and K. L. Cavalca, “Model reduction and dynamic matrices extraction
from state-space representation applied to rotating machines,” Mechanism and Machine Theory, vol. 149,
ep
p. 103804, 2020.
[16] T. Wasfy and A. Noor, “Computational strategies for flexible multibody systems,” Appl. Mech. Rev., vol. 56,
no. 6, pp. 553–613, 2003.
Pr
[17] A. A. Shabana, “Flexible multibody dynamics: Review of past and recent developments.,” Multibody System
Dynamics, vol. 1, pp. 189–222, 1997.
[18] A. Zwölfer and J. Gerstmayr, “A concise nodal-based derivation of the floating frame of reference formulation
for displacement-based solid finite elements: Avoiding inertia shape integrals,” Multibody System Dynamics,
vol. 49, pp. 291–313, 2020.
12
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
[19] A. Zwölfer and J. Gerstmayr, “The nodal-based floating frame of reference formulation with modal reduction:
How to calculate the invariants without a lumped mass approximation,” Acta Mechanica, vol. 232, pp. 835–
851, 2021.
d
[20] A. Zwölfer and J. Gerstmayr, “A unified framework for corotational flexible multibody system dynamics formu-
lations.,” Journal of Structural Dynamics, vol. 2, pp. 51–81, 2023.
we
[21] J. Gerstmayr, “Exudyn – a C++-based Python package for flexible multibody systems,” Multibody System
Dynamics, vol. 60, pp. 533–561, 2024.
[22] Y. Briend, M. Dakel, E. Chatelet, M.-A. Andrianoely, R. Dufour, and S. Baudin, “Effect of multi-frequency
parametric excitations on the dynamics of on-board rotor-bearing systems,” Mechanism and Machine Theory,
ie
vol. 145, p. 103660, 2020.
[23] A. Pechstein, D. Reischl, and J. Gerstmayr, “The applicability of the floating-frame based component mode
synthesis to high-speed rotors,” in Proceedings of the ASME 2011 International Design Engineering Technical
ev
Conferences and Computers and Information in Engineering Conference, vol. 4: 8th International Conference
on Multibody Systems, Nonlinear Dynamics, and Control, Parts A and B, pp. 933–942, 2011.
[24] H. Irschik, M. Nader, M. Stangl, and H.-G. v. Garssen, “A Floating-Frame-of-Reference Formulation for De-
formable Rotors Using the Properties of Free Elastic Vibration Modes,” in Proceedings of the ASME 2009
r
International Design Engineering Technical Conferences and Computers and Information in Engineering Con-
ference, vol. 4: 7th International Conference on Multibody Systems, Nonlinear Dynamics, and Control, Parts
er
A, B and C, pp. 881–888, 2009.
[25] M. Géradin and V. Sonneville, “A two-field approach to multibody dynamics of rotating flexible bodies,” Multi-
body System Dynamics, vol. 60, pp. 375–415, 2024.
pe
[26] J. Sopanen, Studies of rotordynamics using a multibody simulation approach. PhD thesis, Lappeenranta
University of Technology, 2004.
[27] A. Zwölfer and J. Gerstmayr, “Co-rotational formulations for 3D flexible multibody systems: A nodal-based
approach,” in Contributions to Advanced Dynamics and Continuum Mechanics (H. Altenbach, H. Irschik, and
V. P. Matveenko, eds.), pp. 243–263, Cham: Springer International Publishing, 2019.
ot
[29] F. Trainotti, A. Zwoelfer, J. Westphal, and D. J. Rixen, “Rotordynamics Continuum Finite Element Formulations
tn
from a Structural and Multibody Dynamics Perspective,” in IMAC XLII, (Orlando, FL, USA), 2024.
A Gravitational force
rin
The external nodal forces f may include gravity, which can be conveniently calculated as
Z
T
fg = ρ0 S g dV0 (110)
V0
Z
T
= ρ0 S SdV0 Φt g
ep
(111)
V0
= MΦt g (112)
= Mφ g (113)
T
Pr
= Mφ A g , (114)
since SΦt = I (partition of unity), where S is the FE shape function matrix, ρ0 is the body’s density and g is the
gravitational acceleration vector. Note that Mφ is constant and can be precomputed too. Furthermore,
f g = Mφ g , (115)
13
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
since Φt A> = A> Φ and M = Abd M A>
bd t bd
= M due to the fact that the consistent mass matrix and the rigid body
translation modes are composed out of blocks of the identity matrix .
d
B Simple example
we
Some of the to the proposed formulations’ specific matrices and a Campbell plot are shown in this section for the
simple example shown in Fig. 1 to foster understanding and help those interested to check their implementation.
ie
k2z
k1z 2
2 k2 y /2 k2x
k1 y /2 kc
z m2
ev
k2 y /2
y x m1 k2z
k1x k1 y /2
k1z 2
2
r
Figure 1: Hollow inertia-less cylinder rotating around the x axis; masses m1 and m2 are located inside the cylinder and are
er
attached with linear springs ki j (i = 1, 2, j = x, y, z) to its surface. A coupling spring kc acting in the x direction connects
the masses to one another.
The local (cylinder-fixed) coordinate system in Fig. 1 is located at the left base of the cylinder. The masses are
pe
located on the x axis at distances of 1 and 7/4 from the origin of the local coordinate system in their undeformed
configurations, respectively. Hence,
e.g. >
x = 1 0 0 7/4 0 0 , (116)
0 0 0
0 0 −1
e.g. 0 1 0
xe = . (117)
tn
0 0 0
0 0 −7/4
0 7/4 0
The cylinder rotates around the x axis. Therefore, according to Eq. (53), the rotation axis is defined by
rin
e.g. >
e = 1 0 0 (118)
e.g.
e = 0 0 −1 ,
e (119)
0 1 0
0 0 0 0 0 0
0 0 −1 0 0 0
e.g. 0 1 0 0 0 0
e bd =
e . (120)
0 0 0 0 0 0
0 0 0 0 0 −1
0 0 0 0 1 0
14
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
Then the plane orthogonal to the rotation axis is spanned by the columns of
0 0
e.g.
d
N = 1 0 , (121)
0 1
we
i.e., the unit vectors in y and z directions.
The mass matrix of the considered system shown in Fig. 1 reads
m1 0 0 0 0 0 1 0 0 0 0 0
0 m1 0 0 0 0 0 1 0 0 0 0
e.g. 0 0 m1 0 0 0 0 0 1 0 0 0
ie
M = = , (122)
0 0 0 m2 0 0 0 0 0 1 0 0
0 0 0 0 m2 0 0 0 0 0 1 0
0 0 0 0 0 m2 0 0 0 0 0 1
ev
and the local stiffness matrix, assuming that the springs in x, y, z are decoupled, is given by
r
0 0 k1z 0 0 0 0 0 100 0 0 0
e.g.
K = = 104 × . (123)
−kc 0 0 k2x + kc 0 er 0 −128 0 0 272 0 0
0 0 0 0 k2 y 0 0 0 0 0 196 0
0 0 0 0 0 k2z 0 0 0 0 0 100
Hence, the local elastic modes for no rotation speed are given by
pe
0 0 1 0 0 −1
0 0 0 0 1 0
e.g. 1 0 0 0 0 0
Ψ = , (124)
0 0 1 0 0 1
0 0 0 1 0 0
0 1 0 0 0 0
ot
and are associated with the angular eigenfrequencies from det K − ν2 M = 0
tn
e.g. >
ν0 = 1000 1000 1200 1400 1400 2000 . (125)
Hence,4
0 −1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
rin
1 0 0 0 0 0 0 0 −1 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 1 0 0 0 0 −1 0 0 0 −1 0
e e.g.
Ψ = , (126)
0 0 0 0 −1 0 0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 1 0 0 0 0 −1 0 0 0 0 0 0 0 0 −1
0 0 0 0 0 0 0 1 0 −1 0 0 0 0 0 0 1 0
ep
and
>
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Pr
e.g. 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
Iζ ⊗ e = . (127)
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
4
Note that in practical applications not all the elastic modes will be included in the reduction basis Ψ in order to achieve model reduction.
15
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498
Finally, a rigid body translation of the cylinder, and therefore the two masses, would be given by
1 0 0
d
0 1 0
e.g. 0 0 1
Φt = . (128)
1 0 0
we
0 1 0
0 0 1
e.g.
Writing Eq. (105) (no modal reduction) with Eqs. (120), (122), (123) and D = 0 in state space and solving
the homogeneous equation for multiple rotation speeds ϑ̇ according to [29] with Matlab’s eig(...) function,
ie
yields the Campbell diagram shown in Fig. 2. We see that the axial ( x ) modes are not affected by the rotation.
Furthermore, the y and z behaviors, i.e., forward and backward modes, are the same for both masses due to the
chosen stiffnesses. The angular frequencies at a rotation speed of 300 reads (eig(...))
ev
0.834497940852066
0.834497940852066
1.199999999910782
e.g.
ν300 = 103 × . (129)
1.563206060221641
r
1.563206060221641
2.000000000536606
er
pe
ot
tn
rin
ep
Figure 2: Campbell diagram (angular eigenfrequency ν over rotation speed ϑ̇) of the mechanical system depicted in Fig. 1.
Pr
16
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4830498