Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Dynamic modeling and simulations of rigid-flexible coupled systems using quaternion dynamics
(USC Thesis Other)
Dynamic modeling and simulations of rigid-flexible coupled systems using quaternion dynamics
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Dynamic Modeling and Simulations of Rigid-Flexible Coupled Systems Using
Quaternion Dynamics
by
Homin Choi
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MECHANICAL ENGINEERING)
June 2014
Copyright 2014 Homin Choi
ii
Table of Contents
List of Figures …………..……………………….…………….………..…………….iv
List of Tables …………………………………..…………………………….………vii
Nomenclature …………...…………………………………..…………….…..…….viii
Abstract …………………………………………………………………...…..………xi
Chapter 1 Introduction ……………..…………………………….………...………1
1.1 Motivation and Purpose …...……………………………..……………...…….1
1.2 Background ……………………..……………………………………….……2
1.3 Scope and Summary of Study …………………..………………………….…4
Chapter 2 Singularity of Rigid Body Dynamics Using Quaternion-Based Models
…………………………………………….……….………………………8
2.1 Conditions of Torque-Induced Singularity …………………………..………..8
2.2 Quaternion Model of a Rigid Body …………..…………………………...…13
2.2.1 Coordinate Transformation ……………………………….…….…15
2.2.2 Kinetic Energy ………………….………………………………….16
2.2.3 Potential Energy ……………….……………………………...…...17
2.2.4 Generalized Forces ………………………………………………...19
2.2.5 Equations of Motions ……………………………………………...23
2.3 Numerical Example ……………………….…………………………………24
2.3.1 Free Vibration (Case 1) …………….…………………..………….25
iii
2.3.2 Forced Vibration – Torque in the Dual Euler Basis (Case 2) …..….30
2.3.3 Forced Vibration – Torque in the Body Frame (Case 3) ……..……35
Chapter 3 Dynamic Modeling of Rigid-Flexible Coupled System ………...…….39
3.1 Description of the Model ………………….……………….………………...39
3.2 Governing Equation ………………….…………………………...……….…40
3.3 EOM of the Element of the Flexible Body ……….……………..…….…..…45
3.3.1 Nodal Coordinates and Interpolation of the Position ….…………..46
3.3.2 Mass Matrix ……..…………………………………...…………….48
3.3.3 Gyroscopic Matrix ………….………………………………….......53
3.3.4 Stiffness Matrix ………………………………………………..…..58
3.4 EOM of the Rigid Body …………..……………………………...……….….61
Chapter 4 Dynamic Analysis of Rigid-Flexible Coupled System ……….….........63
4.1 Free Vibration (Case 1) ……………….…………………………….........….66
4.2 Forced Vibration (Case 2) ……….……..………………………...……..……73
Chapter 5 Conclusion ……………….…………………………………………..….80
References …………………..…………………………….……….……………….…82
Appendices ………………….……………………………….…………………….....87
iv
List of Figures
Fig. 1 A sequence of Euler rotations …………………….……………………………9
Fig. 2 A constrained rigid body …………………….…..……………………………14
Fig. 3 A rotation sequence with Euler angles …………………..……………………15
Fig. 4 Deformation of the springs ……………………………..……………………..19
Fig. 5 Free vibration (Case 1): translation of the body ………………………………26
Fig. 6 Free vibration (Case 1): rotation of the body in terms of quaternion
components ……………………………………………………………………27
Fig. 7 Free vibration (Case 1): rotational displacements of the body in terms
of Euler angles ………………………………………………………………...28
Fig. 8 Free vibration (Case 1): rotational velocities of the body in terms of
Euler angles ……………………………………………………………………29
Fig. 9 Forced vibration – torque on
3
ˆ g (Case 2): translation of the body ……………31
Fig. 10 Forced vibration – torque on
3
ˆ g (Case 2): rotation of the body in
terms of quaternion components ……………………………………………..32
Fig. 11 Forced vibration – torque on
3
ˆ g (Case 2): rotational displacements of
the body in terms of Euler angles ……………………………………………33
Fig. 12 Forced vibration – torque on
3
ˆ g (Case 2): rotational velocities of the
body in terms of Euler angles ………………………………………………..34
Fig. 13 Forced vibration – torque on
3
ˆ e (Case 3): translation of the body ………….35
v
Fig. 14 Forced vibration – torque on
3
ˆ e (Case 3): rotation of the body in
terms of quaternion components ……………………………………………..36
Fig. 15 Forced vibration – torque on
3
ˆ e (Case 3): rotational displacements of
the body in terms of Euler angles ……………………………………………37
Fig. 16 Forced vibration – torque on
3
ˆ e (Case 3): rotational velocities of the
body in terms of Euler angles ………………………………………………..38
Fig. 17 Elements and nodes of the satellite ………………………………………….40
Fig. 18 Two-nodal beam element ……………………………………………………46
Fig. 19 Profile of the external torque ………………………………………………..65
Fig. 20 Free vibration (Case 1): translational displacements of the body ………...…67
Fig. 21 Free vibration (Case 1): translational velocities of the body ………………...68
Fig. 22 Free vibration (Case 1): rotational displacements of the body in terms
of quaternion components ……………………………………………………69
Fig. 23 Free vibration (Case 1): rotational velocities of the body in terms of
quaternion components ………………………………………………………70
Fig. 24 Free vibration (Case 1): rotational displacements of the body in terms
of Euler angles ……………………………………………………………….71
Fig. 25 Free vibration (Case 1): rotational velocities of the body in terms of
Euler angles ………………………………………………………………….72
Fig. 26 Forced vibration (Case 2): translational displacements of the body ………...74
Fig. 27 Forced vibration (Case 2): translational velocities of the body ……………..75
Fig. 28 Forced vibration (Case 2): rotational displacements of the body in
terms of quaternion components ……………………………………………..76
vi
Fig. 29 Forced vibration (Case 2): rotational velocities of the body in terms
of quaternion components ……………………………………………………77
Fig. 30 Forced vibration (Case 2): rotational displacements of the body in
terms of Euler angles …………………………………………………………78
Fig. 31 Forced vibration (Case 2): rotational velocities of the body in terms
of Euler angles ……………………………………………………………….79
vii
List of Tables
Table 1 Physical parameters of beam elements ……………………………………..64
Table 2 Physical parameters of the rigid body ………………………………………64
viii
Nomenclature
{ } = inertial frame
{ } = body frame
= displacement vector
= velocity vector
= angular velocity vector
= mass moment of inertia
= mass
= stiffness
= coefficient of
N
E
I
m
k
c
r
v
ω
viscous damping
= coefficient of rotating damping
= applied force
= applied torque
= kinetic energy
= potential energy
= work
= Lagrangian
= Lagrange multiplier
r
c
f
T
V
W
L
τ
λ
ix
= constraint
= alternating symbol
lmn
Γ
ε
Subscripts
= origin
= geometriccenter
= mass center
O
C
G
Superscripts
*
= vector
= unit vector
= quaternion
= complex conjugate of quaternion
= first derivative with respect to time
= second derivative with respect to time
= inertial frame
N
E
.
..
= body frame
x
Operators
= quaternion multiplication ⊗
xi
Abstract
In this study, a three-dimensional dynamic modeling and simulation of a satellite that
consists of coupled rigid-flexible bodies is presented. In the modeling of this coupled
system, the satellite is assumed to be a rigid body; and solar arrays, antennas or booms
are assumed to be flexible bodies. Coupling the equations of motions of the rigid body
and the flexible bodies is accomplished by using the finite element method (FEM). Since
the satellite rotates in space, the conventional FEM has a singularity problem caused by
rotational sequences of Euler angles. To solve this problem, quaternion-based FEM is
used in this study. Quaternion is well known as its singularity-free representation of
bodies in three-dimensional motions. The governing equations of motions are derived
from Lagrange equation, and derivatives in the equations are explicitly solved for
accurate results.
Although quaternions are singularity-free in modeling and analysis of rigid bodies in
three-dimensional motion, description of torques may lead to an unbounded response of a
quaternion-based model. In this study, theorems on the conditions of torque-induced
singularity in four coordinate systems – inertial frame, body frame, Euler basis, and dual
Euler basis – are derived. According to the theorems, torques applied in an inertial frame
or a body frame or an Euler basis will never cause unbounded motion. Torques applied in
a dual Euler basis, however, may lead to unbounded motion.
xii
Numerical simulations of free and forced vibrations are conducted for the quaternion-
based dynamic models of one rigid body and rigid-flexible coupled bodies. For the
dynamic model of one rigid body, the simulation results show that there is no singularity
in free vibrations, but when external torque is applied in dual Euler bases, singularity
occurs. However, when the torque is applied in the other coordinate systems, there is no
singularity. For the dynamic model of rigid-flexible coupled bodies, torque is applied
only in the body frame in this simulation, so both free and forced vibrations have no
singularity.
1
Chapter 1
Introduction
1.1 Motivation and Purpose
The multi-body systems that consist of rigid bodies and flexible bodies are commonly
found in many mechanical systems such as spacecrafts, robot arms and radars. In such
systems, the motions of rigid bodies and flexible bodies are coupled, and a lot of studies
have been done for these rigid-flexible coupled systems [1-3]. However, most of the
studies are limited to two-dimensional motions. Since space structures, such as satellites,
are subjected to three-dimensional motions, the two-dimensional studies are not proper to
be applied in such a case. Thus, in this study, the three-dimensional dynamic model of
the satellite consisting of rigid-flexible coupled bodies is developed, and three-
dimensional motions of the coupled system are simulated.
One of the important differences between two-dimensional motions and three-
dimensional motions is the representation of orientations of bodies. Unlike two-
dimensional motions, three-dimensional motions need rotational sequences to prescribe
orientations when Euler angles are used. As it is well known, the rotational sequence
makes a singularity problem, which causes significant error when controlling the attitude
of the satellite. This singularity problem can be solved by using quaternion. Quaternions
2
or Euler parameters form a homogenous set of four numbers which can remove kinematic
singularity. The purpose of this study lies in developing a singularity-free dynamic model
of a rigid-flexible coupled system by using quaternion dynamics, and simulating the
system in three dimensions.
The finite element method (FEM) is used for the flexible bodies of the system.
However, the conventional FEM uses Euler angles for the orientations, and usually is
applied to small deformations. Therefore, it is necessary to develop a quaternion-based
FEM for this study. The quaternion-based FEM not only can avoid singularity problems
but also can be applied to the large deformations.
1.2 Background
A lot of studies have been done on the rigid-flexible coupled systems. Huang et al [1]
applied symplectic geometric algorithms (SGA) to the rotating rigid-flexible coupled
system. SGA can preserve the features of Hamiltonian systems, which become a
dissipative system when integrated by an ordinary numerical method. Bai et al [2]
developed a dynamic equation of a satellite with flexible appendages based on the
coupling deformation field. In this study, the second-order coupling term of axial
displacement caused by transverse is included. Zhang et al [3] analyzed a rigid-flexible
coupled system by establishing the equations of motions according to the flexible multi-
body system’s way, and reduced the system equation dimensions by an introduction of
transition matrixes.
3
In 1843, Hamilton first introduced quaternions to express complex numbers in a
three-dimensional space [4]. Almost 60 years later, Gibbs devised vector analysis [5].
Because it is intuitive and easy to use, vector analysis is a more popular mathematical
tool than quaternions in science and engineering. In modeling the rotational motion of a
rigid body, vector analysis uses Euler angles in a sequence of three independent rotations
[6, 7]. However, because the Euler rotation sequence can cause singularity in dynamic
response of rigid bodies [8], quaternions have been revisited and adopted by many
researchers.
The major advantage of quaternions is that an Euler rotation sequence and related
ambiguity can be eliminated by a simple representation that uniquely describes the
rotation of a body. For instance, quaternions can be used to avoid a gimbal lock that is
related to Euler rotations [9]. Quaternions have been applied to attitude dynamics [10]
and kinematic analysis of robotic systems [11, 12], and more recently to three-
dimensional animation [13-15], signal processing [16, 17], and molecular dynamics [18-
20].
In 1941, the FEM was developed by Hrennikoff for the first time in the field of
structural engineering [21], and it has been improved by many researchers. In 1983,
Shabana and Wehage [22] represented flexible bodies with combined sets of reference
and local elastic generalized coordinates. In this paper, equations of motion and
constraints are formulated in terms of minimum mixed sets of modal and reference
generalized coordinates. Wallrapp and Schwertassek [23] included geometric stiffening
in multibody system simulation, and Meijaard [24] presented a spatial beam element for
static and dynamic problems which involve large displacements and rotations. In a recent
4
paper [25], a quaternion-based special beam element has been developed with Euler-
Bernoulli assumptions.
1.3 Scope and Summary of Study
In Chapter 2, conditions on torque-induced singularity of quaternion-based models
are derived for four coordinate systems. In study of the motion of a rigid body subject to
torques, four coordinate systems are available for torque description: inertial frame, body
frame, Euler basis, and dual Euler basis. Selection of a particular coordinate system
depends on the physics of the problem in consideration, and sometimes on the
investigator’s preference. It is well known that a quaternion-based model of the rigid
body can totally avoid singularities caused by establishment of equations of motion
directly from Euler angles, such as those in Euler’s equations. However, a quaternion-
based model still can have unbounded response due to torque description in a selected
coordinate system. This singularity is different from that of Euler’s equations in both
format and cause. The torque-induced singularity is caused by unbounded generalized
forces that appear on the right-hand side of an equation of motion; the singularity of
Euler’s equations or alike, on the other hand, is characterized by unbounded coefficients
contained in the left-hand side terms of an equation of motion. As shown in [26],
transformations between quaternions and Euler angles alone cause no singularity. In other
words, unbounded response of a quaternion-based model is solely caused by torque
description, but not coordinate transformation. With this understanding, it is necessary to
extend the singularity analysis to all the coordinate systems.
5
The current investigation is focused on the torque-induced singularity. This chapter
consists of three parts. In the first part, through investigation of generalized forces,
conditions on torque-induced singularity are derived for all of the four coordinate systems.
It is found that torque description in a dual Euler basis can render singularity in the
response of a quaternion-based model of rigid bodies, and that torque description in any
other coordinate will never cause an unbounded response of the quaternion model.
In the second part, a quaternion-based dynamic model of a rigid body is developed. In
the third part, numerical examples of free vibration and forced vibration are presented.
For the forced vibration, two cases are simulated. In the first case, torque is applied in the
body coordinates; and in the second case, torque is applied in the dual Euler basis. The
numerical simulations prove that free vibration and the torque applied in the body
coordinate have no singularity, but the torque applied in the dual Euler basis has
singularity.
In Chapter 3, the dynamic model of a rigid-flexible coupled system is developed
based on quaternion dynamics. The satellite model is used as the example. The satellite
body can be seen as a rigid body; and antennas, solar panels or booms can be seen as
flexible bodies. Equations of motions (EOM) of the flexible beams and the rigid body are
derived by Lagrange equation. Since FEM is used for modeling of the flexible beam,
coupling of the EOM of the flexible bodies and the rigid body is accomplished by
superposition. Since there is no relative motion between the rigid body and the affixed
end of the flexible body, it is feasible for them to share nodal coordinates, and equations
of motions can be coupled together.
6
For the flexible beams, quaternion-based FEM developed by Zhao and Ren [25] is
utilized; yet it is improved here to increase accuracy. Since the four components of
quaternions are not independent, the governing differential equation will have singular
mass matrix. To solve this problem, the paper [25] uses the Backward Differential
formulation (BDE) method [27] to integrate the governing differential equation.
However, in this study, an augmented form of the governing differential equation [28,
29] is used to avoid singular matrices. Since an augmented form of EOM has no singular
matrices, it can be simply integrated by the Runge-Kutta Method without singularity.
Another improvement is numerical differentiation. Numerical differentiation is usually
not as accurate as other numerical methods since it differentiates the approximated
interpolation polynomial. Thus, instead of numerical differentiation, all the
differentiations in the governing equation are explicitly solved in this study to increase
accuracy.
In Chapter 4, dynamic analysis of a quaternion-based rigid-flexible coupled system is
conducted. The dynamic model derived from the previous chapter is simulated in three
dimensions in this chapter. Two types of numerical simulations are performed here: free
and forced vibrations. Even if the quaternion can avoid singularity, it doesn’t give us
clear intuition of the orientations of the bodies. Thus, in the simulation, the initial values
are given in Euler angles, and they are converted to quaternions for the singularity-free
simulations. Then, after the simulations, they are converted to Euler angles again for the
final presentation of the simulation.
As it was mentioned earlier, the conversion itself between Euler angles and
quaternions makes no singularity problems. On the forced vibration, the torque is applied
7
in the body coordinate of the rigid body, so the rotational axis of the torque rotates
together with the rigid body of the satellite. The results of the simulations show no
singularity on both free and forced vibrations as the theorem in Chapter 2 prescribed.
8
Chapter 2
Singularity of Rigid-Body Dynamics Using
Quaternion-Based Models
2.1 Conditions of Torque-Induced Singularity
In study of the three-dimensional motion of a rigid body, there are four sets of base
vectors: a set of orthogonal unit vectors of the inertial frame { }
12 3
ˆˆ ˆ ,, nn n ;
a set of
orthogonal unit vectors of the body frame
{ }
12 3
ˆˆ ˆ ,, ee e ; a set of non-orthogonal unit vectors
associated with Euler rotations
{ }
12 3
ˆˆ ˆ , , gg g , which is also known as Euler basis; and the
dual Euler basis
{ }
12 3
ˆˆ ˆ , , gg g [30], which satisfies the conditions
jj
ii
ˆˆ ⋅= δ gg , with
j
i
δ being
the Kronecker delta. The
12 3
ˆˆ ˆ , , gg g in general are non-orthogonal, and they don’t have to
be unit vectors. Let ψ, φ and θ be the yawing, pitching and rolling angles of a rigid body
in rotation, respectively, as shown in Fig. 1. A virtual rotation of the body then can be
expressed as
1 2 3
ˆˆ ˆ δψ + δφ + δθ gg g .
9
Fig. 1 A sequence of Euler rotations: (a) yawing ψ; (b) pitching φ; (c) rolling θ
In quaternion-based modeling, a torque can cause unbounded motion of a rigid body.
This torque-induced singularity is related to the coordinate system chosen for torque
description. For instance, a torque applied in the direction of
3
ˆ g
can render the response
of a rigid body unbounded. In this section, we examine the torque-induced singularity in
quaternion-based modeling of rigid bodies. The analysis will consider torque description
in four coordinate systems: inertial frame, body frame, Euler basis, and dual Euler basis.
Investigation of the boundedness of generalized forces for Lagrange’s equation of motion
yields singularity conditions.
A torque can be either an externally applied load, or an internal one (like
damping/spring torque). The virtual work by a torque τ that is described in any
coordinate system can be expressed as
0 0 11 2 2 3 3
W Q p Q pQ p Q p
τ
δ = ⋅δ + ⋅δ + ⋅δ + ⋅δ
(1)
10
Here
01 2 3
, , and p pp p are the components of the unit quaternion, with the constraint
2222
01 2 3
1 pppp +++ = ;
01 2 3
, , , Q QQ Q are generalized forces, which
can be expressed as
[31]
0
1
1
2
2
3
3
ˆ
ˆ
ˆ
T
Q
Q
Q
Q
⋅
= ⋅
⋅
τn
A τn
τn
(2)
with { }
12 3
ˆˆ ˆ ,, nn n being a set of base vectors of an inertial frame, and
1 0 32
23 0 1
3 21 0
2
p p p p
pp p p
p pp p
−−
= −−
−−
A (3)
The generalized inverses of matrix A is
1
, with
4
T−−
= = A A AA I (4)
where I is the identity matrix [31]. It follows from Eqs. (2) and (4) that
2
2222
01 2 3
4 QQQQ +++ = τ
(5)
So, the generalized forces are bounded if the torque is finite.
Theorem 1
For a rigid body subject to a torque that is specified in an inertial frame,
11 2 2 3 3
ˆˆ ˆ
nn n
= τ + τ + τ nn n τ (6)
where
12 3
, ,
nn n
ττ τ are finite quantities, the dynamic response of a quaternion-based
model of the body has no singularity.
11
Proof: Because
2
22 2
12 3 nn n
= τ + τ + τ τ , which is finite, the generalized forces by Eq. (5)
are always bounded. So the response of the quaternion-based model has no singularity.
Theorem 2
For a rigid body subject to a torque that is specified in a body frame,
11 2 2 3 3
ˆˆ ˆ
ee e
= τ + τ + τ ee e τ (7)
where
12 3
,,
e e e
ττ τ are finite quantities, the dynamic response of a quaternion-based model
of the body has no singularity.
The proof of Theorem 2 follows that of Theorem 1.
Theorem 3
For a rigid body subject to a torque that is specified in an Euler basis
11 2 2 3 3
ˆˆ ˆ
gg g
= τ + τ + τ gg g τ (8)
where
1 2 3
,,
gg g
ττ τ are finite quantities, the dynamic response of a quaternion-based
model of the body has no singularity.
Proof: Because
12 3
ˆˆ ˆ , , gg g
are unit vectors, the torque given in Eq. (8) is finite. By Eq.
(5), the generalized forces are always bounded, indicating that the dynamic response of
the quaternion-based model has no singularity.
Theorem 4
12
For a rigid body subject to a torque that is described in a dual Euler basis
{ }
12 3
ˆˆ ˆ , , gg g ,
a quaternion-based model of the body has the following two cases of singularity:
(C1) If the applied torque is in line with the second base vector
2
ˆ g , the dynamic
response of the quaternion-based model is always bounded.
(C2) If the torque is in line with either the first base vector
1
ˆ g or the third base vector
3
ˆ g , the dynamic response of the quaternion-based model becomes unbounded.
Proof: A torque applied in line with one of the dual Euler base vectors is expressed as
ˆ ,
i
= τg τ
where i = 1 or 2 or 3 and τ is a finite quantity. By Eq. (5),
2
2222 2
01 2 3
ˆ 4
i
QQQQ +++ =τ g
(9)
The magnitudes of
13
ˆˆ and gg vary, depending on the rotation of the rigid body. To see
this, denote a sequence of Euler rotations of the body. Without loss of generality,
consider a set of asymmetric Euler angles,
12 3
ˆˆ ˆ ψ +φ +θ gg g
as shown in Fig. 1. By the
conditions ˆˆ
j j
ii
⋅= δ gg ,
11
11
22
22
33
33
ˆˆ ˆ ˆ cos 1
ˆ ˆ ˆ ˆ 1
ˆˆ ˆ ˆ cos 1
⋅ = φ=
⋅= =
⋅ = φ=
g g g g
g g g g
gg g g
(10)
where φ is the second Euler angle (pitching angle in Fig. 1), and
2
2
ˆˆ = gg has been used.
Because
12 3
ˆˆ ˆ , , gg g are unit vectors, one obtains from Eq. (10) that
1 23
11
ˆ ˆˆ , 1,
cos cos
= = =
φφ
g gg (11)
13
When φ approaches
2
π
± ,
1
ˆ g and
3
ˆ g in Eq. (11) become infinite. It follows from Eq.
(9) that the generalized forces are unbounded, which leads to the result in (C2). On the
other hand,
2
ˆ 1 = g for all values of the Euler angles, which by Eq. (9) yields the result in
(C1). For a set of symmetric Euler angles (such as rotation sequences 1-2-1, 2-3-2 and 3-
1-3), it can be shown that
1
ˆ g and
3
ˆ g go to infinity as φ goes to 0 or π , and that
2
ˆ g
is always bounded for any values of Euler angles.
Discussion
Consider the Euler basis in Fig. 1. Since
1 2
ˆˆ ⊥ gg and
2 3
ˆ ˆ ⊥ g g ,
12 3
ˆˆ ˆ , , gg g form a base
when 90 φ≠ ± ° . However,
12 3
ˆˆ ˆ , , gg g become linearly dependent when 90 φ= ± ° . In other
words, depending on a rotation sequence, at a particular value of φ ,
12 3
ˆˆ ˆ , , gg g become
linearly dependent, and they cannot form a basis in the 3-D space. Because of this, the
tensor transformation by
jj
ii
ˆˆ ⋅= δ gg cannot be performed, and it can cause singularity
in
j
ˆ g .
2.2 Quaternion Model of a Rigid Body
Consider a rigid body constrained by four pairs of springs and dampers; see Fig. 2,
where τ is an external torque applied to the body, C and G is the geometric center and
the center of mass of the body, respectively. The body has six degrees of freedom: three
14
translations and three rotations. The translations can be expressed by u, v and w, which
are the displacement components of the mass center G of the body measured in an inertial
frame. The rotations of the body are conventionally expressed by Euler angles. Fig. 3
shows a typical sequence of Euler rotations from xyz to XYZ, which is a rotation about x-
axis (yawing ψ, from xyz to
1 11
x yz ), followed by a rotation about the new y-axis (pitching
φ, from
1 11
x yz to
2 22
x yz ), and followed by a newer z-axis (rolling θ, from
2 22
x yz to XYZ).
Fig. 2 A constrained rigid body
15
(a) Yawing ψ (b) Pitching φ (c) Rolling θ
Fig. 3 A rotation sequence with Euler angles
2.2.1 Coordinate Transformation
The rotations of the body can also represented by the unit quaternion
01 2 3 0
ˆˆ ˆ
pp p p p = ++ + = + p i jk p
, with the constraint
2222
01 2 3
() 1 0 pppp = +++−= p Γ .
The quaternion can be written as [32]
= ⊗⊗ pg h t
(12)
with
ˆˆ ˆ
cos sin , cos sin , cos sin
2 2 22 22
ψ ψ φφ θθ
=+ =+=+ g i h jt k
(13)
where vectors ,, gh t
describe a sequence of yawing ( ψ), pitching ( φ) and rolling (θ)
rotations, and
ˆ
i ,
ˆ
j ,
ˆ
k are the unit vectors of frame xyz before the rotations ( Fig. 3).
Define two coordinate systems: inertial frame {N} and body frame {E}; see Fig. 2.
Here, superscripts N and E denote quantities in the initial and body frames, respectively.
A vector prescribed in the body frame can be transformed to the inertial frame by unit
16
quaternion multiplication, or vice versa. For instance, a vector
E
r
in the body frame can
be transformed to the corresponding vector
N
r
in the inertial frame by
* NE
=⊗⊗ rp r p
(14)
where p is a unit quaternion representing the coordinate transformation,
*
p is the
complex conjugate of p , and ⊗ is the multiplication operator for quaternions.
As shown in Appendix A, the Euler angels (ψ, φ and θ ) and the quaternion
components ( p
0
, p
1
, p
2
and p
3
) are in the relation
[ ] []
T
p
= RP (15)
where the matrices
22
2 3 1 2 0 31 30 2
22
1 2 03 1 3 23 0 1
22
13 0 2 2 3 01 1 2
cos cos cos sin sin sin cos cos sin cos sin sin
[ ] sin cos sin sin sin cos co
12 2 2 2 2 2
[] 2 2 1 2 2 2 2
2 2 2 2 12 2
p
p p pp p p pp p p
p p pp p p p p pp
p p pp p p pp p p
θ φ θφ ψ θ ψ θφ ψ θ ψ
θ φ θ φψ θ
+ − +
= − − +
− − − +
= + − − −
− + − −
R
P
s sin sin cos cos sin
sin cos sin cos cos
ψ θ φ ψ θ ψ
φ φψ φ ψ
+
−
(16)
For instance, the rolling angle θ is related to the quaternion components by
0 3 12
22
23
22
tan
12 2
p p pp
pp
−
θ=
−−
(17)
2.2.2 Kinetic Energy
The displacement vector
N
OG
r
of the mass center G of the body by Eq. (14) is
17
* NN E
OG OC CG
= +⊗ ⊗ r r p r p
(18)
where
N
OC
r
is the displacement vector of the geometric center C , and
E
CG
r
is the vector
from C to G. The kinetic energy of the body is
{ } { }
2
11
22
T
N E EE
OG G
Tm = +
v ωI ω
(19)
where m is the mass of the body,
N
OG
v
is the velocity
of the mass center,
{ }
E
ω is the
angular velocity of the body, and
E
G
I is the inertia matrix of the body with respect to
the mass center G and in the body frame. The
N
OG
v
is obtained through differentiation of
the
N
OG
r
with respect to time in the inertial frame:
( )
N
NN
OG OG
d
dt
= v r
(20)
The angular velocity
{ }
E
ω , measured in the body frame, is given by [29, 33]
00
2( )
E
pp = − − × ω p p pp
(21)
2.2.3 Potential Energy
The potential energy of the body is produced during the deformation of the four linear
springs shown in Fig. 2. Without loss of generality, gravity is not considered herein. In
Fig. 4, denote the two ends of the body in equilibrium by
1
D and
2
D , at which the springs
and dampers are attached. Denote the positions of the two ends of the body in motion by
18
1
D ′ and
2
D ′ . The displacement vectors of the two ends of the body, denoted by
11 2 2
' '
and
NN
DD D D
rr
, are
11 1 1
22 2 2
*
''
*
''
N N E N
D D OC CD OD
NN E N
D D OC CD OD
= +⊗ ⊗ −
= +⊗ ⊗ −
r r pr p r
r r pr p r
(22)
The elongation of the i th spring is
1, ..., 4
ki ki ki
l ll i ′ ∆= − = (23)
where
ki
l is the initial length of the spring without deformation, and
ki
l ′ the length after
deformation. Write
11
' 11 2 2 3 3
ˆˆ ˆ + +
N
DD
= αα α r nn n
and
22
' 11 2 2 3 3
ˆˆ ˆ + +
N
DD
= ββ β r n n n
. It is from
Fig. 4 that
2 2 2 222
1 1 1 2 3 2 2 1 2 3
2 2 2 222
3 3 2 1 3 4 42 1 3
' ( ) ,' ( )
' ( ) ,' ( )
kk k k
kk k k
ll l l
ll l l
= + α + α + α = + β + β + β
= + α + α + α = + β + β + β
(24)
The potential energy of the body then is
22 22
1 1 22 3 3 44
11 1 1
22 2 2
k kkk
V k l kl k l kl = ∆ + ∆ + ∆ + ∆ (25)
19
Fig. 4 Deformation of the springs
2.2.4 Generalized Forces
Assume that an external torque τ
is applied on the body. Then the virtual work done
by the external torque τ
and damping forces is
d
WW W
τ
δδ δ = + (26)
with
1 1 22 3 3 44
()
d c k ck c k c k
cr
W f lfl f l fl
W
τ
δ δ δδδ
δδ
= ∆ +∆ + ∆ +∆
= +⋅ τ τξ
(27)
20
where
d
W δ and W
τ
δ are virtual works done by the dampers at the two ends of the body
and torque respectively,
cj
f are the damping forces exerted by the dampers, and
cr
τ
is a
damping torque against the rolling motion of the body caused by τ
. Viscous damping is
assumed:
, , 1,2,3,4
kj
cr r cj j
dl
d
c f c j
dt dt
∆
= −= − =
ξ
τ
(28)
with
r
c and
j
c being damping coefficients. Because virtual displacements
ki
l δ∆ and δ ξ
in Eq. (27) can be expressed in terms of the variations of the displacement parameters (u,
v, w) and quaternion parameters, one can write
1 2 3 4 0 51 6 2 7 3
W Q u Qv Q w Qp Q p Q p Qp δ = δ +δ + δ +δ + δ + δ +δ (29)
where the generalized forces
i
Q are functions of the quaternion components. From Eqs.
(23) and (24),
ki
l δ∆ is expressed as follows.
1 1 1 1 10 0 11 1 12 2 13 3
2 2 2 2 20 0 21 1 22 2 23 3
3 3 3 3 30 0 31 1 32 2 33 3
4 4 4 4 40 0 41 1 42 2 4
ku v w
ku v w
ku v w
ku v w
l uv w p p p p
l uv w p p p p
l uv w p p p p
l uv w p p p
δ γ δ γ δ γ δ γδ γ δ γδ γ δ
δ γ δ γ δ γ δ γδ γ δ γδ γ δ
δ γ δ γ δ γ δ γδ γ δ γδ γ δ
δ γ δ γ δ γ δ γδ γ δ γδ γ
∆= + + + + + +
∆= + + + + + +
∆ = + + + + + +
∆= + + + + + +
33
p δ
(30)
with,
11 2
11
222 222
11 2 3 11 2 3
3 1 1 2 21 3 0
1 10
222 222
11 2 3 11 2 3
11 3 2 0 11 0 2 3
11 12
222 22
11 2 3 11 2
,,
() ()
() 2
,,
() ()
() ()
,
() ()
k
uv
kk
k
w
kk
kk
kk
l
ll
ll p l p l p
ll
ll p l p ll p l p
l l
α α
γγ
α αα α αα
α αα α
γγ
α αα α αα
α α α α
γγ
α αα α αα
+
= =
+ ++ + ++
−+ + −
= =
+ ++ + ++
−+ + −+ −
= =
+ ++ + ++
2
3
,
21
1 1 1 2 2 33 2 1
13 2
222 2 2 2
11 2 3 2 1 2 3
3 2
22
222 222
21 2 3 21 2 3
21 2 2 1 3 0 21 3 2 0
20 21
222 2
21 2 3 21
() 2
,,
() ( )
,,
() ()
() 2 ()
,
() ()
kk
u
k k
vw
k k
kk
kk
ll p l p l p l
ll
ll
ll p l p l p ll p l p
l l
α α α β
γγ
α αα β β β
β β
γγ
β ββ β ββ
β β β β β
γγ
β ββ β
−+ − − +
= =
+ ++ + + +
= =
+ ++ + ++
+ −+ + −
= =
++ + ++
22
2 3
,
ββ +
21 0 2 3 21 1 2 2 3 3
22 23
222 222
21 2 3 21 2 3
32 1
33
222 222
32 1 3 32 1 3
() () 2
,
() ()
,,
() ()
k k
k k
k
uv
kk
ll p l p ll p l p l p
l l
l
ll
β β ββ β
γγ
β ββ β ββ
α α
γγ
α αα α αα
+ + ++ +
= =
+ ++ + ++
+
= =
+ ++ + ++
(31)
3 3 2 1 1 2 3 0
3 30
222 222
32 1 3 32 1 3
32 0 1 3 32 3 1 0
31 32
222 222
32 1 3 32 1 3
3 2 2 1 1 3 3
33 4
222
32 1 3
() 2
,,
() ()
() ()
,,
() ()
() 2
,
()
k
w
kk
kk
kk
k
u
k
ll p l p l p
ll
ll p l p ll p l p
l l
ll p l p l p
l
α αα α
γγ
α αα α αα
α α α α
γγ
α αα α αα
α αα
γγ
α αα
+ −−
= =
+ ++ + ++
+ − −+ −
= =
+ ++ + ++
−+ − −
= =
+ ++
1
222
42 1 3
,
()
k
l
β
β ββ + ++
42 3
44
222 222
42 1 3 42 1 3
42 1 1 2 3 0 42 0 1 3
40 41
222 222
42 1 3 42 1 3
42 3 1 0 42 2 1 1
42 43
222
42 1 3
,,
() ()
() 2 ()
,,
() ()
() ()
,
()
k
vw
k k
kk
kk
k k
k
l
ll
ll p l p l p ll p l p
ll
ll p l p ll p l p
l
ββ
γγ
β ββ β ββ
ββ β β β
γγ
β ββ β ββ
ββ β β
γγ
β ββ
+
= =
+ ++ + ++
−+ + + −+ +
= =
+ ++ + ++
+ + + ++
= =
+ ++
33
222
42 1 3
2
()
k
lp
l
β
β ββ + ++
Virtual displacement δ ξ
is expressed in terms of quaternions as follows from [29] and Eq.
(15).
11 2 2 3 3
11 2 2 3 3
12 3
ˆˆ ˆ
ˆˆ ˆ
ˆˆ ˆ
nn n
ee e
δ δξ δξ δξ
δξ δξ δξ
δψ δφ δθ
= ++
= + +
= ++
ξn n n
ee e
gg g
(32)
with
22
1 23 32 0 1 1 0
2 3 1 2 0 1 3 0 2
3 1 2 0 3 3 0 21
2( )
2( )
2( )
n
n
n
p p p p p p pp
p p p p pp p p
pp p p p p p p
δξ δ δ δ δ
δξ δ δ δ δ
δξ δ δ δ δ
= − + −
= − −+
= + −−
(33)
and
1 32 23 1 0 0 1
2 1 3 0 2 31 2 0
3 21 3 0 0 3 1 2
2( )
2( )
2( )
e
e
e
p p p p pp p p
pp p p p p p p
p p p p p p pp
δξ δδ δδ
δξ δ δ δ δ
δξ δδ δδ
= − −+
= + −−
= −+ −
(34)
and
22
1 12
0 2 22 2
1 2 2 3 01
22
0 0 2 01 1 2 3
1 2 22 2
1 2 2 3 01
22
3 3 2 31 0 1 2
2 2 22 2
1 2 2 3 01
22
2 12
22
12
2 (1 2 2 )
(1 2 2 ) ( 2 2 )
24 4 8
(1 2 2 ) ( 2 2 )
24 4 8
(1 2 2 ) ( 2 2 )
2 (1 2 2 )
(1 2 2
p pp
p
p p p p pp
p pp pp p p p
p
p p p p pp
p pp pp p p p
p
p p p p pp
p pp
pp
δψ δ
δ
δ
−−
=
− − +− +
−+ −
+
− − +− +
− − + +
+
− − +− +
−−
−
−−
3 22
2 3 01
31 1 3 2 0 0 2
22 22
1 3 01 2 3 0 2
22
3 23
0 2 22 2
2 3 12 0 3
22
2 23
1 2 22 2
2 3 12 0 3
22
1 1 2 13
) (2 2 )
2( )
14 8 4
2 (1 2 2 )
(1 2 2 ) ( 2 2 )
2 (1 2 2 )
(1 2 2 ) ( 2 2 )
24 4 8
p
p p pp
p p pp p p p p
pp p p p p pp
p pp
p
p p pp p p
p pp
p
p p pp p p
p pp pp
δ
δ δ δ δ
δφ
δθ δ
δ
+− +
+ + +
=
−− −
−−
=
− − +− +
−−
−
− − +− +
− − + +
+
0 2 3
2 2 22 2
2 3 12 0 3
22
0 0 2 0 3 1 2 3
3 2 22 2
2 3 12 0 3
(1 2 2 ) ( 2 2 )
24 4 8
(1 2 2 ) ( 2 2 )
pp p
p
p p pp p p
p pp pp p p p
p
p p pp p p
δ
δ
− − +− +
−+ −
+
− − +− +
(35)
where,
i
ˆ n ,
i i
ˆˆ and eg are unit vectors of inertial coordinate, body coordinate and Euler
basis respectively. In this study, torque is applied in the dual Euler basis and the body
23
frame for the simulations. For the first case, the axis of the external torque τ
is in line
with the unit vector
3
ˆ g of dual Euler basis. Then the virtual work W
τ
δ is
3
12 3
ˆ ˆˆ ˆ ( )( )
()
cr
cr
W
τ
δ τ τ δψ δφ δθ
τ τ δθ
= +⋅ + +
= +
g gg g
(36)
where,
cr
τ and τ are the absolute values of
cr
τ
and τ
respectively. For the second case,
the axis of the external torque τ
is in line with the unit vector
3
ˆ e of the body frame {E}.
In other words, the torque always follows the motion of the body. Then the virtual work
W
τ
δ is
3 11 2 2 3 3
3
ˆˆ ˆ ˆ ( )( )
()
cr e e e
cr e
W
τ
δ τ τ δξ δξ δξ
τ τ δξ
= +⋅ + +
= +
ee e e
(37)
2.2.5 Equations of Motions
To derive the equations of motion, define the seven-vector
{ } { } { }
12 7 0 1 2 3
TT
q q q uv w p p p p = = q (38)
With the constraint
2222
01 2 3
() 1 0 pppp = +++−= p Γ , Euler-Lagrange equations for the
body read [34]
, 1,2,...,7
n
nn n
dL L
Qn
dt q q q
∂∂ ∂
− −λ = =
∂∂ ∂
Γ
(39)
where the Lagrangian L TV = − , and λ is a Lagrange multiplier. Because the constraint
is holonomic, we can use the Pfaffian form () 0
d
dt
= p Γ [34], which in the matrix form is
24
( ) [ ( )]{ } 0
d
dt
= = p Aq q Γ (40)
with
[ ]
0 1 2 3
[ ( )] 000 p pp p = Aq . It follows that the equations of motion are
given in an augmented form [28]
{ ( , )} { } { ( )} { ( )} [ ( )] [ ( )]
() 00 [ ( )] 0
T
c
D
+ +=
− λ
Dq q q K q Qq M q Aq
q Aq
(41)
where [ ( )] M q is a mass matrix, { ( , )} Dq q describes gyroscopic forces, { ( )} Kq is
related to elastic restoring forces, { ( )} Qq is a generalized force vector, and [ ( )]
T
λ Aq is
due to the constraint (40), and
2222
01 2 3
()
c
D pppp = +++ q . Equation (41) can be solved
via numerical integration.
2.3 Numerical Example
The results presented in the previous sections are illustrated on the rigid body shown
in Fig. 2. For numerical simulation, the values of the system parameters are chosen as
follows: for the body,
3
density 2800 / kg m ρ= , length 0.3 l m = , radius 0.02 rm = ,
mass 1.056 m kg = , geometric center C located at the origin of the initial frame, mass
center G with coordinates 0.003 , 0.003 , 0.008
G GG
x m y mz m = = = , and the matrix of
inertia
42
72 0.532 1.419
0.532 72 1.419 10
1.419 1.419 2.856
E
G
kg m
−
− −
= − −× ⋅
− −
I ;
25
and for the supports, initial length of springs and dampers 1
i
lm = , spring coefficients
0.1 /
i
k Nm = , viscous damping coefficients 0.1 /
i
c Ns m = ⋅ , 1,2,3,4 i = .
Three cases of vibration are considered in the simulation:
Case 1. The body is in free vibration that is caused by initial rotational velocities
(0) ψ= 1 rad/s, (0) φ=
10 rad/s and (0) θ=
1 rad/s. All other displacement and velocity
parameters are initially zeros. In this case, the body has no torsional damping ( 0
r
c = ) .
Case 2. The body is in forced vibration that is excited by an external torque (Fig. 2) of
constant magnitude 0.1 Nm τ= ⋅ . The axis of the external torque τ
is in line with the unit
vector
3
ˆ g of the dual Euler basis. The body is initially at rest and has rotational damping
with 0.01 s/
r
c N m rad = ⋅⋅ .
Case 3. The body is in forced vibration that is excited by an external torque of constant
magnitude 0.1 Nm τ= ⋅ like Case 2. However in this case, the axis of the external torque
τ
is in line with the unit vector
3
ˆ e of the body frame. The body is initially at rest and has
rotational damping with 0.01 s/
r
c N m rad = ⋅⋅ .
2.3.1 Free Vibration (Case 1)
The dynamic response of the body in the free vibration (Case 1) is plotted in Figs. 5-8.
Figs. 5 and 6 shows that both the translation and rotation of the body are bounded. The
body rotation can also be expressed in terms of Euler angles by coordinate transformation
(15). Figure 7 shows that the rotational response of the body is always bounded even
26
when the pitching angle φ is near /2 ±π rad, at which values an Euler angle-based model
would have singularity. The free response is bounded because no torque is applied to the
body ( ( )0
cr
τ + τ δθ = ). The free vibration of the body could be unbounded if there exists
torsional damping
cr
τ . This has been observed in our numerical simulation.
Fig. 5 Free vibration (Case 1): translation of the body
27
Fig. 6 Free vibration (Case 1): rotation of the body in terms of quaternion components
28
Fig. 7 Free vibration (Case 1): rotational displacements of the body in terms
of Euler angles
29
Fig. 8 Free vibration (Case 1): rotational velocities of the body in terms of
Euler angles
30
2.3.2 Forced Vibration – Torque in the Dual Euler Basis (Case 2)
The forced vibration simulation (Case 2), however, exhibits singularity around t =
0.161 s; see Fig. 9-12. This is predicted by Theorem 4 in Section 2.1. It is commonly
believed that quaternion-based models of rigid bodies can totally eliminate the singularity
experienced by those models with Euler rotations. Our numerical simulation herein
reveals that singularity still occurs when torques are applied to the body. It is easy to see
from Eqs. (35) and (36) that δθ in Eq. (36) becomes infinite when the denominator in Eq.
(35) vanishes at particular values of quaternion components. This renders the generalized
forces
4 5 6 7
,, , Q QQ Q
in Eq. (29) unbounded at a finite time. This means that application
of a torque can lead to unbounded response of a quaternion-based model like Eq. (41).
31
Fig. 9 Forced vibration – torque on
3
ˆ g (Case 2): translation of the body
32
Fig. 10 Forced vibration – torque on
3
ˆ g (Case 2): rotation of the body in terms of
quaternion components
33
Fig. 11 Forced vibration – torque on
3
ˆ g (Case 2): rotational displacements
of the body in terms of Euler angles
34
Fig. 12 Forced vibration – torque on
3
ˆ g (Case 2): rotational velocities of the
body in terms of Euler angles
35
2.3.3 Forced Vibration – Torque in the Body Frame (Case 3)
The forced vibration simulation (Case 3) exhibits no singularity; see Fig. 13-16. Since
the torque is applied in the body frame, the dynamic response of the body is bounded
according Theorem 2 in Section 2.1. We can see this obviously from Eqs. (34) and (37).
Unlike Eq. (35), Eq. (34) has no denominators and is always bounded. Therefore, Eq.
(37) never goes infinite.
Fig. 13 Forced vibration – torque on
3
ˆ e (Case 3): translation of the body
36
Fig. 14 Forced vibration – torque on
3
ˆ e (Case 3): rotation of the body in terms of
quaternion components
37
Fig. 15 Forced vibration – torque on
3
ˆ e (Case 3): rotational displacements
of the body in terms of Euler angles
38
Fig. 16 Forced vibration – torque on
3
ˆ e (Case 3): rotational velocities of the
body in terms of Euler angles
39
Chapter 3
Dynamic Modeling of Rigid-Flexible Coupled
System
3.1 Description of the Model
The satellite model considered in this paper is shown in Fig. 17. It consists of two
parts: the satellite body that is assumed to be a rigid cube and the two antennas that are
assumed to be flexible beams with a square cross-section. The centerlines of each beam
and the center of the rigid body lie on the same line when there is no deformation. FEM
is used for dynamic modeling of this rigid and flexible bodies-combined system. Each
element of the flexible body has two nodes at the ends, and the rigid body has one node
as shown in Fig. 17.
Since the two beams are firmly fixed on the rigid body, we can assume that there is
no relative motion between the rigid body and the attached ends of the beams. Therefore,
they share one node, k as shown in Fig. 17. In Fig. 17, the italic numbers on the top of the
bodies and the bold numbers below the bodies represent the number of nodes and the
number of elements respectively, and nodes are laid on the centerline of the beam. The
40
last node number n should be an odd number for each beam to have the same numbers of
elements, and node k is the node in the middle.
Fig. 17 Elements and nodes of the satellite
3.2 Governing Equation
The governing equation of this system is accomplished by the following steps. First,
EOM of the beam elements and the rigid body are derived by a Lagrange equation.
Second, the EOM are added together using a method of superposition. And finally, the
augmented form of the EOM is built up for non-singular simulation. The full derivations
of EOM of the beam elements and the rigid body are given in Sections 3.3 and 3.4,
respectively, so we will focus on the steps, constraints, superposition and augmented
matrices in this section.
First of all, the nodal coordinate of the system should be defined before we find the
governing equation. The nodal coordinates on i th node is expressed as
41
{}
{ } { } , 1,2,...,
i
i i
i
in
ε
= =
r
qp (42)
Here, {}
i
r is the position vector of node i, {}
i
p is four-tuple quaternion that prescribes
the orientation of node i, and
i
ε is normal strain at neutral axes [25]. Thus the nodal
coordinates {}
i
q has 8 components.
Next, EOM of the elements of the flexible beam is to be obtained. The i th element of
the flexible body yields
1
[ ]{ } [ ]{ } [ ]{ } [ ] { } { },
{}
{ } , 1,2,..., 1
{}
T
ii ii ii i i i
i
i
i
in
+
++ + =
= = −
ee ee ee e e e
e
M q Vq K q A ΛF
q
q
q
(43)
where superscript i and subscript i represent node number and element number
respectively, []
i e
M is a mass matrix, [ ]
i e
V describes gyroscopic forces, []
i e
K is related
to elastic restoring forces, {}
i e
F is a external force or damping force vector, and
[ ]{ }
T
ii ee
A Λ is a term for constraints. Since quaternion has four elements for three DOF,
we have the following constraint for each quaternion.
22 2 2
01 2 3
( ) ( ) ( ) ( ) 1 0, 1,2,...,
ii i i i
p pp p i n Γ = + + + −= = (44)
Then, the constraint for each element yields
1
{ } , 1,2,..., 1
i
i
i
in
Γ
Γ
+
= = −
e
Γ (45)
By differentiating Eq. (4) with respect to {}
i e
q , we get [ ]
i e
A [34],
42
{}
[ ] , 1, 2,..., 1
{}
i
i
i
in
∂
= = −
∂
e
e
e
Γ
A
q
(46)
and {}
i e
Λ is a Lagrange multiplier vector for i th element [34].
1
{ } , 1, 2,..., 1
i
i
i
in
λ
λ
+
= = −
e
Λ (47)
Note that the right end of k-1 element and the left end of k element don’t coincide with
their nodal coordinates, {}
k
q . While the orientations and the normal strains of the three
points are the same, the positions are different. Therefore Eq. (43) for k-1 and k element
must be derived by the modified position vector {}
k ′
r as follows.
{ } { } { } for element
2
{ } { } { } for element
2
kk k
kk k
h
h
′
′
= −
= +
r r l k -1
rr l k
(48)
where {}
k
l is the unit normal vector of the cross-section of the beam at node k.
Next, EOM of the rigid body is expressed as
[ ]{ } [ ]{ } [ ] { }, ( 1) / 2
k k kT k k
kn λ + += =+
bb
M q Vq A F (49)
where node k locates at the center of the rigid body, {}
k
q is the nodal coordinate at node
k, []
b
M is a mass matrix, [ ]
b
V is a gyroscopic matrix of the rigid body, and {}
k
F is an
external force vector on node k. []
k
A in the constraint term is obtained by
[]
{}
k
k
k
Γ ∂
=
∂
A
q
(50)
43
and
k
λ is the Lagrange multiplier for k th node. Since gravity is not considered here and
there is no elastic constraint between the rigid body and the attached ends of flexible
bodies, there is no elastic term for this EOM.
For the next step, the EOM of the whole system is accomplished by combining EOM
of each flexible element and the rigid body using a method of superposition. The EOM of
the whole system is expressed as
1
2
[ ]{ } [ ]{ } [ ]{ } [ ] { } { },
{}
{}
{ }
{}
T
n
++ + =
=
M q V q K q A ΛF
q
q
q
q
(51)
If we write []
i e
M with sub-matrices as follows,
( ) ()
1
( ) ( )
[ ][ ]
{}
[ ]{ }
[ ][ ]
{}
i
i aa i ab
ii
i
i ba i bb
+
=
ee
ee
ee
M M
q
Mq
MM
q
(52)
by super position, the mass matrix of the whole system yields
44
1( ) 1( )
88
1( ) 1( ) 2( ) 2( )
2( ) 2( ) 3( )
88
1( ) ( )
( )
[ ] []
[ ] [ ][ ] [ ]
[ ] [ ][ ]
[]
[ ][ ][ ]
[]
aa ab
ba bb aa ab
ba bb aa
k bb k aa
n bb
×
×
−
+
+
=
++
e e
e ee e
e e e
e eb
e
MM 0
M MM M
0 M MM
M
M MM
M
(53)
Note that []
b
M is added to the sum of
1( )
[]
k bb − e
M and
( )
[]
k aa e
M . [] V and [ ] K can be
obtained in the similar manner. The constraint for the whole system is
1
2
{}
n
Γ
Γ
Γ
=
Γ
(54)
by differentiating Eq. (54) with respect to { } q , we get matrix [] A as follows,
{}
[]
{ }
∂
=
∂
Γ
A
q
(55)
and {} Λ is the Lagrange multiplier vector for the whole system.
1
2
{}
n
λ
λ
λ
=
Λ
(56)
45
Finally, the EOM of the system in Eq. (51) are given in an augmented form in order
to make a non-singular matrix [28]
11
[ ]{ } { } [ ] []
[ ]{ } { }
[]
{} {}
T
c
nn nn ×× ×
++ =
Kq F MA
Vq q
00 A0
V Λ
(57)
with
1 21 21 21 2
01 2 3
22 22 22 22
01 2 3
2 2 2 2
01 2 3
( )( )( )( )
() () () ()
{}
() () () ()
c
n n n n
p pp p
p pp p
pppp
++ +
++ +
=
+++
V
(58)
We will derive EOM from the element of the flexible body and EOM from the rigid body
in the following sections.
3.3 EOM of the Element of the Flexible Body
In this section, Zhao and Ren’s Euler-Bernoulli beam model [25] will be used for the
flexible body. In this paper, however, differentiations are explicitly solved rather than
using numerical differentiation, which has relatively less accuracy compared with
numerical integration or the Runge-Kutta method that are used here. Thus, in this section,
we will get [ ],[ ] and [ ]
ii i ee e
MV K in Eq. (43) from Zhao and Ren [25] and solve the
differentiation.
46
Fig. 18 Two-nodal beam element
3.3.1 Nodal Coordinates and Interpolation of the Position
Let’s consider the element between node I and II as shown in Fig. 18. The
coordinates of the element are expressed with nodal coordinates I and II as follows:
e
{}
{} ,
{}
{} { }
{} {} and { } { }
I
II
I II
I I II II
I II
εε
=
= =
q
q
q
rr
qp q p
(59)
The cross-section of the beam is denoted as a grey square in Fig. 18. We assume that the
cross-section is rigid and its shape and size are consistent throughout the element; its
geometric center, which lies on the centerline, and center of gravity are identical; and its
47
material reference frame is represented by orthonormal vectors
[ ] { },{ },{ } l mn . The vector
{ } l is normal to the cross-section and parallel with the unit tangent {} t of centerline. The
material reference frame of the cross-section is expressed in terms of { } p as follows:
2 222
01 2 3
12 0 3
13 0 2
12 0 3
2 222
0 21 3
2 3 01
13 0 2
2 3 01
2 222
0 31 2
{ } 2( ) ,
2( )
2( )
{} ,
2( )
2( )
{ } 2( )
p ppp
pp p p
pp p p
pp p p
p ppp
p p pp
pp p p
p p pp
p ppp
+ −−
= +
−
−
= + −−
+
+
= −
+ −−
l
m
n
(60)
The vectors { },{ },{ }
I II
l mn and { },{ },{ }
II II II
l mn are material reference frames of
the cross-sections at node I and II, respectively. The EOM of the element is accomplished
by integrating the EOM of the cross-section in Fig. 18 from node I to node II, and the
EOM of the cross-section is accomplished by the Lagrange equation with the position
{} r and the orientation { } p of the cross-section. The position vector {} r is given by [25]
1 2 3 4
{ ( )} ( ){ } ( ) (1 ){ } ( ){ } ( ) (1 ){ }
I I I II II II
N NL N NL ξ ξ ξε ξ ξε = ++ + ++ r r lr l (61)
with shape function [35],
23 2
12
22
34
() 1 3 2 () (1 )
() (3 2 ) () ( 1)
NN
NN
ξ ξ ξ ξξ ξ
ξξ ξ ξξ ξ
=−+ =−
= −=−
(62)
Here, xL ξ = , and x is the arc length of centerline at initial stress free configuration and
L is the total length of element. Orientation { } p is expressed in quaternion form as
follows from [25].
48
( )
01 2 3
() , , , () ()
I
p p p p ξ ξξ = =⊗⊗ pp π b
(63)
with
( )
( )
( ) ( )
( )
01 2 3
01 2 3
01 2 3
,,,
() , , ,
1 { ( )} { } { ( )} { } { ( )} { }
, 0, ,
2
2 1 { ( )} { } 2 1 { ( )} { }
() , , ,
cos , sin ,0,0
2 2
I IIII
TI T I T I
TI TI
pppp
b bb b
ξ π ππ π
ξξ ξ
ξ ξ
ξ
αξ αξ
=
=
+−
=
++
=
=
p
π
t l tn tm
tl tl
b
(64)
Here, characters with tildes mean quaternions, ⊗ is quaternion multiple, { ( )} ξ t is unit
tangent, and α is defined as follows [25]:
1
22
01
2arcsin
µ
α
µµ
=
+
(65)
with
0 00 11 22 33
1 01 1 0 2 3 3 2
I II I II I II I II
I II I II I II I II
pp pp pp pp
pp pp pp pp
µ
µ
= +++
= −−+
(66)
3.3.2 Mass Matrix
Mass matrix is obtained by kinetic energy in an Lagrange equation,
1
{} [ ]{}
2
T
T =
e ee
q Mq . By Zhao and Ren [25], the kinetic energy of the element is given
by
1
0
1
{ } { } 4{ } [ ] [ ][ ]{ }
2
T T T
T LA d ρ ξ = +
∫
r r p E JE p (67)
49
Here, ρ, L and A are density, length, area of the cross-section of the element, respectively,
[] E is the rotation matrix [29] defined as
10 3 2
2 30 1
3 2 10
[]
pp p p
p pp p
p p pp
−−
=−−
− −
E (68)
and [] J is an inertia matrix given by [] [ , , ]
x y z
diag J J J = J , where , and
xy z
JJ J are
moments of inertia of the cross-section with respect to x, y and z axes respectively. By
using a Lagrange equation and the following relationship,
{} { }
{ } {}, { } {}
{} {}
∂∂
= =
∂∂
ee
ee
rp
r qp q
qq
(69)
the mass matrix []
e
M has been derived from Eq. (67).
1
0
{} {} { } { }
[ ] 4 [] [ ][]
{} {} {} {}
TT
T
md ξ
∂∂ ∂ ∂
= +
∂∂ ∂ ∂
∫
e
ee e e
rr p p
M E JE
qq q q
(70)
Here, m is the mass of the element, and [ ] [ ]/ A = JJ .
Now we have to find { }{} and { }{} ∂∂ ∂ ∂
e e
r q pq in Eq. (70). By differentiating Eq.
(61) with respect to
e
{} q , we have
12 2 3 4 4
{}
[ ] (1 )[ ] { } [ ] (1 )[ ] { }
{}
I I I II II II
N NL N L N N L N L εε
∂
=+ +
∂
ll
e
r
I Δ lI Δ l
q
(71)
where [] I is identity matrix, and []
i
l
Δ is defined as
01 2 3
3 21 0
23 0 1
{}
[] 2
{}
ii i i
i
i i i i i
i
i i i i
pp p p
pp p p
pp p p
−−
∂
= =
∂
−−
l
l
Δ
p
(72)
To get { }{} ∂∂
e
pq , Eq. (63) is re-written in the following matrix form:
50
0 0 1 2 30 01 2 3
1 1 0 3 21 1 0 3 2
2 2 3 0 12 23 0 1
3 3 2 1 03 3 21 0
{ } [ ][ ]{ }
I
IIII
II II
II I I
II II
pb pppp
pb pp pp
pb pp p p
pb pp pp
π π π π
ππ π π
π π π π
π ππ π
=
− − − −−−
− −
=
− −
− −
pP Π b
(73)
By differentiating both sides of Eq. (73), we can get { }{} ∂∂
e
pq . However, since
differentiating the right side of Eq. (73) with product rule results in tensors of order 3
(TO3), we use indices [36] for differentiation. The elements of { }{} ∂∂
e
pq is obtained by
the following index notation:
( )
, , , ,
,
p P b P bP bP b ,
, , 1,2,..., 4 and 1, 2,...,16
I II I
i l ij jk k ij l jk k ij jk l k ij jk k l
l
i j k l
= Π = Π + Π + Π
= =
(74)
The matrix form of this equation is given in Appendix C1. Eq. (74) follows ‘summation
convention’ in [36], and commas in the equation represent differentiation with respect to
{}
e
q . Then, to solve Eq. (74), we have to find
,, ,
P , and b
I
ij l jk l k l
Π , which are expressed as
[ ] {}, [ ] {} and { }{}
I
∂ ∂ ∂ ∂ ∂∂
ee e
Pq Π q b q respectively in matrix form. [ ] {}
I
∂∂
e
Pq is
determined by differentiating []
I
P in Eq. (73). However, since [ ] {}
I
∂∂
e
Pq is TO3, we
differentiate each row of []
I
P separately to write it in TO2, and substitute them into Eq.
(74). Then [ ] {}
I
∂∂
e
Pq is expressed as
43 49
[]
[ ] , 1,2,3,4
{}
I
i
i
i
××
∂
= =
∂
r
e
P
0G 0
q
(75)
with
51
12
34
1 000 0 1 0 0
0 1 0 0 1 00 0
[] , [ ] ,
0 0 1 0 000 1
0 0 0 1 0 01 0
0 01 0 0 0 01
000 1 0 0 1 0
[] , [ ]
1 00 0 0 1 00
0 1 0 0 1 000
−
= =
− −
−
−
= =
−
GG
GG
(76)
Here, []
I
i r
P denotes i th rows of matrices []
I
P , thus [ ] {}
I
i
∂∂
re
Pq is 4 by 16 matrices.
[ ] {}
i
∂∂
re
Πq and { }{} ∂∂
e
bq need derivatives of π ’s and b’s. The derivatives of π ’s
are
0
2
2
32 22
3
32 2
1 1
{} {} {}
81
1 11
{} {} {} {} {}
21 21
1
{} {}
21
ll
l
n l l nn
l l
m ll
l
π ββ γ
γγ
β
γ
βγ β β β β π γγ
γγ γ γ
β β
γ γ
π βγ β β
γγ
β
γ
∂∂ ∂
= −
∂ ∂∂
+
∂∂ ∂ ∂∂
= − − −
∂ ∂∂ ∂ ∂
+ +
∂∂
= −−
∂∂
+
e ee
e ee e e
ee
q qq
q qq q q
qq
2
11
{} {} {}
21
mm
l
ββ γγ
γγ
β
γ
∂ ∂∂
+ −
∂ ∂∂
+
e e e
q qq
(77)
with following variables:
{ } { }, { } { },
{ } { },
TI T I
lm
TI
n
ββ
βγ
′′ = =
′′ = =
r l rm
rn r
(78)
where
12 3 4
{ } { } (1 ){ } { } (1 ){ }
I I I II II II
N NL N NL εε
′′ ′ ′
′= ++ + ++ rr l r l (79)
52
In Eq. (79), apostrophes represent the derivative with respect to ξ . See Appendix B1 for
{}, {}, {} and {}
l m n
β β β γ ∂ ∂ ∂ ∂ ∂ ∂ ∂∂
e ee e
q qq q in Eq. (77). Then, [ ] {} ∂∂
e
Π q is
expressed as follows:
0
1 16
2
3
{}
[]
[ ] , 1,2,3,4
{}
{}
{}
i
i
i
π
π
π
×
∂∂
∂
= =
∂
∂∂
∂ ∂
e
r
e
e
e
q
0
Π
G
q
q
q
(80)
The derivatives of b’s are expressed as
0 1
3 2
1 16 1 16
sin , cos ,
{} 2 2 {} {} 2 2 {}
,
{} {}
b b
b b
αξ ξ α αξ ξ α
××
∂ ∂ ∂∂
= − ⋅⋅ = ⋅⋅
∂ ∂∂ ∂
∂ ∂
= =
∂∂
e ee e
ee
q qq q
00
qq
(81)
with
0 1
01 22
01 e e
2
{ } {} {}
µ µ α
µµ
µµ
∂ ∂ ∂
= −
∂ +∂ ∂
e
q qq
(82)
where
0
0 1 2 3 01 2 3
1 3 1 4
1
1 0 3 2 1 03 2
1 3 1 4
0
{}
0
{}
IIIIIIII IIII
II II II II I I I I
p p p p pppp
p p p p ppp p
µ
µ
××
××
∂
=
∂
∂
= − − − −
∂
e
e
00
q
00
q
(83)
Then, { }{} ∂∂
e
bq is given by
0
1
1 16
1 16
{}
{}
{ }
{}
b
b
×
×
∂∂
∂∂
∂
=
∂
e
e
e
q
q
b
0
q
0
(84)
53
3.3.3 Gyroscopic Matrix
We will find the gyroscopic term []
e
V in Eq. (43). The gyroscopic term is given by
[25] as follows:
( ) ( ) [ ]{ } [ ]{ }
1
[]
{} 2 {}
T
∂∂
= −
∂∂
ee
e
ee
M q M q
V
qq
(85)
From Eq. (70) and (85), we have the following expression for gyroscopic term:
( ) ( )
1
0
1
[ ] [ ] 4[ ] [ ] 4[ ]
2
T
md ξ
= + − +
∫
e 12 12
V VV VV (86)
with
ee e e
{} {} { } { }
{ } [] [ ][] { }
{} {} {} {}
[] , [ ]
{} {}
TT
T
∂∂ ∂ ∂
∂∂
∂∂ ∂ ∂
= =
∂∂
ee
12
ee
rr p p
q E JE q
qq q q
VV
qq
(87)
2
[ ] and [ ]
1
VV can be solved by using indices. The elements of []
1
V with indices yield
( )
1 ,, e , , e ,, e
,
V rr q r r q r r q ,
, , 1,2,...,16 and 1,2,3
il ji j k k jil j k k ji j kl k
l
ik l j
= = +
= =
(88)
See Appendix C2 for the matrix form of this equation. In order to solve this equation, the
second derivative of {} r ,
,
r
j il
(or
,
r
j kl
) has to be determined, and it can be obtained by
differentiating Eq. (71) with respect to {}
e
q . Since the second derivative of {} r is TO3,
we find
2
ee
{} {}
i
r ∂
∂∂ qq
as follows so that we can substitute them into Eq. (88).
54
( )
( )
2
ee
3 3 3 4 31 3 3 3 4 31
22
4 3 4 4 4 1
1 3 1 4
3 3 3 4 31
44
{} {}
2 (1 )[ ] [ ]
00
2 (1 )[ ] [ ]
.0
1,2,3
i
T
II
ii
T
II II
ii
r
N L N L
N L N L
symm
i
ε
ε
× × ×× × ×
×× ×
××
×× ×
∂
=
∂∂
+
+
=
lr
lr
qq
0 0 00 0 0
H Δ0 0 0
00
00 0
H Δ
(89)
with
1
2
3
10 0 0
01 0 0
[] ,
00 1 0
00 0 1
000 1
0 01 0
[] ,
01 0 0
1 000
0 0 10
000 1
[]
1 000
01 0 0
=
−
−
=
−
=
−
H
H
H
(90)
In the same manner, the element of [ ]
2
V is acquired by using indices as follows:
55
( )
2 , ,e
,
, ,e , , ,e
, , ,e , , e
V p E JE p q
p E JE p q p E JE p q
p E JE p q p E JE p q ,
, , 1,2,...,16, , 1,2,...,4, and , 1,2,3
io j i kj kl lm m n n
o
j io kj kl lm m n n j i kj o kl lm m n n
j i kj kl lm o m n n j i kj kl lm m no n
i no j m k l
=
= +
++
= = =
(91)
The matrix form of this equation is given in Appendix C3. To solve this equation,
e
[]
{}
∂
∂
E
q
(
,
E
kj o
or
,
E
lm o
in Eq. (91)) and
2
ee
{ }
{} {}
∂
∂∂
p
qq
(
,
p
j io
or
,
p
m no
in Eq. (91)) have to be
determined. By differentiating Eq. (68) with respect to {}
e
q ,
e
[]
{}
∂
∂
E
q
is expressed as
e e
[] { }
[ ] , 1,2,3
{} {}
i
i
i
∂ ∂
= =
∂∂
r
E p
H
qq
(92)
with
1
2
3
0 10 0
1 0 00
[] ,
0 0 01
0 0 10
00 1 0
00 0 1
[] ,
10 0 0
01 0 0
000 1
0 01 0
[]
0 10 0
1 000
−
=
−
−
−
=
−
=
−
H
H
H
(93)
2
ee
{ }
{} {}
∂
∂∂
p
qq
is obtained by differentiating Eq. (74) with respect to {}
e
q as follows with
index notation:
56
( )
,
,
, ,, , , ,
,, , , , , ,
p Pb
P bP bP bP b
P b P bP bP b ,
, and 1,2,...,4 and and 1,2,...,16
I
i lm ij jk k
lm
I III
ij lm jk k ij l jk m k ij m jk l k ij jk lm k
II II
ij jk l k m ij m jk k l ij jk m k l ij jk k lm
i j k l m
= Π
= Π + Π + Π + Π
+ Π + Π + Π + Π
= =
(94)
The matrix form of this equation is given in Appendix C4. In this equation,
2
ee
[]
{} {}
I
∂
∂∂
P
qq
(
,
P
I
ij lm
),
2
ee
[]
{} {}
∂
∂∂
Π
qq
(
, jk lm
Π ) and
2
ee
{ }
{} {}
∂
∂∂
b
qq
(
,
b
k lm
) have to be determined. The second
derivatives of []
I
P is
2
ee
[]
{} {}
I
∂
=
∂∂
P
0
qq
, and
2
ee
[]
{} {}
∂
∂∂
Π
qq
and
2
ee
{ }
{} {}
∂
∂∂
b
qq
are obtained by
differentiating Eq. (80) and Eq. (84) respectively as follows:
2
0
ee
1 16
2
2
2
e
ee
2
3
ee
{} {}
[]
[ ] , 1,2,3,4 and 1,2,...,16
{}
{} {}
{} {}
j
i
i
j
j
j
ij
q
π
π
π
×
∂
∂∂
∂
= = =
∂
∂∂
∂∂
∂
∂∂
r
r
r
r
qq
0
Π
G
q
qq
qq
(95)
For matrices
22 2
03 2
ee ee ee
,,
{} {} {} {} {} {}
ππ π ∂∂ ∂
∂∂ ∂∂ ∂∂ qq qq qq
, see Appendix B2. Since
2
ee
{ }
{} {}
∂
∂∂
b
qq
is OT3, we find TO2,
2
ee
{} {}
i
b ∂
∂∂ qq
so that we can substitute into Eq. (94).
57
2 22
0
2 22
1
2
2
16 16
2
3
16 16
sin cos
{} {} 2 2 {} {} 2 4 {} {}
cos sin
{} {} 2 2 {} {} 2 4 {} {}
{} {}
{} {}
T
T
b
b
b
b
αξ ξ α αξ ξ α α
αξ ξ α αξ ξ α α
×
×
∂ ∂ ∂∂
= − ⋅⋅ − ⋅
∂∂ ∂∂ ∂ ∂
∂ ∂ ∂∂
= ⋅⋅ − ⋅
∂∂ ∂∂ ∂ ∂
∂
=
∂∂
∂
=
∂∂
ee ee e e
ee ee e e
ee
ee
qq qq q q
qq qq q q
0
qq
0
qq
(96)
In this equation,
{}
α ∂
∂
e
q
is given in Eq. (82), and
2
{} {}
α ∂
∂∂
ee
qq
yields
( )
2
2 2
0 0 0 11 1
01 22
01
00 11
01 0 1 2
22
01
{} {}
2
{} {} {} {} {} {} {} {}
4
{} {} {} {}
TT
T
α
µµ µ µµ µ
µµ
µµ
µµ µµ
µµ µ µ
µµ
∂
∂∂
∂∂ ∂ ∂∂ ∂
= − +−
+ ∂ ∂ ∂ ∂ ∂∂ ∂∂
∂∂ ∂∂
− −+
∂∂ ∂ ∂
+
ee
e e e e ee ee
ee e e
qq
q q q q qq qq
qq q q
(97)
where
58
3 3 3 4 31 3 3 3 4 31
4 4 4 1 4 3 4 1
2
0 1 3 1 4
3 4 31
4 1
3 3 3 4 31 3 3 3 4 31
4 4 4 1 4 3 4 1
2
1 1 3 1 4
ee
3 4 31
4 1
[]
00
{} {}
.0
[]
00
{} {}
.0
symm
symm
µ
µ
× × ×× × ×
× ×× ×
××
××
×
× × ×× × ×
× ×× ×
××
××
×
∂
=
∂∂
∂
=
∂∂
ee
μ
0 0 00 0 0
0 00 I 0
00
qq
0 0
0
0 0 00 0 0
0 00 T 0
00
qq
00
0
(98)
with
0 10 0
00 0
[]
01
0 symm
=
−
μ
T (99)
3.3.4 Stiffness Matrix
We will find the stiffness term []
e
K in Eq. (43), now. The stiffness term is given by
[25] as follows:
1
3 12
2 3
0
[ ]{ }
{} {} {} {}
xy z
L EA GJ EJ EJ d
κ κκ ε
ε κ κ ξ
∂ ∂∂ ∂
= + + +
∂∂ ∂ ∂
∫
ee
ee e e
Kq
qq q q
(100)
59
In this equation, we get [ ]{ }
ee
Kq instead of []
e
K . and EG are modulus of elasticity and
shear modulus of elasticity respectively, and
1 2 3
, , and ε κ κ κ are normal strain, torsion
and two bending curvatures expressed as follows [25]:
1 23 2 22
1 1 1
1, { } { } , { } { }, { } { }
T TT
L L LL L
γα
εκ κ κ ′′ ′′ ′′ = −= + = − = rc rn rm (101)
where
12 3 4
{ } { } (1 ){ } { } (1 ){ }
I I I II II II
N NL N NL εε
′′ ′′ ′′ ′′
′′= ++ + ++ rr l r l (102)
and
[ ] {}
{ }
1 { } {}
IX
IT
=
+
lt
c
lt
(103)
In Eq. (103), []
IX
l is the skew-symmetric matrix of
{ }
1 23
{}
T
I III
lll = l defined as
32
31
21
0
[] 0
0
II
IX I I
II
ll
l l
ll
−
= −
−
l (104)
Now the derivatives
3 12
, , and
{} {} {} {}
κ κκ ε ∂ ∂∂ ∂
∂∂∂ ∂
e e e e
qqq q
in Eq. (100) need to be determined.
By differentiating Eq. (101) with respect to {}
e
q , we have the following expressions.
1
2
2
2
3
2
1
{} {}
1 { } { } 1
{ } { }
{} {} {} {}
1 { } { }
{ } { }
{} {} {}
1 { } { }
{ } {}
{} {} {}
TT
TT
TT
L
LL
L
L
ε γ
κ α
κ
κ
∂∂
= ⋅
∂∂
′′ ∂ ∂ ∂∂
′′ = + + ⋅
∂ ∂ ∂ ∂
′′ ∂ ∂∂
′′ = −+
∂ ∂∂
′′ ∂ ∂∂
′′ = +
∂ ∂∂
ee
e e ee
e ee
e ee
qq
r c
cr
q q q q
rn
nr
q qq
r m
mr
q qq
(105)
60
Here,
{}
γ ∂
∂
e
q
and
{}
α ∂
∂
e
q
are given in Appendix B1 and Eq. (82) respectively, and
{ }
{}
′′ ∂
∂
e
r
q
is
12 2 3 4 4
{ }
{}
[ ] (1 )[ ] { } [ ] (1 )[ ] { }
I I I II II II
N NL NL N NL NL εε
′′ ∂
=
∂
′′ ′′ ′′ ′′ ′′ ′′
++
e
l l
r
q
I Δ lI Δ l
(106)
and
{ }
{}
∂
∂
e
c
q
is given by
( )
( )
2
[ ]{ }
{ } 1 [ ] { }
{} {} {} {}
XI
XI
l
l
l
β γ
γβ
γβ
′ ∂
′ ∂ ∂∂
= − − +
∂ +∂ ∂ ∂
+
e e ee
rl
c rl
q q qq
(107)
where
( )
1 34 4
31
[ ]{ }
{}
[] [ ] [ ] [] (1 )[] [ ] [ ] { }
XI
IX X I IX II IX II IIX I
N N NL NL ε
×
′ ∂
=
∂
′ ′′ ′
′ −+
e
ll
rl
q
lr Δ 0 l lΔ l l
(108)
{ }
{}
∂
∂
e
n
q
and
{}
{}
∂
∂
e
m
q
in Eq. (105) are found as follows:
{ } {} { } { } { } { }
,
{} { } {} {} { } {}
∂ ∂∂ ∂ ∂ ∂
= =
∂ ∂∂ ∂ ∂ ∂
e e e e
n np m m p
q pq q p q
(109)
with
2 3 0 1 32 1 0
1 032 0 1 2 3
0 1 23 1 0 3 2
{ } { }
2 ,2
{ } { }
p p p p p p p p
p p p p p pp p
p p pp p p p p
−−
∂∂
=− − = − −
∂∂
−−
nm
pp
(110)
61
3.4 EOM of the Rigid Body
The EOM of the rigid body is given by Eq. (49), in which, []
b
M and [ ]
b
V will be
derived in this section in the same way as Section 3.3. However, the generalized
coordinates of the rigid body have 8 components since the rigid body has only one node,
and there is no shape function for the rigid body. The generalized coordinates of the rigid
body {}
k
q are defined as
{ }
{} {} , ( 1) / 2
k
kk
k
kn
ε
= = +
r
qp (111)
where k and n are the node numbers on the center and on the right end on Fig. 17
respectively. As we did in Section 3.3, we can get the mass matrix from the kinetic
energy of the rigid body, which is similar to Eq. (67) as follows:
1
{ } { } 2{ } [ ] [ ][ ]{ }
2
k T k k T k T k k
b
Tm = +
b
r r p E JE p (112)
where []
k
E is given by Eq. (68) with {}
k
p , []
b
J is mass moment of inertia of the rigid
body with respect to material coordinate at node k, and
b
m is mass of the rigid body.
Then, the mass matrix []
b
M is obtained as
{ } { } { } { }
[ ] 4 [] [ ][]
{} {} {} {}
TT
k k k k
kT k
b k k k k
m
∂∂ ∂ ∂
= +
∂∂ ∂ ∂
bb
r r p p
M E JE
qq q q
(113)
where
3 8 4 3 4 1
{ } { }
[] , []
{} {}
kk
kk
× × ×
∂∂
= =
∂∂
rp
I 0 0I 0
qq
(114)
62
The gyroscopic term is obtained from the following equation [25].
( ) ( )
[ ]{ } [ ]{ }
1
[ ]
{} 2 {}
T
k k
k k
∂∂
= −
∂∂
bb
b
Mq Mq
V
qq
(115)
63
Chapter 4
Dynamic Analysis of a Rigid-Flexible
Coupled System
The satellite model in Chapter 3 is simulated in this section. Even though a
quaternion-based model has the merit of avoiding singularity, it is not appropriate to
show the physical orientation of bodies. Therefore, in this simulation, Euler angles are
partially used for the presentation of orientations. The initial conditions are given with
Euler angles, and they are converted to quaternions for simulation. After the simulation,
the results with quaternions are converted again to Euler angles to show the orientations.
Note that the conversion between Euler angles and quaternions makes no singularity
issue. It has been proved in [26]. The Euler angle sequence in this simulation is
( ) () ( ) x yz ψ φθ →→ . The physical properties used in this simulation are shown in Table
1 and 2, and two beam elements are used for each arm.
64
Property Symbol Value
Length L 1m
Cross-Section A
32
2.5 10 m
−
×
Density ρ
33
2.8 10 / kg m ×
Moments of inertia of cross-
section
x
J
6 4
1.04 10 m
−
×
y
J
74
5.2 10 m
−
×
z
J
74
5.2 10 m
−
×
Modulus of elasticity E
10 2
7.17 10 / N m ×
Shear modulus of elasticity G
10 2
2.68 10 / N m ×
Table 1 Physical parameters of beam elements
Property Symbol Value
Edge h 0.2m
Density ρ
33
2.8 10 / kg m ×
Mass moments of inertia
[]
b
J
2
0.149[ ]kg m ⋅ I
Table 2 Physical parameters of the rigid body
65
Two cases of vibrations are considered in the simulation:
Case 1. The satellite is in free vibration that is caused by initial rotational velocities
(0) ψ= 0.1 rad/s, (0) φ=
0.5 rad/s and (0) θ=
0.1 rad/s. All other displacement and
velocity parameters are initially zeros.
Case 2. The satellite is in forced vibration that is excited by an external torque given by
0
0
sin 2 1 , 0
() 22
0,
with 15 and 1
t
tT
t T
Tt
Nm T s
τ π
π
τ
τ
− + ≤<
=
≤
= =
(116)
The profile of the torque in Eq. (116) is plotted in Fig. 19. The satellite has initial
rotational velocities of (0) ψ= 0.1 rad/s, (0) φ=
0.1 rad/s and (0) θ=
0.1 rad/s, and all
other displacement and velocity parameters are initially zeros.
Fig. 19 Profile of the external torque
The dynamic response of the satellite is plotted in Figs. 20-31. In the legends,
subscripts represent numbers of nodes specified in Fig. 17. Each beam has 2 elements, so
66
node 1 designates left end of the satellite, node 3 for the center of the rigid body, and
node 5 for the right end of the satellite. u, v and w in Figs. 20 and 26 are the
displacements in the direction of x, y and z specified in Fig. 17. Displacements in Figs. 20
and 26 are the displacements from their initial point, not from the origin. The orientations
of the satellite in Euler angle in Figs. 24, 25, 30 and 31 are respectively converted from
Figs. 22, 23, 28 and 29 that are the quaternion representations of orientations.
4.1 Free Vibration (Case 1)
The dynamic response of the satellite in the free vibration (Case 1) is plotted in Figs.
20-25. Fig. 24 shows that Euler angle φ reaches 90 ° at 2.5 sec and 270° at 6.7 sec. When
Euler angle φ reaches 90 ° or 270°, we can observe a jump in the profile of θ. This is not
a singular point. This jump is caused by the definition of the Euler angle. Fig. 22 shows
smooth curves. Since Fig. 24 is converted from Fig. 22, there cannot be singularity.
However, if Euler angles were used in the simulation, singularity would occur when φ
reached 90° or 270°. In Fig. 24, we can see that the angular displacements of all of the
nodes are the same. This means there is no bending on the beams, which is a reasonable
result since no external force is applied in this simulation.
67
Fig. 20 Free vibration (Case 1): translational displacements of the body
68
Fig. 21 Free vibration (Case 1): translational velocities of the body
69
Fig. 22 Free vibration (Case 1): rotational displacements of the body in terms of
quaternion components
70
Fig. 23 Free vibration (Case 1): rotational velocities of the body in terms of quaternion
components
71
Fig. 24 Free vibration (Case 1): rotational displacements of the body in terms of Euler
angles
72
Fig. 25 Free vibration (Case 1): rotational velocities of the body in terms of Euler angles
73
4.2 Forced Vibration (Case 2)
In the forced vibration simulation (Case 2), the torque prescribed in Eq. (116) is
applied on the center of the rigid body. The torque is applied about the material
coordinate { }
k
m on node k, so the rotational axis of the torque rotates together with the
rigid body of the satellite. The result is plotted in Figs. 26-31. Fig. 30 shows that the
Euler angle φ reaches 90 ° at 5 sec.
Like Case 1, no singularity is observed, and there is a jump in the profile when the
Euler angle φ reaches 90 °. However, there is bending in this case. The angular
displacements of the nodes in Fig. 30 are not identical because the beams are bent. Also,
Fig. 30 shows that angular displacements of nodes 1 and 5 are almost same. It is because
the satellite is bilateral symmetric; thus the angular displacements of both ends are the
same. Unlike Case 1, we can see vibration in Case 2. When torque is applied on the rigid
body, the system vibrates because of flexible beams attached on the satellite.
74
Fig. 26 Forced vibration (Case 2): translational displacements of the body
75
Fig. 27 Forced vibration (Case 2): translational velocities of the body
76
Fig. 28 Forced vibration (Case 2): rotational displacements of the body in terms of
quaternion components
77
Fig. 29 Forced vibration (Case 2): rotational velocities of the body in terms of
quaternion components
78
Fig. 30 Forced vibration (Case 2): rotational displacements of the body in terms of Euler
angles
79
Fig. 31 Forced vibration (Case 2): rotational velocities of the body in terms of Euler
angles
80
Chapter 5
Conclusion
In investigation of the motion of a rigid body subject to torques, four coordinate
systems are available for torque description: inertial frame, body frame, Euler basis, and
dual Euler basis. It is known that quaternions are non-singular, and that transformations
between Euler angles and quaternions are singularity-free. Nevertheless, the motion of a
quaternion-based model of the body subject to torques still can be unbounded, depending
upon which coordinate system is chosen for torque description. The singularity analysis
in all the coordinate systems yields conditions on torque-induced singularity, which are
summarized as follows:
(i) A torque described in an inertial or body frame or an Euler basis will never yield
unbounded response of a quaternion-based model of rigid bodies.
(ii) A torque applied in the direction of the first or the third dual Euler base vector can
cause unbounded response of a quaternion-based model of rigid bodies. Otherwise there
is no torque-induced singularity in the motion.
The torque in the above singularity conditions can be either an external one or an internal
one (like a spring/damping torque).
81
The dynamic model of the quaternion-based rigid-flexible coupled system has been
established in this study. Both free and forced vibration simulation is conducted for
comparison. The results of this modeling and simulation are summarized as follows:
(i) For the rigid-flexible coupled systems, FEM turns out to be a good solution. FEM
allows the EOM of the rigid body and flexible bodies to be coupled easily by super
position.
(ii) To avoid the singularity problem that usually occurs when Euler angles are used,
quaternion is used for the orientations of the system, and quaternion-based FEM is
derived. Quaternion-based FEM has not only the merit of singularity-free representations,
but also can be applied for large deformations.
(iii) Numerical differentiation has less accuracy compared with numerical integration
or ODE. Thus, differentiations in the governing equation are explicitly solved, which
renders the result solid.
(iv) Numerical simulations prove that the quaternion-based model has no singularity
in both free and forced vibrations unless the torque is applied on dual Euler angles.
Conversion between quaternions and Euler angles makes no singularity.
This singularity-free model that consists of rigid and flexible bodies can be usefully
applied on many applications. It can be applied on the systems that bend and rotate such
as flexible surgical robot arms and snake-arm robots, and also it can be applied on rotors
and radars.
82
References
[1] Huang, Y. A., Deng, Z., and Yao, L., 2007, "An Improved Symplectic Precise
Integration Method for Analysis of the Rotating Rigid-Flexible Coupled System,"
Journal of Sound and Vibration, 299(1-2), pp. 229-46.
[2] Bai, S., Huang, X., and Liu, Y., 2008, "Dynamic Modeling and Simulation of a
Flexible Satellite," 2008 Asia Simulation Conference - 7th International
Conference on System Simulation and Scientific Computing, ICSC 2008, pp.
1068-1072.
[3] Zhang, G., Lu, N., and Che, R., 2011, "Dynamic Analysis on Rigid-Flexible
Coupled Multi-Body System with a Few Flexible Components," ICQR2MSE
2011 - Proceedings of 2011 International Conference on Quality, Reliability, Risk,
Maintenance, and Safety Engineering, pp. 1010-1015.
[4] Hamilton, W. R., 1853, Lectures on Quaternions, Hodges and Smith, Dublin.
[5] Gibbs, J. W., 1901, Vector Analysis, Yale University Press, New Haven.
[6] Gray, A., 1918, A Treatise on Gyrostatics and Rotational Motion, Macmillan,
London.
[7] Pio, R. L., 1966, “Euler Angle Transformations,” IEEE Transactions on
Automatic Control, Oct., Vol. AC-11, No. 4, pp. 707-715.
83
[8] Stuelpnagel, J., 1964, “On the Parameterization of the Three-Dimensional
Rotation Group,” SIAM Rev., Vol. 6, No. 4, pp. 422-430.
[9] Mitchell, E. E. L. and Rogers, A. E., 1965, “Quaternion Parameters in Simulation
of Spinning Rigid Body,” Simulation, June, Vol. 4, No. 6, pp. 390-396.
[10] Bar-Itzhack, I.Y. and Oshman, Y., 1985, “Attitude Determination from Vector
Observations Quaternion Estimation,” IEEE Transactions on Aerospace and
Electronic Systems, Vol. AES-21, No. 1, pp. 128-36.
[11] Taylor, R.H., 1979, “Planning and Execution of Straight Line Manipulator
Trajectories,” IBM Journal of Research and Development, July, Vol. 23, No. 4,
pp. 424-36.
[12] Gu, Y. L., 1988, “Analysis of Orientation Representations by Lie Algebra in
Robotics,” Proceedings of the 1988 IEEE International Conference on Robotics
and Automation (Cat. No.88CH2555-1), Vol. 2, pp. 874-9.
[13] Kyung, M., Kim, M. and Hong, S., 1996, “A New Approach to Through-the-Lens
Camera Control,” Graphical Models and Image Processing, May, Vol. 58, No. 3,
pp. 262-85.
[14] Kavan, L., Collins, S., Zara, J. and O'Sullivan, C., 2008, “Geometric Skinning
with Approximate Dual Quaternion Blending,” ACM Transactions on Graphics,
Oct., Vol. 27, No. 4, pp. 105 (23 pp.).
[15] Ham, W., Im, K. and Khurelbaatar, Ts., 2009, “Computer Graphic Animation for
Helicopter Control,” 2009 IEEE International Symposium on Industrial
Electronics (ISIE 2009), pp. 1937-42.
84
[16] Hsiao, S. and Delosme, J., 1994, “Parallel Processing of Complex Data Using
Quaternion and Pseudo-Quaternion CORDIC Algorithms,” Proceedings of the
International Conference on Application Specific Array Processors (Cat.
No.94TH0687-4), pp. 324-35.
[17] Yamamoto, H. and Aoshima, N., 2002, “Three Dimensional Signal Processing
Based on Quaternion,” SICE 2002. Proceedings of the 41st SICE Annual
Conference (Cat. No.02TH8648), Vol. 3, pp. 1420-4.
[18] Markey, B. R. and Walmsley, S. H., 1982, “Quaternion formulation of lattice
dynamics of molecular crystals,” Chemical Physics Letters, Oct. Vol. 92, No. 3,
pp. 257-61.
[19] Rapaport, D.C., 1985, “Molecular dynamics simulation using quaternions,”
Journal of Computational Physics, Sept., Vol. 60, No 2, pp. 306-14.
[20] Miller III, T. F., Eleftheriou, M., Pattnaik, P., Ndirango, A., Newns, D. and
Martyna, G. J., 2002, “Symplectic quaternion scheme for biophysical molecular
dynamics,” Journal of Chemical Physics, May, Vol. 116, No. 20, pp. 8649-8659.
[21] Hrennikoff, A., 1941, “Solution of problems of elasticity by framework method,”
Journal of Applied Mechanics, 8(4), pp. 169-175.
[22] Shabana, A., and Wehage, R. A., 1983, “Variable Degree-of-Freedom Component
Mode Analysis of Inertia Variant Flexible Mechanical Systems,” ASME J. Mech.,
Transm., Autom. Des., 105, pp. 371–378.
[23] Wallrapp, O., and Schwertassek, R., 1991, “Representation of Geometric
Stiffeningin Multibody System Simulation,” Int. J. Numer. Methods Eng., 32, pp.
1833–1850.
85
[24] Meijaard, J. P., 1996, "Validation of Flexible Beam Elements in Dynamics
Programs," Nonlinear Dynamics, 9(1-2), pp. 21-36.
[25] Zhao, Z., and Ren, G., 2012, "A Quaternion-Based Formulation of Euler-
Bernoulli Beam without Singularity," Nonlinear Dynamics, 67(3), pp. 1825-1835.
[26] Liu, D., Xia, Q., and Wen, Q., 2010, "Attitude solution of vertical penetrative
trajectory based on quaternion," CMCE 2010, Vol. 1, pp. 293-296.
[27] Hairer, E., and Wanner, G., 1996, Solving Ordinary Differential Equations II:
Stiff and Differential-Algebraic Problems, 2nd ed., Springer, Berlin.
[28] Nikravesh, P. E., 1998, Computer-Aided Analysis of Mechanical Systems,
Prentice-Hall, New Jersey
[29] Shabana, A. A., 2005, Dynamics of Multibody Systems, Cambridge University
Press, New York.
[30] O’Reilly, O. M., 2007, “The dual Euler basis: Constraints, potentials, and
Lagrange’s equations in rigid body dynamics,” ASME Journal of Applied
Mechanics, 74(2), pp. 256–258.
[31] Senan, N. A. F., and O’Reilly, O. M., 2009, “On the Use of Quaternions and
Euler Rodrigues Symmetric Parameters with Moments and Moment Potentials,”
International Journal of Engineering Science, 47(4), pp. 595–609.
[32] Kuipers, J. B., 1999, Quaternions and Rotation Sequences, Princeton University
Press, New Jersey.
86
[33] Shuster, D. M., 1993, "A Survey of Attitude Representations," The Journal of the
Astronautical Sciences, 41(4), pp. 439-517.
[34] Howland, R. A., 2006, Intermediate Dynamics: A linear Algebraic Approach,
Springer Science+Business Media, Inc., New York.
[35] Thomson, W. T., and Dahleh, M. D., 1997, Theory of Vibration with Applications,
Prentice-Hall, New Jersey
[36] Mal, K. A. and Singh, J. S., 1991, Deformation of Elastic Solids, Prentice Hall,
New Jersey.
87
Appendices
A1. The Proof of Equation (15)
Equation (14) can be written in the matrix form [32]
{ } [ ]{ }
NE
= r Pr (A1)
where vectors {}
N
r and {}
E
r contain the components of r
in frame {N} and frame {E},
respectively, and [] P is given in Eq. (16). On the other hand, by coordinate
transformation
{ } [ ]{ }
EN
p
= rR r
(A2)
where []
p
R is given in Eq. (16), and it is a direction cosine matrix describing the
sequence of rotations shown in Fig. 3. Because []
p
R is orthogonal, substituting Eq. (A1)
into Eq. (A2) yields
1
[ ] [] []
T
pp
−
= = PR R , which is the same as Eq. (15).
88
B1.
( )
( )
1 13 4 2
34 4
1 13 3
1 5
1 13
{ } { } { } (1 ){ } [ ]
{}
{} (1 ){ } [ ] { } {}
{} { } { } [ ] 0 {} [ ]
{}
{ } {}
{}
IT IT IIT II IIT I l
I T II I T II II T I
IT IT IIT I IT m
IT IT n
N N N NL NL
N NL NL
N NN N
N NN
β
ε
ε
β
β
×
∂
′ ′′ ′ ′
= + ++
∂
′′ ′
+
∂
′ ′′ ′
= +
∂
∂
′′
= +
∂
l
e
l
m
e
e
l rr l Δ
q
l l Δ ll
m rr Δ m0
q
nr
q
( ) 3
1 5
{ } [ ] 0 { } []
1 { }
{ }
{} {}
IIT I IT
T
N
γ
γ
×
′′
′ ∂∂
′ =
∂∂
n
ee
r Δ n0
r
r
qq
(B1)
where,
12 3 4
12 2 3 4 4
{ } { } (1 ){ } { } (1 ){ }
{ }
[ ] (1 )[ ] { } [ ] (1 )[ ] { }
{}
I I I II II II
I I I II II II
N NL N NL
N NL NL N NL NL
εε
εε
′′ ′ ′
′= ++ + ++
′ ∂
′′ ′ ′′ ′
=++
∂
ll
e
rr l r l
r
I Δ lI Δ l
q
(B2)
and
3 21 0
0 12 3
1 03 2
2 3 01
1 032
0 1 23
{}
[] 2 ,
{}
{}
[] 2
{}
ii i i
i
i i ii i
i
i iii
iii i
i
i i i i i
i
i i i i
pp p p
p pp p
p ppp
ppp p
p p pp
p p pp
− −
∂
= = − −
∂
∂
= =−−
∂
−−
m
n
m
Δ
p
n
Δ
p
(B3)
89
B2.
Set the following variable.
2
ee
2
e e
2
e e
1
{} {}
1
{} {}
1
{} {}
ll
mm
nn
ββ γ
γγ
β β γ
γγ
ββ γ
γγ
∂ ∂
= −
∂∂
∂ ∂
= −
∂∂
∂ ∂
= −
∂∂
l
m
n
F
qq
F
qq
F
qq
(B4)
where,
designate a row matrix. Then,
2
2
ee
3
11
{} { } { } {} {} {} {}
{ } { }
3 []
{} {} {} {}
TT
l l l
TT
l
ββ β γγ
γγ
β γγ
γ
∂ ∂ ∂ ∂ ∂∂
=− +
∂ ∂ ∂ ∂∂ ∂∂
′′ ∂∂ ∂∂
+ −−
∂∂ ∂∂
l
e e e e e
e e e e
F
q q q qq qq
rr
T
qq qq
(B5a)
2
2
ee
3
11
{} { } { } {} {} {} {}
{ } { }
3 []
{} {} {} {}
TT
mm m
TT
m
ββ β γγ
γγ
β γγ
γ
∂ ∂ ∂ ∂ ∂∂
=− +
∂ ∂ ∂ ∂∂ ∂∂
′′ ∂∂ ∂∂
+ −−
∂∂ ∂∂
m
e e e e e
e e e e
F
q q q qq qq
rr
T
qq qq
(B5b)
2
2
ee
3
11
{} { } { } {} {} {} {}
{ } { }
3 []
{} {} {} {}
TT
nn n
TT
n
ββ β γγ
γγ
β γγ
γ
∂ ∂ ∂ ∂ ∂∂
=− +
∂ ∂ ∂ ∂∂ ∂∂
′′ ∂∂ ∂∂
+ −−
∂∂ ∂∂
n
e e e e e
e e e e
F
q q q qq qq
rr
T
qq qq
(B5c)
with
2 22
3 12
12 3
ee ee ee
[]
{} {} {} {} {} {}
r rr
rr r
′ ′′
∂ ∂∂
′ ′′
= ++
∂∂ ∂∂ ∂∂
T
qq qq qq
(B6)
90
Here,
12 3
,, rr r
′ ′′
are components of { } ′ r , and
2
ee
{} {}
i
r
′
∂
∂∂ qq
can be obtained from Eq. (89)
by replacing
j
N
′
to
j
N .
2
ee
{} {}
l
β ∂
∂∂ qq
,
2
ee
{} {}
m
β ∂
∂∂ qq
and
2
ee
{} {}
n
β ∂
∂∂ qq
in Eq. (B5) are
obtained as follows.
2
ee
{} {}
l
β ∂
∂∂ qq
is
2
ee
2
3 3 31 3 3 3 4 31
34 4
4 1
1 3 1 4
3 3 3 4 31
4
{} {}
[]
[ ] [] (1 )[] [ ] [] { }
00
[ ] [ ]{ }
.0
l
I
I T II I T II I T II
II T I
N
N NL NL
NL
symm
β
ε
× ×× × ×
×
××
×× ×
∂
=
∂∂
′
′′ ′
+
′
l
l1 l l l l
l2 l
qq
0 Δ0 0 0 0
T0 Δ ΔΔ Δ l
00
00 0
T Δl
(B7)
with
13 4 1
13 4 2
13 4 3
13 4 4
{ } { } (1 ){ } [ ]
{ } { } (1 ){ } [ ]
[ ]2
{ } { } (1 ){ } [ ]
{ } { } (1 ){ } [ ]
T
I II II II
T
I II II II
T
I II II II
T
I II II II
N N NL
N N NL
N N NL
N N NL
ε
ε
ε
ε
′′ ′
+ ++
′′ ′
+ ++
=
′′ ′
+ ++
′′ ′
+ ++
l1
r r lG
r r lG
T
r r lG
r r lG
(B8a)
91
4 1
42
43
44
(1 ){ } [ ]
(1 ){ } [ ]
[ ]2
(1 ){ } [ ]
(1 ){ } [ ]
T
II I
T
II I
T
II I
T
II I
NL
NL
NL
NL
ε
ε
ε
ε
′
+
′
+
=
′
+
′
+
l2
lG
lG
T
lG
lG
(B8b)
where,
1
2
3
4
1 000
[ ] 0 0 0 1,
0 0 10
01 0 0
[ ] 0 0 1 0,
000 1
0 0 10
[ ] 0 1 0 0,
1 000
000 1
[ ] 1 00 0
01 0 0
=
−
=
−
=
−
−
=
G
G
G
G
(B9)
2
ee
{} {}
m
β ∂
∂∂ qq
is
1
3 3 3 1 3 3 3 5
3
2 4 1 4 5
1 3 1 5
ee
3 3 3 5
5 5
[]
[] [ ]
0
{} {}
.
I
I T
m
N
N
symm
β
× ×× ×
××
××
××
×
′
′
∂
=
∂∂
m
mm
0 Δ0 0 0
T0 Δ0
00
qq
00
0
(B10)
with
92
13 4
13 3
13 2
13 1
{} { } [ ]
{} { } [ ]
[ ]2
{} { } [ ]
{} { } [ ]
T
I II
T
I II
T
I II
T
I II
NN
NN
NN
NN
′′
+
′′
−+
=
′′
+
′′
−+
m
r rG
r rG
T
r rG
r rG
(B11)
2
ee
{} {}
n
β ∂
∂∂ qq
is
1
3 3 3 1 3 3 3 5
3
2 4 1 4 5
1 3 1 5
ee
3 3 3 5
5 5
[]
[] [ ]
0
{} {}
.
I
IT
n
N
N
symm
β
× × × ×
××
××
××
×
′
′
∂
=
∂∂
n
n n
0 Δ0 0 0
T0 Δ0
00
qq
00
0
(B12)
with
13 3
13 4
13 1
13 2
{} { } [ ]
{} { } [ ]
[ ]2
{} { } [ ]
{} { } [ ]
T
I II
T
I II
T
I II
T
I II
NN
NN
NN
NN
′′
−+
′′
−+
=
′′
+
′′
+
n
r rG
r rG
T
r rG
r rG
(B13)
Then, the second derivatives of π ’s are expressed as
2
0
3
2
ee
14
{} {} { }
81
81
T
l
l
π
β
β
γ
γ
∂ ∂
= −
∂∂ ∂
+
+
l
l l
e
F
FF
qq q
(B14a)
93
2
2
3
2
ee
5
2
1
{} {} { }
21
3 1
{}
21
21
T
n
l
T
n
l
l
β π
γ
β
γ
β
γ
β
β
γ
γ
∂ ∂
= +
∂∂ ∂
+
∂
−−
∂
+
+
l
l n
e
n
l l
e
F
FF
qq q
F
FF
q
(B14b)
2
3
3
2
ee
5
2
1
{} {} { }
21
31
{}
21
21
T
m
l
T
m
l
l
π β
γ
β
γ
β
γ
β
β
γ
γ
∂ ∂
= −−
∂∂ ∂
+
∂
++
∂
+
+
l
l m
e
m
l l
e
F
FF
qq q
F
FF
q
(B14c)
94
C1.
1
11
2
22
3
33
4
44
[] { }
{ } [ ] [ ][ ] [ ][ ]
{} {}
[] { }
{ } [ ] [ ] [ ] [ ] [ ]
{} {}
{ }
{}
[] { }
{ } [ ] [ ][ ] [ ][ ]
{} {}
[] {
{ } [ ] [ ] [ ] [ ] [ ]
{}
I
TT I I
I
TT I I
I
TT I I
I
TT I I
∂ ∂
+ +
∂∂
∂ ∂
++
∂∂
∂
=
∂
∂ ∂
++
∂∂
∂ ∂
++
∂
r
rr
ee
r
r r
ee
e
r
r r
ee
r
r r
e
P b
b Π PD PΠ
qq
P b
b Π PD PΠ
qq
p
q
P b
b Π PD PΠ
qq
P b
b Π PD PΠ
q
1
2
3
4
,
}
{}
[]
{ }
{}
[]
{ }
{}
with [ ]
[]
{ }
{}
[]
{ }
{}
T
T
T
T
∂
∂
∂
∂
∂
=
∂
∂
∂
∂
e
r
e
r
e
r
e
r
e
q
Π
b
q
Π
b
q
D
Π
b
q
Π
b
q
(C1)
where, []
I
i r
P and []
i r
Π are i th rows of matrices []
I
P and [] Π respectively.
95
C2.
2
e
1 e 1
2
e
2e 2
2
e
16 e 16
2
1
e
e
{} {}
{ }[ ] [ ]
{}
{} {}
{ }[ ] [ ]
[] ,
{}
{} {}
{ }[ ] [ ]
{}
{}
{
with [ ]
T
q
TT
ee
T
q
TT
ee
T
q
TT
ee
T
q q
qq
qq
r
∂∂
+
∂∂ ∂
∂∂
+
=
∂∂ ∂
∂∂
+
∂∂ ∂
∂
∂
=
v
v
1
v
v
r r
qr D
q
rr
qr D
V
q
rr
qr D
q
q
q
D
e
2
2
e
ee
2
3
e
ee
}{ }
{}
{} {}
{}
{} {}
T
T
r
r
∂
∂
∂∂
∂
∂∂
q
q
qq
q
qq
(C2)
96
C3.
2
e
1 e
11 1
2
e
2e
22
{ }
{ } [ ][ ][ ][ ]
{}
{ } { } { }
[ ] [ ] [ ][ ] [ ] [ ][ ][ ]
{ }
{ } [ ][ ][ ][ ]
{}
[ ]
{ } { }
[]
q
TT T T
e
TT T
TT
ee e
q
TT T T
e
TT
ee
q
qq q
q
qq
∂
+
∂∂
∂∂ ∂
++
∂∂ ∂
∂
+
∂∂
=
∂∂
+
∂∂
ab c
2
a
p
q pE J E
q
pp p
D EJ D EJ E D
p
q pE J E
q
V
pp
D
2
2
e
16 e
16 16 16
{ }
[ ] [ ][ ] [ ] [ ][ ][ ]
{ }
{ } [ ][ ][ ][ ]
{}
{ } {} { }
[ ] [ ] [ ][ ] [ ] [ ][ ][ ]
T
TT
e
q
TT T T
e
TT T
T T
e e e
q
q
qq q
∂
+
∂
∂
+
∂∂
∂∂ ∂
+ +
∂∂ ∂
bc
ab c
p
EJ D EJ E D
p
q pE J E
q
pp p
D EJ D EJ E D
(C3)
with
97
1
e
e
2
e
e
3
e
e
4
e
e
1
e
e
2
e
e
[]
{ } [ ][ ][ ]
{}
[]
{ } [ ][ ][ ]
{}
[] ,
[]
{ } [ ][ ][ ]
{}
[]
{ } [ ][ ][ ]
{}
[]
{ }[ ]
{}
[]
[ ] { }[ ]
{
q
TT T T
q
TT T T
q
TT T T
q
TT T T
q
TT
q
TT
∂
∂
∂
∂
=
∂
∂
∂
∂
∂
∂
∂
=
∂
c
c
a
c
c
r
r
b
E
q pE J
q
E
q pE J
q
D
E
q pE J
q
E
q pE J
q
E
qp
q
E
D qp
q
3
e
e
2
0
e
ee
2
1
e
ee
2
2
e
ee
2
3
e
ee
,
}
[]
{ }[ ]
{}
{}
{} {}
{}
{} {}
[]
{}
{} {}
{}
{} {}
q
TT
T
T
T
T
p
p
p
p
∂
∂
∂
∂∂
∂
∂∂
=
∂
∂∂
∂
∂∂
r
c
E
qp
q
q
qq
q
qq
D
q
qq
q
qq
(C4)
where, []
i c
E and []
i r
E are i th columns and i th rows of matrix [] E respectively.
98
C4.
2
1
ee
1e
1
2
2e
ee
16
16 e
{} {}
[] []
{ }
{}
{ }[ ] [ ]
[] []
{ } { }[ ]
[] { }
{} []
{} {}
{ }[ ]
[] []
{ }
{}
i
I T
T i
e
TI
I T
T
T i T
I
i
e
T
I T
T i
e
p
q
q
q
−
∂
=
∂∂
∂ ∂
∂ ∂
∂ ∂
∂ ∂
∂∂ ++ +
∂∂
∂ ∂
∂∂
r
pr
r
p
r
p
r
qq
P Π
b
q
bD P
P Π
b bD
P b
q Π
qq
bD
P Π
b
q
1
2
16
2
1 e
1
2
2
2e
e e
16
2
16 e
[ ]
[ ][ ]
[ ][ ]
{ }
[ ][ ]
{}
[ ][ ]
{ }
[ ][ ] [ ][ ]
[] { }
{} []
{} {}
[ ][ ]
{ }
[ ][ ]
{}
i
I
i
I
i
I
i
e
I
i
T
I I
IT
i i T i
e
I
i
I
i
e
q
q
q
∂
∂∂
∂
∂ ∂
∂∂ + ++
∂∂
∂
∂∂
π
r π
r π
r
rb
r rb
r
rb
r
D
PD
PD
b
P Π
q
PD
b
P Π PD
P b
q Π
qq
PD
b
P Π
q
(C5)
with
1
e
2
e
4
e
[ ] []
{}
[ ] []
{} []
[ ] []
{}
I
i
ej
I
i
ej j
I
i
ej
q
q
q
∂∂
∂∂
∂∂
∂∂ =
∂∂
∂∂
rc
rc
p
rc
P Π
q
P Π
q D
P Π
q
(C6)
99
2
11
ee
2
22
ee
2
16 16
ee
[] [] { }
{ }
{} {}
[] [] { }
{ }
{} {} []
[] [] { }
{ }
{} {}
T
ej ej
T
ej ej j
T
ej ej
qq
qq
q q
∂∂ ∂
+
∂ ∂ ∂∂
∂∂ ∂
+
∂ ∂ ∂∂ =
∂∂ ∂
+
∂ ∂ ∂∂
rr
r r
π
rr
ΠΠ b
b
qq
ΠΠ b
b
qq D
ΠΠ b
b
qq
(C7)
1
e
2
e
4
e
[] { }
{}
[] { }
{} []
[] { }
{}
T
ej
T
ej j
T
ej
q
q
q
∂ ∂
∂∂
∂ ∂
∂∂ =
∂ ∂
∂∂
r
r
b
r
Π b
q
Π b
q D
Π b
q
(C8)
where, 1,2,...,4 and 1,2,...,16 ij = = .
Abstract (if available)
Abstract
In this study, a three‐dimensional dynamic modeling and simulation of a satellite that consists of coupled rigid‐flexible bodies is presented. In the modeling of this coupled system, the satellite is assumed to be a rigid body
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
An approach to experimentally based modeling and simulation of human motion
PDF
Modeling and dynamic analysis of coupled structure-moving subsystem problem
PDF
On the synthesis of dynamics and control for complex multibody systems
PDF
Flexible formation configuration for terrain following flight: formation keeping constraints
PDF
Dynamic analysis and control of one-dimensional distributed parameter systems
PDF
Mathematical model and numerical analysis of a shape memory polymer composite cantilever beam for space applications
PDF
Transient modeling, dynamic analysis, and feedback control of the Inductrack Maglev system
PDF
Control of spacecraft with flexible structures using pulse-modulated thrusters
PDF
Terrain-following motion of an autonomous agent as means of motion planning in the unknown environment
PDF
An approach to dynamic modeling of percussive mechanisms
PDF
Active delay output feedback control for high performance flexible servo systems
PDF
Nonlinear control of flexible rotating system with varying velocity
PDF
An analytical dynamics approach to the control of mechanical systems
PDF
Periodic motion analysis in spacecraft attitude dynamics
PDF
Modeling and vibration analysis of wheelchair-users
PDF
The role of rigid foundation assumption in two-dimensional soil-structure interaction (SSI)
PDF
Modeling, analysis and experimental validation of flexible rotor systems with water-lubricated rubber bearings
PDF
Medium and high-frequency vibration analysis of flexible distributed parameter systems
PDF
Closed-form transient analysis of one-dimensional distributed dynamic systems
PDF
Control of two-wheel mobile platform with application to power wheelchairs
Asset Metadata
Creator
Choi, Homin
(author)
Core Title
Dynamic modeling and simulations of rigid-flexible coupled systems using quaternion dynamics
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Aerospace and Mechanical Engineering (Dynamics and Control)
Publication Date
06/20/2014
Defense Date
05/05/2014
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
Euler angle singularity,Euler angles,multi-body dynamics,OAI-PMH Harvest,quaternion singularity,quaternion-based FEM,quaternions,rigid body dynamics,rigid-flexible coupled systems,satellite orientation,torque induced singularity
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Yang, Bingen (Ben) (
committee chair
), Flashner, Henryk (
committee member
), Wang, Chunming (
committee member
)
Creator Email
homincho@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-421237
Unique identifier
UC11286053
Identifier
etd-ChoiHomin-2568.pdf (filename),usctheses-c3-421237 (legacy record id)
Legacy Identifier
etd-ChoiHomin-2568.pdf
Dmrecord
421237
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Choi, Homin
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
Euler angle singularity
Euler angles
multi-body dynamics
quaternion singularity
quaternion-based FEM
quaternions
rigid body dynamics
rigid-flexible coupled systems
satellite orientation
torque induced singularity