close

Вход

Забыли?

вход по аккаунту

?

6.2011-6277

код для вставкиСкачать
AIAA 2011-6277
AIAA Modeling and Simulation Technologies Conference
08 - 11 August 2011, Portland, Oregon
Reduced Order Multi-body Problems
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Peter M. Thompson1
Systems Technology, Inc., Hawthorne, CA, 90250
A multi-body problem has six degrees of freedom for each body. Constraints are used to
restrict motion, for example by defining joints that tie the bodies together. These constraints
reduce the total degrees of freedom. Two methods are described to create a transformation
matrix used to create a reduced order, unconstrained multi-body problem. One method
starts from the constraint function, computes the Jacobian of the constraint function, and
then the desired transformation matrix is the null space of the Jacobian. The second method
is more direct and creates the same transformation matrix starting from the velocity
relationships that define the joints. Equivalence of the two methods is proved. Closed form
expressions are given of the transformation matrix, which can be used either for numerical
simulation or as a starting for further symbolic expansion. Methods for linearization of the
multi-body problem are discussed. This paper is tutorial and starts by defining a general and
precise notation. The Newton-Euler dynamics problem is reviewed, a catalog is presented of
constraint functions used for joints, and then a catalog of generalized forces appropriate for
aerospace problems. Examples are presented using combinations of rigid bodies including
an overview of satellite, rotorcraft, and wind turbine problems. The examples in this paper
are considered complete when the transformation matrix used to create the reduced order
problem is defined.
I.
Introduction
M
ulti-body dynamics are used for problems with connected bodies. Classic examples are a rotorcraft with a hub
and blades1,2,3, tilt rotors4,56,7, wind turbines8,9, a satellite with appendages10,11, and bio-dynamic problems12.
Multi-body methods use algebraic constraints that define joints, which are then used to automatically determine the
equal and opposite forces and torques that tie together the different bodies. Using algebraic constraint functions is a
more straightforward approach than the alternative of directly determining the constraint forces and torques.
The multi-body reference by Haug13 has had the most influence on this work. His high level description of a full
order, constrained multi-body system is adopted here, which is:
Mq  Tq   Q
The terms M , q, and Q are respectively the generalized mass, acceleration, and force. The term q is Jacobian of
the constraint function, and the combined term Tq  is the generalized constraint force. Haug’s work, which is
typical of many approaches to multi-body problems, is to develop and implement simulation methods that work
directly with the constrained, full order system.
The objective here is different, which is to extract a reduced order, unconstrained version of the multi-body
problem. If a transformation matrix H2 can somehow be found such that  q H 2  0 and q  H 2 q2 , where q2 is a
reduced order velocity vector, then the following reduced order multi-body system can be created:
( H 2T MH 2 )q2  ( H 2T Q)  ( H 2T MH 2 )q2
Two methods for creating the H2 transformation matrix are described in this paper and then demonstrated. The
first approach starts with  q H 2  0 and finds the null space. A numerical solution is not sufficient because the
1
Chief Scientist, Systems Tech., Inc., 13766 S. Hawthorne Blvd. Hawthorne, CA 90250. Senior Member AIAA.
Copyright © 2011 by Systems Technology, Inc. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
derivative of H2 is also needed. Closed form expressions for the null space, from which the derivative is computed,
are found for individual problems, including significant problems such as rotorcraft systems, but a more general
method is needed. The second approach arranges the velocity transformations at the joints into the matrix equation
H a q  Hb q2 from which it follows that H 2  H a1Hb and H 2  H a1 ( Hb  H a H 2 ) . Both the transformation
matrix H2 and its derivative are available, and so the second approach generalizes to arbitrarily large multi-body
systems. The difficulty is that it needs to be shown (as is done in Section IV) that the H2 matrix constructed in this
way is also the null space of the Jacobian.
This paper is tutorial and has a second objective of describing and promoting a general notation for multi-body
problems. This is done in Section II followed by comparisons with notation used by Haug13, Kane14, and Britting15.
The recommended notation includes more information than is typically done in the descriptions of position and
velocity vectors, in particular the relative points and frames that used to define the vectors, the coordinate frame, and
the frame in which the derivatives are computed. Most treatments assume this information is known from context,
which in the multi-body problem quickly leads to confusion. The notation includes treatment of Euler angles,
direction cosine matrices, and perturbations.
Section III reviews the single-body and multi-body dynamic problems, and then describes the main contribution
of this paper, which are the methods for computing the H2 transformation matrix. The linear dynamic problem is
also discussed.
Sections IV and V respectively give descriptions of constraints and generalized forces needed for aerospace
problems. The algebraic description of the constraints used for joints is needed in order to prove that the H2
transformation matrix created from the joint velocity relationships is also the null space of the constraint Jacobian.
Section VI gives a set of examples, starting with two familiar single body examples that are included in order to
build confidence. Different connections of multiple bodies are then treated, those being multiple appendages off of
one body, and connected joints. These two arrangements of multiple bodies are then used to describe the rotorcraft
and wind turbine problems. The examples are not carried all the way through to numerical simulations, but rather
are considered “complete” when a closed form description is given of the H2 transformation matrix.
Quaternions using the Euler-Rodrigues formulation are presented in the Appendix, using vector and matrix
notation consistent with that presented in Section II.
II.
Notation
A general notation is defined for single-body and multi-body dynamics problems. Discussion, examples, and
comparison with other notational schemes are included.
Points, frames, and rigid bodies:
A point is a location in three-space. A frame is a point and a right-handed axis system with its origin at the point.
An inertial frame does not translate or rotate. A rigid body is defined by a point fixed with respect to the body,
which is usually but not necessarily at the center-of-mass (cm) of the rigid body. The notation is:
b  frame, point, or rigid body (letter, subscripted letter, or number)
e  inertial frame (the letter e is used here for the inertial frame)
Diagrams of a frame, a position vector, and a rigid body are shown in Figure 1.
xb
b
b
b
yb
zb
e
reb / e
rbp / b
reb / e
p
e
b) Position vector
a) Frame
c) Rigid Body
Figure 1: Frame, Point, and Rigid Body
2
Vectors:
A position vector represents three quantities: a point, a second point to which the first point is relative, and a
frame in which the Cartesian coordinates of the first point are defined. Vector notation is defined below:
rxb/ a  position of point b relative to point a in x-frame Cartesian coordinates (31)
vbx / a  velocity of point b relative to point a in x-frame Cartesian coordinates (3 1)
abx / a  acceleration of point b relative to point a in x-frame Cartesian coordinates (31)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
xb/ a  angular rate of frame b relative to frame a in x-frame Cartesian coordinates (31)
Fxname  force with given name in x-frame Cartesian coordinates (3 1)
nxname  torque with given name in x-frame Cartesian coordinates (3 1)
I  [e1 e2 e3 ]  unit vectors respectively along x, y, and z-axes (each 31)
xb , yb , zb  unit vectors at point b respectively along x, y, and z-axes (each 31)
f xa , g xa , hxa  vectors at point a (unit norm, usually but not necessarily orthogonal, each 31)
In the vector notation superscripts are reserved for points and names, and subscripts for coordinate frame. If the
subscript is omitted the notation becomes coordinate free, but the preference here is to always include the subscript.
The slash used to separate points and frames can be read “relative to.” If a transpose is needed surround the term in
parentheses, for example (rxb/a)T. To help the reader the dimensions of the various terms are included in parentheses.
The following names of forces and torques are suggested: the name of the body followed by a single letter:
Fxbt , nbt
x  total force, torque on rigid body b (3 1)
Fxbw  weight force on rigid body b (3 1)
Fxba , nba
x  aero force, torque on rigid body b (3 1)
Fxbe , nbe
x  elastic force, torque (spring and damping) on rigid body b (3 1)
nbbg  bb / e Jbbb / e  gyroscopic torque on rigid body b (3 1)
Vector components: no general rule is given, so components should be defined before being used. Typical
components for aerospace problems are:
reb/ e  [ x y z]T  position of body b in earth frame coordinates (3 1)
vbb/ e  [u v w]T  velocity of body b in body frame coordinates (3 1)
bb/ e  [ p q r ]T  angular velocity of body b in body frame coordinates (31)
Fbaero  [ X Y Z ]T  aero force in body frame coordinates (3 1)
nbaero  [ L M
N ]T  aero torque in body frame coordinates (3 1)
Direction cosine matrices:
The orientation of one frame relative to another is defined using a direction cosine matrix:
Aa /b  direction cosine matrix of frame a relative to frame b (3  3)
Single axis rotations, where ca  cos(a ), sa  sin(a) :
0
1 0
ca 0  sa 
 ca sa 0 




Rx (a)  0 ca sa  , Ry (a)   0 1 0  , Rz (a)    sa ca 0 
0  sa ca 
 sa 0 ca 
 0
0 1 
R f (a)  f f T  f sa  f f ca  single axis rotation about unit norm vector f
Direction cosine matrices are used to change the coordinate frame of a vector, for example: rax / y  Aa /b rbx / y .
3
Cross products:
a3
 0

a   a3
0
 a2 a1
a2 
skew symmetric matrix version
a1  
of cross product, where ab  a  b
0 
Cross product identities include:
aa  0
tilde(ab)  ab  ba  baT  abT
aab  aba
a a0
abc  (cb  bc )a  (a c)b  (a b)c
aT ba  0
ab  ba
ab  baT  (aT b) I  ba  (baT  abT )
babc  b (bT a )c
aaa  a(aT a)
tilde( Aa)  AaAT
baab  abba
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
T
T
T
The following combination of a cross-product and a change of reference frame occurs frequently and is given a
special name:
bxa /b   Ae/ x rxa /b (3  3)
Note the sub- and super-scripts of b and r are the same. This term converts force at point a to torque about point b:
nxname  (bxa /b )T Fename  rxa /b Fxname (3 1)
Derivatives:
The derivatives of vectors and direction cosine matrices are defined below:
rxa /b  derivative computed in x-frame (3 1)
rxa /b  two derivatives each computed in x-frame (3 1)
vea /b  rea /b  the letter v is reserved for derivatives computed in the inertial frame (31)
aea /b  vea /b  rea /b  the letter a is reserved for acclerations computed in the inertial frame (31)
Aa /b  Aa /bbb/ a  ab/ a Aa /b  derivative of a direction cosine matrix (3  3)
Aa /b  Aa /b (bb/ a  bb/ abb/ a )  (ab/ a  ab/ aab/ a ) Aa /b (3  3)
The frame in which a derivative is computed needs to be defined. The letters v and a are reserved for derivatives
computed in an inertial frame. If the subscript of v or a is changed, this means the derivative is first computed in the
inertial frame, and then the coordinate frame is changed. When the dot notation is used the subscript serves double
duty by defining the frame in which the derivative is computed and the coordinate frame used to represent the
derivative. A consequence of this double use is the subscript of a dotted quantity cannot be changed. The rules to
follow are:
vxa /b  Ax / e vea /b (subscript of v, a, , F , and n vectors can be changed)
rya /b  Ay / x rxa /b (subscript of a dotted vector cannot be changed)
Velocity computed in inertial and rotating frames are related as follows:
vxa /b  rxa /b  xx / e rxa /b  velocity in a rotating frame
The derivation uses the chain rule for differentiation:
vxa /b  Ax / e
4
d ( A r a /b ) 
dt e / x x
Ax / e ( Ae/ x rxa /b  Ae / x rxa /b )  rxa /b  xx / e rxa /b
Acceleration is simple in the inertial frame:
aeb /e  reb /e
But is complicated in general:
abb /e  rbb /e  2bb /e rbb /e  bb /e rbb /e  bb /ebb /e rbb /e  rbb /e  bb /e rbb /e  bb /e rbb /e  bb /e (rbb /e  bb /e rbb /e )
Coriolis
Tangential
Centripetal
vbb /e
vbb /e
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Derivatives of a single axis direction cosine matrix are:
R f (a)   R f (a) f a   f aR f (a)
R f (a)  R f (a)( f a  f f a 2 )  ( f a  f f a 2 ) R f (a)
Euler angles:
Euler angles are single axis rotations sequentially about one, two, or three axes. The number of rotations, the
axes, the order, and the names depend on context and should be defined before being used. A general notation is
defined below that relates angular velocity and Euler angle rates:
bb/ a  Bbb/ ab/ a  angular velocity (3 1), where:
Bbb/ a  matrix whose columns are unit-norm axes of rotation in the b-frame (3  L)
b/ a  vector with corresponding Euler angles (L 1)
L  number of Euler angles
Some examples: 1) rotation by β about the y-axis, 2) rotation by  and β respectively about the z and y-axes, 3)
the usual aerospace convention:
Ab / a  R y (  )
Ab / a  R y (  ) Rz ( )
Ab / e  Rx ( ) R y ( ) Rz ( )
bb / a  Bbb / ab / a , where:
bb / a  Bbb / ab / a , where:
bb / e  Bbb / eb / e , where:
Bbb / a  e2
Bbb / a  R y [e2
Bbb / e  Rx [e1 e2
b / a  
b / a  [   ]T
e3 ]
R y e3 ]
b / e  [   ]T
The orthogonal complement to the axes of rotation plays a role in subsequent analysis, defined by:
Ybb/ a  matrix whose unit-norm columns are orthogonal to Bbb/ a (3  (3  L)), where:
(Ybb/ a )T Bbb / a  0 ((3  L)  L)
b/e
When there are three Euler angles there is gimbal lock condition where Bb
is not invertible. For three
b/ e
rotations β respectively about uvw the transformation matrix is Bb  Rw ( )[w v Rv (  )u] . The matrix
Bbb / e is not invertible if the vector Ry (  )u is in the plane formed by vectors w and v, equivalently if
(vw)T Ry (  )u  0 , and equivalently if   atan2( wT vu, wT vvu) , where atan2 if the four-quadrant arctangent. For
rotations  respectively about zyx the gimbal lock condition is  = /2. See the Appendix for a discussion of the
use of quaternions as an alternative to Euler angles.
Perturbations:
All of the vectors and angles above are total quantities. A variable name preceded by  is a perturbation. The
total quantity equals the equilibrium value plus a perturbation:
x  x0   x
5
The equilibrium value is denoted by a subscript or superscript zero; use parentheses if other scripts are present:
x0 , (reb/ e )0 , ( Ab/ e )0  equilibrium values
( X )0q  partial of X with respect to q evaluated at equilibrium
The perturbation of a direction cosine matrix13 is defined using a perturbation angle, also called a variation angle,
which is defined as a small rotation:
 Aa /b  ( Aa /b )0  bb/ a , where:
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
In the limit as  t  0,  bb/ a  (bb/ a )0  t
Variation angles are typically replaced by Euler angle perturbations, where:
 bb/ a  ( Bb / a )0 b/ a
Here is an example with an expanded version of a perturbation:
 ( Aa /b rbx / y )  ( Aa /b )0  rbx / y  ( Aa /b rbx / y Bb/ a )0 b/ a
Mass properties:
mb  mass of body b (11)
Jb  moment of inertia matrix of body b at the cm in body frame coordinates (3  3)
Jbp /b  Jb  mrbp /b rbp /b  moment of inertia matrix of body b at point p in b-frame coordinates (3  3)
J ap / b  Aa / b Jbp /b Ab / a  similarity transformation (3  3)
The moment of inertia of a rigid body is constant in its own body-fixed coordinates, which is why the body-fixed
frame is used for rotation.
Multibody dimensions:
A multi-body system is defined with N rigid bodies. Any names or numbers can be given to the bodies. The
dimensions are:
N  # rigid bodies
6 N  # unconstrained degrees of freedom
L  # constrained degrees of freedom
P  6 N  L  # unconstrained degrees of freedom
M  # inputs (such as actuator positions)
Generalized mass, velocity and force
Generalized terms combine mass and moments, translations and rotations, and forces and torques:
M  generalized mass (6N  6N )
q  generalized velocity (6N 1)
Q(q, q, u)  generalized force (6N 1)
u  inputs (M 1)
The generalized force Q may also depend on acceleration, in which case the acceleration should be added to the list
of dependencies. Expanded versions are listed below, where the rigid bodies are numbered 1 to N:
6
 m1I


M 



J1
mN I
e 
 v1/
 Fe1t


 e 



 11/ e 
 n11t  11/ e J111/ e






, q  
 , Q ( q, q, u )  


 N /e 
 Nt


 ve 
 Fe

J N 
 N / e 
 n Nt   N / e J  N / e 
N
N N 
 e 
 N
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Semicolon notation (as in Matlab code) can be used to save space when defining tall columns: for example:
e
1/ e
q  [v1/
e ; 1 ;
veN / e ; eN / e ]
The vector q
The vector q is a surrogate name for positions and orientation. The orientation can be either Euler angles or
quaternions, with the former used here. The partial with respect to q is more important than the vector q itself, which
is well defined by x  xq q . Examples such as the following occur frequently in the Jacobian (see Figure 1c):
Given: rep / e  reb / e  Ae / b rbp / b (3  1)
And: q  [veb / e ; bb / e ] (6 1), then:
d r p / e  (r p / e ) q
e
q
dt e
(rep / e )q  [(rep / e )r
(3 1), where:
(rep / e ) ]  [ I
 Ae / b rbp / b ]  [ I
bbp / b ] (3  6)
If the partials with respect to translation and rotation need to be separately called out then use subscript r for
translation and  for rotation.
Constraint function
The constraint function and related terms are:
(q)  constraint function ( L 1)
q  Jacobian ( L  6 N )
Tq   generalized constraint force (6 N 1)
  constraint vector ( L 1)
H1  right inverse of q (6 N  M )
H 2  null space of q (6 N  P)
[ H1 H 2 ]  invertible (6 N  6 N )
A more general version13 of the constraint function is (q,u), where u is an input variable. This more general
form is motivated by kinematics problems, not dynamic problems, and so is not used here. A technical condition on
the validity of the constraint function is that the Jacobian q should have full row rank, which means there are no
redundant constraints.
Reduced order system:
The terms used to represent the reduced order multi-body system are:
Mˆ 22  H 2T MH 2  generalized mass (P  P)
Nˆ 22  H 2T MH 2  gyroscopic term ( P  P)
Qˆ 2  H 2T Q  generalized force (P 1)
7
Shortcut Notation:
Using shorter versions of the names for a given problem is always acceptable but should be defined before being
used, for example:
B  Bbb/ e  axes of rotation (B is shortcut name, 3  3)
Comparison with Haug
The notation used here is compared to Haug13:
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
vb  veb/ e  velocity at cm of rigid body in inertial-frame coordinates (3 1)
vb  vbb/ e  velocity at cm of rigid body in body-fixed coordinates (3 1)
b  eb/ e  rotational rate rigid body b in inertial-frame coordinates (3 1)
b  bb/ e  rotational rate of rigid body b in body-fixed coordinates (3 1)
b  tilde(b )  cross-product matrix (same as used here) (3 3)
Ab  Ae/b  direction cosine matrix (3 3)
Haug’s notation is very efficient when defining motion at the cm of a rigid body and very efficient for defining
multi-body systems. The use of the letters q, M, Q, J, , F, and n and the tilde operator for cross products follows
from Haug. His notation is not as good for defining points and motion at locations not at the cm, and for defining
relative motion between rigid bodies. The preference here is that the subscript should always to be the coordinate
frame, and so the prime notation is not adopted.
Comparison with Kane
A comparison is made with Kane14:
A, B, C, etc  frames of reference or points
p  rxP /O  position vector (3 1)
v 
B P
a 
B P
A
B
dp
 Ax / B rBP /O  velocity (3 1)
dt
B B P
d v
 Ax / B rBP / O  acceleration (3 1)
dt
ωB  xB / A  angular velocity (3 1)
Kane’s notation is coordinate free. Cross products use the traditional × notation. A derivative in Kane’s notation
indicates the frame in which the derivative is computed but not the coordinate system in which the vector is
represented. In the above comparison the coordinate frame x is unspecified. The two points that define a position
vector need to be understood from context, as does the relative point for a velocity vector. Direction cosine matrices
are not introduced in Kane except as a problem in Chapter 5. Kane’s work is a classic and beyond reproach, but
some find it … clumsy.
Comparison with Britting
Textbooks in inertial navigation such as Britting15 are very precise in their notation. Here is a comparison:
a, b, c  points or frames of reference
i, I  operational, absolute inertial frames
r x  rxa /b  position vector (3 1)
r x  rxa /b  velocity (3 1)
r x  rxa /b  acceleration (3 1)
8
C yx  Ax / y  direction cosine matrix (3  3)
x
ωab
 xa /b  angular velocity (3 1)
x
Ωab
 xa /b  skew symmetric cross product matrix (3  3)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
The superscript is used for the coordinate frame, and the superscript in derivatives serves double duty by
indicating the frame in which the derivative is computed. In Britting the tilde operator indicates a measured quantity,
and capital omega or superscript * indicates transformation to the skew-symmetric cross product matrix. There is no
attempt in Britting to generalize to multi-body systems. No complaints – the differences in the notation used here are
personal preferences.
Personal preferences
Notation inevitably includes personal preference. The preferences here are 1) subscripts are reserved for
coordinate frames, 3) superscripts are reserved for points and names, 3) relative points and frames are not assumed
but included as part of the notation with a slash to help distinguish, 4) no use of bold fonts or arrows, which is
considered old fashioned, 5) no superscripts or subscripts on the left, which is considered intimidating, 6) skew
symmetric matrices are used in place of cross-products, 7) flexibility is allowed in naming components, forces and
torques, and in using shortcuts, but the user is asked to define these terms before using them. Some purists may
consider the use of a subscript for the coordinate frame at best unnecessary or at worst a waste of space. The author
is willing to admit the use of this subscript is a crutch, but it is not without justification, because the ultimate goal is
to get things onto a computer.
The examples in Table 1 help to explain and also to justify the notation.
Table 1: Notation Examples
rda /b

rca / b
 rca / d
Changing coordinates multiplies subscripts: d  (d / c)  c . The
superscripts do not change. Coordinates can be changed for r, v, a,
, F, and n vectors.
Ad / c rca /b
Adding multiplies superscripts: (a / b)  (a / d )  (d / b) . The
subscripts must be the same. This rule does not apply to F and n,
since their superscripts refer to names, not frames.
 rcd /b
rca /b  rcb / a
Reversing superscripts changes the sign.
Aa /b  AbT/ a
Transposing transformations reverses the subscripts.
Aa / b  Aa / d Ad / b
Multiplying transformations multiplies subscripts:
(a / b)  (a / d )  (d / b)
Aa /b  Aa /bbb/ a  ab / a Aa /b
For rotational velocity, reverse the subscript of A and match the
frame on the side of A where the rotation is written.
rxa / b  Ax / y rya /b Ay / x
Similarity transformation changes the subscript from y to x.
rya / b  Ay / x rxa / b
Subscript of a dotted vector cannot be changed
III.
Newton Euler Dynamics
A. Unconstrained single rigid body
The Newton-Euler form of the dynamic equation for a single rigid body is shown below 13:
 mb I
 0

9

0   veb / e  
Febt


 (6 1)

J b   b / e    b / e J  b / e  nbt 
b 
 b   b b b
(1)
The point b is at the cm of the rigid body, and the velocities, the moment of inertia and the torques are defined at
that point. The coordinate systems for the translations and rotations are respectively in the inertial and body-fixed
frames. This is the usual convention used for multi-body problems. The position vector and Euler angles updates are
computed from:
0   reb / e   veb / e 
I
   b / e  (6 1)
0 B  
b / e  b / e  b 

(2)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Quaternions can be used as an alternative to Euler angles as described in the Appendix.
Unconstrained single rigid body (alternative forms):
In aerospace problems the aerodynamics forces and the sensor responses are more naturally expressed in the
body-fixed frame, in which case:
 mb I
 0

0   vbb / e   mbbb / e vbb / e  Fbbt 


 (6 1)
J b   b / e    b / e J  b / e  nbt 
b 
 b   b b b
(3)
For appendages attached at a hinge point, such as rotor blades, centering the dynamics at a point p not at the cm
is a good alternative:

mb I

 mb r p / b Ab / e
b

mb Ae / b rbp / b   vep / e  mbbb / ebb / e rbp / b  Febt 

 (6 1)

J bp / b  bb / e   bb / e J bp / bbb / e  nbpt 
(4)
Where:
vep / e  veb / e  Ae /bbb / e rbp /b  velocity at point p (3 1)
nbpt  nbbt  rbp / b Fbbt  total torque at point p (3 1)
This alternative though good is not adopted here because the generalized mass matrix is no longer constant.
Generalized notation and coordinate transformations:
The Newton-Euler equation (1) is written as follows in generalized form:
Mq  Q (6 1)
Alternative versions (3) and (4) are derived using the he velocity transformation q  H 2 q2 , where:
A
H 2   e/b
 0
I
0
, H2  

I
0
Ae / b rbp / b   I bbp / b 

 (each 6  6)
I
I 
 0
A more general version of this H2 transformation is used in the multi-body problem to compute the unconstrained
dynamic equations.
Constrained multi-body dynamics:
The motion of the rigid bodies is constrained using a constraint function. The constraint function and its
derivatives are:
(q)  0 (L 1)
  q q  0 (L 1)
(5)
  q q  q q  0 ( L 1)
(6)
The constrained dynamic equation for a set of rigid bodies13 is:
M q  Tq   Q (6 N 1)
10
(7)
The combined dynamic equation (7) and constraint function (6) is:
M

 q
Tq   q   Q 
   
 (6 N  L) 1
0      q q 
(8)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Multi-body dynamics software programs work directly with equation (8), which is a combined differential and
algebraic equation, numerically inverting the left hand side using the LU matrix decomposition at each step of the
simulation. Artificial damping is sometimes added to improve numerical stability.
Unconstrained multi-body dynamics
The following matrix decomposition of the Jacobian is used to create the reduced order system:
q [ H1 H 2 ]  [ I 0]  Jacobian matrix decomponsition ( L  6N )
(9)
Use the matrix terms to define the velocity transformation:
q 
q  [ H1 H 2 ]  1  (6N 1), where:
 q2 
q1  0  constrained velocity (L 1)
(10)
(11)
q2  H 2 q  unconstrained velocity ( P 1)
(12)
Equation (11) follows from   q q  q1  0 , and so q can be recreated just from just q2 , and similarly the
acceleration:
q  H2q2 (6N 1)
(13)
q  H 2 q2  H 2 q2 (6N 1)
(14)
T
Substitute (14) into (7) and multiply on the left by H 2 , which cancels the generalized constraint force, and results
in the following reduced order, unconstrained, multi-body dynamic equation:
Mˆ 22 q2  Qˆ 2  Nˆ 22 q2 (P 1)
(15)
Equation (15) depends on the generalized mass M, the generalized force Q, and the transformation matrix H2. The
matrices H1 and H2 are not unique. The H2 matrix can be multiplied on the right by any invertible P×P matrix. This
extra freedom can be used for example to change the coordinate frame of a rigid body velocity vector from the
inertial frame to a body-fixed frame, as is usually done in aerospace problems.
A more direct way to compute the H2 matrix:
The multi-body problem is considered solved if a closed form expression can be found for the H2 matrix. The
above method for doing so is summarized as follow: 1) define the constraint function, 2) create the Jacobian, and
then 3) create the null space of the Jacobian. A more direct way to create the H2 matrix is to arrange the velocities
for the rigid body system into a matrix equation of the form:
H a q  Hb q2
(16)
From which it follows that:
H 2  H a1Hb
H 2  H a1 ( Hb
(17)
 Ha H2 )
(18)
Equations (15) through (18) are the main results of this paper. It needs to be shown 1) how to build the Ha and Hb
matrices and 2) the resulting H2 matrix satisfies qH2 = 0. This will be done in Section IV after defining joint
constraints, but for now continue with the discussion of the dynamic equations.
11
The constraint vector:
The matrix H1 is not needed to create the unconstrained subsystem. It is needed for the more general problem
where the constraint function has the form (q,u), and also if the constraint vector  is desired. The constraint vector
is found by:
  H1T (Q  M q)  constraint vector ( L 1)
(19)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Implicit form (#1)
Equations (7) and (14) can be combined into the following implicit form:
 H 2T M

 I
0   q   H 2T Q 
   
 (6 N  P) 1
 H 2   q2   H 2 q2 
(20)
This is considered an implicit form because it is not reduced to the lowest order. It is not efficient for nonlinear
simulation but can be a convenient starting point for linearization. It is convenient because the first set of equations
can be left in terms of the original generalized coordinates. Convert to equation (7) by multiplying on the left by
[I –H2TM]. There may be common entries in q and q2 , in which case a lower order version of (20) can be defined.
Implicit form (#2)
Another way to combine (7) and (14)is the following:
M

 I
Tq
0
q
0    Q 
    
 (12N 1)
 H 2   q   H 2 q2 
 2
(21)
This form is useful if closed form expressions are available for the constraint force Tq  . Return to equation (7) by
multiplying on the left by H2T [I –M]. This is the approach used by Stoneking10.
Linear dynamics
Linear dynamics are perturbations about an operating point. Multi-body problems that have large motions may
not be suitable for linear analysis; nevertheless there are many problems that can be linearized about stationary,
constant velocity, rotating, or periodic operation conditions. The operating conditions are solutions of (7) or (15).
The perturbation form of (15) is:
(Mˆ 22 )0  q2  (Qˆ 2  Mˆ 22 q2  Nˆ 22 q2 )0q  q2  (Qˆ2  Mˆ 22 q2  Nˆ 22q2 )0q  q2  (Qˆ2 )u0  u (P 1)
2
2
(22)
The acceleration at the operating condition in general may not be zero and so is retained in the partials. Similar
perturbation equations can be written starting from the explicit forms (20) or (21).
Closed form expansions of the partials in (22) are not in general available and so this equation is of limited
usefulness. An alternative approach (used for aircraft and rotorcraft problems) is to expand (15) or (20) in
component form using symbolic manipulation. The result is a set of equations of the form:
 A1
 x 
A2      A3
 y 
 x 
A4   
 u 
Where:
 x  state
 y  outputs
 u  inputs
A1 , A2 , A3 , A4  partials evaluated at operating condition
12
(23)
There needs to be as many equations as states and outputs. Multiplication on the left brings the linear equations into
state space form:
 x 
 y    A1
 
1
A2 
 A3
 x 
A4   
 u 
(24)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
The outputs can include auxiliary variables, which are “extra terms” that trade off more equations for simpler
equations. The final selection of outputs can be a subset of those used in the calculations. Problems such as
rotorcraft that have periodic motion result in a linear time varying solution, which can be converted to linear time
invariant by using multi-blade coordinates.
IV.
Constraint catalog
The constraints of interest define joints with one, two, or three rotary degrees of freedom. The axes of rotation
are defined by vectors fixed on the bodies, and the angles of rotation are Euler angles. Diagrams of the three lower
level constraints are shown in Figure 2
fc
fc
e
e
rbp / b
ps
rcs / c
b

b
b

e
c
gb

gb
ps
ps
c
c
hb
a) Spherical
b) Perpendicular
c) Parallel
Figure 2: Spherical, Perpendicular, and Parallel Constraints
The following points and vectors are used to define the constraints:
b, c  points at cm of rigid bodies b, c
p, s  points fixed respectively on rigid bodies b, c
rbp / b , rcs / c  vectors fixed respectively on rigid bodies b, c (each 3 1)
gb , hb , fc  vectors fixed respectively on rigid bodies b, b, and c (each 3 1)
The generalized velocity for rigid bodies b and c is:
q  [veb/ e ; bb / e ; vec / e ; cc / e ]
For each constraint the , q , H1, H 2 , and H 2 terms are listed in Table 2. The interested reader should verify
equations (5), (9), (12), and (13) for each constraint. More details on each of the low level constraints are listed after
the table. Section IV finishes by showing how to construct the H2 for a general multi-body system and then showing
that the H2 matrix constructed in this way is a null space for the Jacobian.
Table 2: Constraint Catalogue
Spherical:
  rep / s
I 
0
H1   
I 
 
0
I
0

H2  
I

0
(123)
13
q  [I
(3  1)
0
I
b pp / b
0



s/c 
bc


I
0
0
(129)
bbp / b
I
bcs / c ] (3  12)
 I 0 0 0
 
H 2 = 0 I 0 0 
0 0 0 I 
 veb / e 


q2  bb / e 


cc / e 


(912)
(91)
Perpendicular:   gbT Ab / c fc
I 0
0 I
H2  
0 0

0 Ac / b
0
0
H1   
0
 
B
0 Y T ] (1  12)





0 B  
 veb / e 


bb / e 
q2  

 vec / e 


 c / b 
(1112)
(111)
0
0
I
(1211)
  [ gb
Parallel:
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
0
I
0
0

I
0
H 2 = 

0
0
0



B
0  B Ac / b
0
0
I
0
(121)
T
hb ] Ab / c fc
I 0
0
0 I
0
H1    H 2  
0 0
0

 
B
0 Ac / b
(122)
0
0
I
0
(2  1)
0
0 
0

B



bcs / cY 

Y 
0
0
 r p/ s 

 e
 gbT Ab / c f c 
I 0
0 I
H2  
 I bc / b
b

0 Ac / b
(124)
0 Y ] (2  12)
0
I
0
I
H 2 = 
0
0

0  B  Ac / b
0
0
I


[ gb
0
0
H1  
 I

 0
I 0
0 I
H2  
 I bc / b
b

0 Ac / b



bcs / cY 

Y 
0
0
(4  1)



bcs / c B 

B

0
0
(128)
Joint: Spherical
+ Parallel
rep / s
hb ]T
(125)
0
0
0
 q  [0 Y Ac / b
T
T





0 B  
 veb / e 


bb / e 
q2  

 vec / e 


c / b 
(1012)
(101)
(1210)
Joint: Spherical
+ Perpendicular
0
0
H1  
 I

 0
 q  [0 Y T Ac / b
(1  1)
0
0
0
I
bbp / b
q  
0 Y T Ac / b
I
I
0
 
H 2 = 0
I


0  B Ac / b
0 

0 

0 B 
 veb / e 


q2  bb / e 


c / b 


(812)
(81)
0
bcs / c 
 (4  12)
Y T 
0
0
I

bbp / b
 I bcs / c 
 (5  12)
 (5  1)  q  
0 Y T Ac / b 0
Ab / c fc 
Y T 
0

 veb / e 
I
0
0 0 



0

 
b/e 


H
=
0
I
0
0
q




2
2
b


bcs / c B 




c / b 

0

B
A
0
B
c
/
b




B

(127)
(712)
(71)
Spherical constraint (common point)
Points p and s fixed respectively on rigid bodies b and c are constrained to be equal. Three translational degrees
of freedom are constrained. The constraint function is:
  rep / s (3 1)
The expanded form is:
  (reb / e  Ae/b rbp /b )  (rec / e  Ae / c rcs / c )
The Jacobian is:
q  [ I bbp / b
Where:
14
I
bcs / c ] (3 12)
bbp /b   Ae /b rbp /b  cross product (3  3)
bcs / c   Ae/ c rcs / c  cross product (3  3)
bbp / b   Ae/bbb/ e rbp / b  derivative needed for H 2 (3  3)
bcs / c   Ae/ ccc / e rcs / c  derivative needed for H 2 (3  3)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
See Table 2 to complete the description of the spherical constraint. There are several equivalent ways to think of a
spherical constraint, each of which can be used as a starting point:
1.
2.
3.
A common point (as done above)
A generalized constraint force
A velocity constraint
The generalized constraint force imposed by the spherical constraint is interpreted as follows:
 e  constraint force applied at the common point p (inertial frame)
 p /b 
 r b  torque at b due to the constraint force
T
 q e   b

e  equal and opposite reaction force


 r s / c   torque at c due to reaction force
 c c
The derivative of the constraint function provides a constraint on the velocity:
  (veb / e  bbp /bbb / e )  (vec / e  bcs / ccc / e )  0
A more general version of the spherical constraint (which is no longer a joint) is used to restrict motion to either a
line or a plane:
  X T rep / s
Perpendicular constraint (two Euler angles)
A vector fc fixed in body c is constrained to be parallel to a vector gb fixed in body b. One degree of freedom is
constrained. The Euler angle interpretation is rotation about gb followed by rotation about fc. The constraint function
for the perpendicular constraint is:
  gbT Ab / c fc (11)
The minimal information needed to define the rotation is:
Ac /b  R fc ( ) Rgb ( )  direction cosine matrix (3  3)
B  Bcc /b  [ fc
Ac /b gb ]  axes of rotation (B is shortcut name, 3  2)
c /b  [  ]T  Euler angles (2 1)
Derived terms are:
cc /b  cc / e  Ac /bbb / e  Bcc /bc /b  relative rotational rate (3 1)
Y  Ycc /b  fc Ac /b gb  vector orthoganol to fc and gc (Y is shortcut name, 3 1)
B  [0 cc /b Ac /b gb ]  derivative needed for H 2 (3  2)
The Jacobian is:
q  [0 Y T Ac /b
15
0 Y T ] (112)
See Table 2 to complete the description of the constraint. The generalized constraint force that results from a
perpendicular constraint force is interpreted as follows:
0

 no constraint force
  A Y   equal and opposite constraint torque
Tq c   b / c c  

 no constraint force
0


 Y c  constraint torque about vectors orthogonal to axes of rotation
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
The derivative of the constraint function is:
  Y T cc /b  0
If the vectors gb and fc are perpendicular, as shown in Figure 2c , then B+ = BT. A similar constraint can be defined
where these vectors have a fixed but non-orthogonal angle:
  gbT Ab/ c fc  cos  (11)
Use B+ rather than BT in the H2+ matrix (as done in Table 2), and everything else about the development is the same.
Parallel constraint (one Euler angle)
The vectors gb and hb fixed in rigid body b are perpendicular to vector fc fixed in rigid body c. The vector fc is the
single axis of rotation between these two bodies. Two rotary degrees of freedom are constrained by this constraint.
The constraint function for the parallel constraint is:
  [ gb hb ]T Ab / c fc (2 1)
The minimal information needed to define the rotation is:
Ac /b  R fc ( )  direction cosine matrix (3  3)
B  Bcc /b  fc  axis of rotation (B is shortcut name, 3 1)
From which:
cc /b  cc / e  Ac /bbb / e  Bcc /bc /b  relative rotational rate (3 1)
c /b    Euler angle (11)
Y  Ycc /b  fc Ac /b [ gb hb ]  vectors orthoganol to fc (Y is shortcut name, 3  2)
B  0 (derivative needed for H 2 , 3 1)
Y   fccc /b Ac /b [ gb hb ]  derivative needed for H 2 (3  2)
The Jacobian is:
q  [0 Y T Ac /b
0 Y T ] (2 12)
The parallel constraint is equivalent to two perpendicular constraints, but the Euler angle interpretations are
different, and for this reason the preference here is to separately define the parallel constraint. See Table 2 to
complete the description of the constraint.
Joint constraint (spherical plus perpendicular or parallel constraint)
Combine these constraints to eliminate 4 or 5 degrees of freedom. The constraint function and its Jacobian are:
 r p/s 

 or 
 e
 gbT Ab / c fc 
[ gb
16

 (4 1 or 5 1)
hb ]T Ab / c f c 
rep / s
I
bbp / b
q  
0 Y T Ac / b
I
0
bcs / c 
 (4 12 or 5 12)
Y T 
See Table 2 to complete the description of the constraint. The combination of these two constraints introduces a new
cross-product term in the H2 matrix:
bbc /b   Ae/b rbc /b  bcc / s Ac /b  bbb/ p  cross product (3  3)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
bbc /b  bcc / s Ac /b  bcc / s Ac /b  bbb / p  derivative needed for H 2 (3  3)
The interpretation of this cross product is that force at the cm of one body is transmitted across the constrained axes
of rotation at the joint and creates a torque at the cm of the other body.
Building the H2 matrix:
The H2 matrix is used to create the reduced order multi-body system. The velocity relationships for a rigid body
system are used to construct matrices Ha and Hb, from which H2 is computed. Equations (16) and (17) are repeated
below:
H a q  Hb q2
H 2  H a1Hb
It needs to be shown that the H2 matrix created in this way is the null space of the Jacobian:
q H 2  0
This is done by constructing a third matrix Hc so that:
q  H c H a
And then showing that:
 q H 2  H c Hb  0
The dimensions are such that an H2 matrix constructed in this way will be a basis for the entire null space.
The construction process is now described. Define the rigid body system with the following dimensions:
i  1, , N  index for rigid bodies
j  1, , M  index for constraints
L  i Li  number of constraints (per body and total)
P  i Pi  reduced order dimension (per body and total)
Where for each body Li + Pi = 6 and in total L + P = 6N. The dimensions of the various matrices are:
Ha
Hb
Hc
H2
has dimension (6N  6 N )
has dimension (6N  P)
has dimension (L  6 N )
has dimension (6N  P)
q has dimension (L  6 N )
Define the following rows and blocks (blocks not defined are zero):
( Ha )i  rows corresponding to body i (6  6 N )
( Ha )ii  rows and columns corresponding to body i (6  6)
( Hb )i  rows corresponding to body i (6  P)
17
( Hb )ii  rows and columns corresponding to body i (6  Pi )
( Hc ) j  rows corresponding to joint j (Li  6 N )
( Hc ) ji  rows for joint j and columns for body i (Li  6)
(q ) j  rows corresponding to joint j (Li  6)
The velocities for each rigid body are:
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
q  [(q)i ;
q2  [(q2 )i ;
( q) N ]
(q2 ) N ]
Some of the rigid bodies will be unconstrained. If body i is unconstrained then:
(q)i  ( Hb )ii (q2 )i
And so:
( H a )ii  I (6  6)
( Hb )ii  identity or invertible transformation (6  6)
The other bodies are constrained. If body i is constrained by joint j that connects with body k, with common points pi
and sk, then the velocities relationships are:
vepi / sk  0  spherical constraint (3 1)
ii / k  Bii / ki / k  relative angular velocity (3 1)
(Yii / k )T ii / k  0  rotary constraint (Li 1)
(q2 )i  i / k  unconstrained velocity ( Pi 1)
The velocity relationships are arranged into the matrices as follows:
(vepi / sk )q 
0 
 0 
I
 , ( H b )ii  
( H a )i  
 , ( H c ) ji  
i
/
k
i
/
k T
i
/
k
 ( ) 
 Bi 
0 (Yi ) 
i
q


The corresponding rows of the constraint function are:
(q ) j  ( Hc ) ji ( Ha )i
If there are no rotary constraints for body i then B will be 3×3 and Y collapses to null. A key point is that the matrix
B (the axes of rotation in the i-frame) and Y are orthogonal, and so YTB = 0 (or collapses to null).
Special case
It is easier to show the format by example than explain it in general. For two bodies b and c connected at respective
points p and s the velocity relationships are:
b/e
0   ve   I 0


0  bb / e  0 I


I bcs / c   vcc / e  0 0


0
I   c / e  0 0
 c 
0 
0 0 I
q  Hc H a  
H
c /b T  a
0 0 0 (Yc ) 
0
I
0
I

 I b p /b
b

0  Ac / b
0
0
  b/e 
  ve 
  b / e 
 b 
 c / b 

Bcc / b  
0
0
0
Lines are included to help define the blocks. It follows that:
18
0
0 0

q H 2  H c Hb  
0
c /b T c /b 
0 0 (Yc ) Bc 
H2 is in the null space of q and the dimensions are such that H2 is a basis for the entire null space. Hence, for this
special case the H2 matrix can be used to create the reduced order multi-body system.
General case
The H2 matrix is in the null space of the jth constraint because:
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
( q ) j H 2  ( H c ) j Hb  [0 ( H c ) ji ( Hb )ii
0


0]  0, where ( H c ) ji ( Hb )ii   i / k T i / k   0
(Yi ) Bi 
This is true for each constraint, hence H2 is in the null space of the entire Jacobian, and the dimensions are such that
H2 is a basis for the entire null space. Hence, in general, the H2 matrix constructed from equations (16) and (17) can
be used to create the reduced order multi-body system.
V.
Generalized Force Catalog
A small number of generalized forces are defined, but large enough to include weight, aero, and spring-damper
forces. Diagrams are shown in Figure 3.
Fbname
fb
s

b
f( )
p
c
b
n( )
b
rep / s
ps
c
p
a) Point force
b) Translational spring-damper
c) Rotational spring-damper
Figure 3: Generalized Forces
Point force:
A vector force is applied at point p fixed on the rigid body b. Define the variables:
rbp / b  point where force is applied, fixed in rigid body b (3 1)
reb/ e  cm of rigid body b (3 1)
Fbname  applied force (3 1)
The generalized force is:
 Fename 
 (6 1)
Q
 r p / b Fbname 
b

The force due to weight is an example of a point force. The vector force is applied at the cm, and so does not create
a torque about the cm. Define g = acceleration due to gravity, and gb = 3×1 unit vector in the direction of gravity:
 m gg 
Q   b b   generalized force due to weight (6 1)
 0 
Aerodynamic force and torque are integrated across the body and modeled as a generalized force either at the cm
point b or at some other point p:
19
 F ba 
Q   e   generalized aerodynamic force applied at cm (6 1)
 nbba 


Fepa
  generalized aerodynamic force applied at point p (6 1)
Q
 n pa  r p / b Fepa 
b
 b

Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Translational spring damper:
A translational spring-damper works in a line connecting points fixed on rigid bodies. Define the problem using
the following variables:
rbp /b , rcs / c  points fixed on rigid bodies that connect translational spring-damper (3 1)
d  rep / s  line along which spring damper is defined (3 1)
 length of line (11)
u  d /  unit vector along the line (3 1)
f ( )  k (  0 )  b  spring-damper force (11)
q  [veb/ e ; bb / e ; vec / e ; cc /e ]  generalized coordinate for bodies b and c (12 1)
The generalized force is:
Q  dqT fu (12 1), where:
dq  [ I
 Ae/b rbp /b
I
Ae / c rcs / c ] (3 12)
Rotational spring-damper (about one axis)
A rotational spring-damper is defined using rotation of an Euler angle. Define the problem using the following
variables:
  fcT cc /b  transformation from relative rotational velocity to Euler angle rate (11)
n( )  k (  0 )  b  spring-damper force (11)
q  [veb/ e ; bb / e ; vec / e ; cc /e ]  generalized coordinate for bodies b and c (12 1)
The generalized force is:
Q   qT n (12 1), where:
q  fcT [0  Ac /b 0 I ] (112)
Multiple Euler Angles
The rotary spring damper can be generalized to multiple Euler angles. An example with three Euler angles is
presented. There is one spring about each Euler angle. Define the problem using:
cc /b  Bc /b (3 1)
c /b  [   ]T (3 1)
 n

n   n
n

20
  k (  u )  b  

 
   k (  u )  b   (3 1)
  k (  )  b  
u
 
  
The generalized force is:
Q  qT n (12 1), where:
q  B [0  Ac /b 0 I ] (3 12)
For structural problems the first bending mode in one or more axes can be approximated as rigid body rotation about
a point and opposed by a rotational spring-damper. Assigning successive Euler angles is arbitrary and it suffices to
use B = I.
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
VI.
Examples
The notation and methodology are demonstrated in several examples. The examples are considered “solved”
when H2 matrix is defined, a precursor either to simulation or to further symbolic expansion. The examples start
with familiar single body examples: aircraft and pendulum, in order to gain confidence with the methodology.
Diagrams for the examples are in Figure 4.
xe
e
pe
xe

ye
xb
ze
b
yb into page
a) Rigid body aircraft
b
p1
b
e
k1
zb
b) One dof pendulum
k2
p2
ze
 length
m  mass
J  diag( J xx J yy J zz )
b
xb
yb
zb
ye into page
reb / e
e
rbs / b
p3
s
h
k3
d) Connected joints
c) Four body problem
pi
h
ki
p
ki
t
pi
n
h
s
b
t
e) Rotorcraft
g
f) Wind turbine problem
Figure 4: Examples
21
k
Rigid body aircraft
The aircraft problem has a single unconstrained rigid body as shown in Figure 4a. The dynamic equation with
velocity in the body-fixed frame and including the Euler angle rates is:
 mI
0

 0
0
J
0
bw
 ba

0   v   Fb  Fb  m v 

0      nbba   J 


B    



Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Where:
v  vbb/ e  velocity (v is shortcut name, 3 1)
  bb/ e  rotational rate ( is shortcut name, 3 1)
  b/ e  Euler angles ( is shortcut name, 31)
B  Bbb/ e  rotational axes (B is shortcut name, 3  3)
m, J  mass properties (11 and 3 3)
Fbba , nbba  aerodynamic force, torque (each 3 1)
Fbbw  weight (3 1)
A linear version is expanded about non-acceleration equilibrium condition. Perturbations are:
 Fbba  ( Fbba )0r  v  ( Fbba )0   ( Fbba )u0  u
 nbba  (nbba )0r  v  (nbba )0   (nbba )u0  u
 Fbbw   ( Ab/ e Febw )   bb/ e ( Fbbw )0  ( Fbbw )0  bb/ e  ( Fbbw B)0 b / e
 (mv)  m0 v  mv0
 ( J )  ( J -tilde( J ))0
 (B )  B0  (B )0
The partial derivatives are grouped together into the linear equation:
 mI
0

 0
0
J
0
 ba 0
0    v  ( Fb )r  m0
0      (nbba )0r

B0    
0

( Fbba )0  mv0
(nbba )0  ( J -tilde( J  ))0
I
( Fbbw B)0    v  ( Fbba )u0 



     ( nbba )u0   u
0
  

( B )0     0 



This completes the aircraft problem.
Pendulum
A pendulum swings about one axis as shown in Figure 4b. The points:
e  origin of inertial frame
b  cm of pendulum
p  rotation point, fixed with respect to b and for convenience set equal to e
The generalized mass, velocity and force are:
M  diag(mI
J )  mass properties (6  6)
q  [veb/ e ; bb / e ]  generalized velocity (6 1)
Q  [ Febw ; nbbe  nbbg ]  generalized force (6 1)
22
The individual forces and torques are:
Febw  mge3  weight (3 1)
nebe  e2 (b   k  )  rotary spring-damper (3 1)
nebe  bb/ e J bb/ e  0  gyroscopic torque, zero because bb/ e is eigenvector of J (31)
The minimal information needed to define the joint is:
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Ab/ e  Ry ( )  direction cosine matrix (3  3)
B  Bbb/ e  e2  axis of rotation (B is shortcut name, 3 1)
rbp / b   e3  pivot point (3 1)
b/ e    Euler angle (11)
Derived terms are:
bb/ e  Bb / e  rotational rate (3 1)
Y  Ybb/ e  e2 Ab/ e [e1 e3 ]  orthoganol axes (Y is shortcut name, 3 1)
b  bbp / b   Ae/b rbp / b  cross product (b is shortcut name, 3  3)
Constraints define the pivot point and the axis of rotation, reducing the problem to one degree of freedom:


rep / e
  constraint function (5 1)

[e1 e3 ]T Ae / b e2 
I b 
q  
 Jacobian (5  6)
T
0 Y 
 bB 
H2  
  null space (6 1)
 B 
q2    unconstrained generalized velocity (11)
The unconstrained dynamic equation is:
( J yy  m 2 )  b   (k   m g sin  )  0
The linear dynamics about the straight down equilibrium is:
( J yy  m 2 )  b   (k  m g )  0
The problem reduces to a Galilean pendulum if the rigid body is reduced to a point mass (J = 0) and the rotary
damping and spring are removed (kβ = 0 and bβ = 0):
  ( g / )  0
This completes the pendulum example.
Four Body Problem
The problem has a main body b with three appendages k1 to k3 attached respectively at joints p1 to p3. A diagram
is in Figure 4c. This problem could be a satellite with three attachments. The problem can also be a rotor hub
attached to a test stand with three blades. The generalized mass, velocity, and force are:
M  diag( mb I ,
23
Jb ,
mki I ,
J ki ,
) (24  24)
k /e
q  [veb / e ; bb / e ;
veki / e ; k i ;
] (24 1)
i
Q  [ Febt ; nbbt  nbbg ;
kt
kg
i
i
Fekit ; nk i  nk i ;
] (24 1)
The attitude and angular velocity of the main body are defined by:
Ab/ e  direction cosine matrix (3  3)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Bb  Bbb/ e  axes of rotation in the b-frame (Bb is shortcut name, 3  3)
b/ e  Euler angles (3 1)
Each joint can have one to three rotational degrees of freedom, for this example say each has two. The joints are
defined by:
Aki /b  direction cosine matrix (3  3)
k /b
Bi  Bk i
i
 axes of rotation in the ki -frame (Bi is shortcut name, 3  2)
ki /b  Euler angles (2 1)
p /b
rb i
p / ki
rk i
i
 joint location in rigid body b (3 1)
 joint location in rigid body ki (3 1)
The cross products that will be needed are:
p /b
 cross product (bb is shortcut name, 3  3)
p / ki
 cross product (bi is shortcut name, 3  3)
bb  bb i
bi  bk i
i
k /b
bb i
 cross product (3  3)
The velocity relationship for each joint is:
k /e
vepi / pi  (veb / e  bbbb / e )  (vepi / e  bik i )  0
i
kki /b
i
  Aki / bbb / e
k /e
 k i
i
 Biki /b
The velocity relationships are grouped into matrix form H a q  Hb q2 , from which:
I
0

I

0

H 2 q2   I

0

I

0
0
0
0
I
0
0
k /b
bb1
b1B1
B1
0
0
Ak1 / b
k /b
0 b2 B2
0
B2
k /b
0
0
bb 2
Ak2 / b
bb 3
Ak3 / b
0
0
0
0 
 vb / e 
0  e 

b/e
0   b 

 
 k /b 
0  1 
 
0   k2 / b 
 k3 / b 
b3 B3 
 (121)
B3 
(2412)
The lines are included to help show the pattern of terms. This completes the four body problem.
24
Connected Joints
The problem has a main body b connected successively to bodies h and k, as shown in Figure 4d. This problem
could be a rotorcraft fuselage, hub, and single blade. It also could be an upper arm, forearm, and hand. The
generalized mass, velocity, and force are:
M  diag( mb I ,
q  [veb/ e ;
Jb , mh I ,
bb/ e ;
veh / e ;
J k ) (18 18)
J h , mk I ,
hh / e ;
vek / e ;
kk / e ] (18 1)
Q  [ Febt ; nbbt  nbbg ; Feht ; nhht  nhhg ; Fekt ; nkkt  nkkg ] (18 1)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
The attitude and angular velocity of body b are defined by:
Ab/ e  direction cosine matrix (3  3)
Bb  Bbb/ e  axes of rotation in the b-frame (Bb is shortcut name, 3  3)
b/ e  Euler angles (3 1)
The connecting joint and body h are defined by:
Ah/b  direction cosine matrix (3  3)
Bh  Bhh/b  axes of rotation in the h-frame (Bh is shortcut name, 3  Lh )
h/b  Euler angles (Lh 1)
rbs /b , rhs / h  joint location in bodies b, h (each 3 1)
The connecting joint and body k are defined by:
Ak / h  direction cosine matrix (3  3)
Bk  Bkk / h  axes of rotation in the k -frame (Bk is shortcut name, 3  Lk )
k / h  Euler angles (Lk 1)
rhp / h , rkp / k  joint location in bodies h, k (each 3 1)
The velocity relationships for each joint are grouped into matrix form H a q  Hb q2 , from which:
 Ae / b
 0

A
 e/b
H 2 q2   0

 Ae / b
 0

0
I
0
0
bbh / b
Ah / b
bhh / s Bh
Bh
bbk / b
Ak / b
bhk / s Bh
Ak / h Bh
0
 b/e 
0   ve 
 b/e 
0   b 

0   k1/ b 
 k 2/ b 

bkk / p Bk  
k 3/ b 

Bk 
(18(6 Lh  Lk ))
((6  Lh  Lk )1)
Note that the Ae/b term is used in this example to convert the velocity of body b to the b-frame. This completes the
connected joint problem.
Rotorcraft problem:
The rotorcraft problem in Figure 4e is a combination of the above multi-appendage and connected-joint
problems. A minor simplification is that points s and h are equal (i.e., the hub connection point is at the cm of the
hub). Each of the n blades has L = 1 dof if just flap is modeled, or L = 2 if both flap and lag are modeled. It follows
immediately that:
25
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
0
0
0
 Ae / b
 0
I
0
0

h
/
b
A
bb
0
0
 e/b
Bh
0
 0 Ah / b

k /b
k /h
k / p1
bh1 Bh bk1 Bk1
 Ae / b bb1
1
H 2 q2  
Bk1
 0 Ak1 / b Ak1 / h Bh



kn / b
k /h
bh n Bh
0
 Ae / b bb

0
 0 Akn / b Akn / h Bh
0
0 
b/e
0  vb 


0  bb / e 


0  h / b 
 

0   k1 / b 




 kn / b 
kn / pn
bk
Bk 3 
n
 ((7  nL )1)
Bk 3 
((12 6 n )(7  nL ))
This completes the rotorcraft problem.
Wind turbine problem
The wind turbine in Figure 4f has three blades facing into the wind. The rigid bodies and connecting points are:
e  origin of inertial frame
t  cm of tower (2 rotary dof)
n  cm of nacelle (1 rotary dof)
h  cm of rotating hub (1 rotary dof)
g , s  connection tower to ground, nacelle
k1, k2 , k3  cm of blades (each has just flap dof )
p1, p2 , p3  connection of blades to rotating hub
For this problem the velocity relationship H a q  Hb q2 is listed below using just one blade:
 I
btg / t
0
0 0
0

I
0
0 0
0
 0

p /t
p/n
I
bn
0
0
  I bt
 0 A
0
I
0
0
n/t

 0
0  I bhh / n
I
0

0
0
0

A
0
I

h/ n

0 0
0  I bhp / h
 0
 0
0 0
0 0  Ak / h

0
0
0
0
0
0
I
0
0

0

0
0

0

0

bkp / k 
I 
(2424)
26
 vet / e 

  0
0
0
0
 tt / e  
0
0
0 

  Bt
n/e
 ve   0
0
0
0
 n/e  

0
0
n   0 Bn

 h/e   0
0
0
0
 ve  

 h/e   0
0 Bh
0

 h  
0
0
0
 k /e   0

 ve   0
0
0 Bk 
 k / e 
 k 
(245)
(241)
 t / e 
 
 n /t 
h / n 


k / h 
(51)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
From which it follows that:
 bt / g B
t
 t
B

t
 n/ g
 bt Bt
 A B
 n/t t
 bh / g B
t
 t
 Ah / t Bt
 k1 / g
b
Bt
H 2 q2   t
 Ak1 / t Bt
 k /g
bt 2 Bt

 Ak2 / t Bt

 btk3 / g Bt

 Ak / t Bt
 3
0
0
0
0
0
0
0
0
bnn / p Bn
Bn
0
0
0
0
0
0
bnh / p Bn
Ah / n Bn
0
Bh
0
0
0
0
bnk1 / p Bn
bh1 Bh
bk1 Bk1
0
Ak1 / n Bn
Ak1 / h Bh
Bk1
0
bnk2 / p Bn
bh 2 Bh
0 bk 2
Bk2
Ak2 / n Bn
Ak2 / h Bh
0
Bk2
bnk3 / p Bn
bh 3 Bh
0
0
Ak3 / n Bn
Ak3 / h Bh
0
0
k /h
k /h
k /h
k /p
1
k /p
2
0

0

0
0
 
t /e 
0 

 n / t 
0 


 h / n 
0 
 k /h 
 1
0  k / h 
 2 


0  k3 / h 

(71)
0

k /p
bk 3 Bk3 
3

Bk3 
(367)
This completes the wind turbine problem.
VII.
Conclusions
Two methods are presented for deriving reduced order, unconstrained multi-body dynamics. A closed form
description is given of the transformation matrix H2, which is used to create the reduced order problem. This
transformation matrix can be used as a starting point either for numerical simulation or for symbolic expansion of
the dynamic equations. Symbolic expansions can be used, for example, to track variables and parameters in the
dynamic equations, or for deriving linear equations.
The model reduction methods are intended to be used for problems with a modest number of rigid bodies, less
than a dozen or so. The examples presented here have one to six rigid bodies, but even this is enough for a wide
array of aerospace problems. Starting with six dof dynamics for each rigid body and then whittling them down using
constraints is a good way to make sure that non-intuitive gyroscopic and inertia coupling terms are properly
included, as opposed to the more usual approach, which is not to include these terms and hope they do not matter
(and fortunately, usually, they do not).
The reduced order methods have been developed and used by the author on various projects involving human
bio-dynamics, rotorcraft, tilt-rotor, and wind turbine problems. Matlab-coded versions of both numerical and
symbolic methods have been developed for example problems, but a general purpose utility is not yet available.
Appendix: Quaternions
The orientation of a frame can be defined using quaternions. The advantage relative to using Euler angles is the
update equation does not have a gimbal lock singularity condition. This disadvantage is the transformation from
quaternions to Euler angles is complicated. A multi-body problem can be formulated entirely using quaternions,
thereby avoiding altogether the use of any Euler angles.
There are good reasons to use both Euler angles and quaternions. The aviation community has stubbornly held
on to Euler angles because for manned flight the physical intuition of Euler angles is important, and as a result many
handling quality requirements are written using Euler angles. The gimbal lock problem does not occur for rotations
constrained to one or two axes, and for these constrained transformations Euler angles can be used without any
problem with singularities.
A summary of quaternion vectors is presented in this Appendix, mainly to introduce notation consistent to that
used in the main body of this work. For background the work of Phillips 16 is recommended. The matrix notation
used here is based on the work of Nikravesh.17
27
Euler’s Theorem and Definitions
The underlying theorem that justifies the use of quaternion vectors is Euler’s theorem 17, which states “The
general displacement of a body with one point fixed is a rotation about some axis.” Using this theorem, the direction
cosine matrix for any 3-axis transformation between frames is equal to a rotation about a single axis:
Ab/ a  Ru ( )  uuT  u sin   uu cos  (3  3)
Where:
  right-handed rotation angle (11)
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
u  unit-norm axis of rotation (3 1)
For Euler angle rotations the axes of rotation are fixed in some frame. This is not the case for the quaternion vector.
There is no frame in which the derivative of u is zero.
The Euler-Rodrigues formulation16 is the following 41 vector called a quaternion:
 cos   e 
2   0 (4  1)
pb / a  
 
u sin 2   e 
Where:
e0  cos 2  Euler-parameter (11)
e  [ex ey ez ]T  u sin 2  Euler parameters (3 1)
Following the work of Nikravesh17 it is convenient to place the Euler parameters into the following matrices:
Gb / a  [e eo I  e] (3  4)
Lb / a  [e eo I  e] (3  4)
T

 
T
Pb / a  eo e
 =  pb / a Gb / a  (4  4)
e
e
I

e
0


The subscripts are added here to be consistent with the notation used in the main body of this work. Another type
of notation is used to emphasize the axis and angle of rotation:
pu ( )  quaternion vector with right-handed rotation  about unit-norm axis u (4 1)
Pu ( )  tranformation matrix created with u and  (4  4)
Properties
Transpose properties: The quaternion vector and related matrices obey the following transpose properties:
uT u  1 (1 1)
GT G  LL  I  ppT (4  4)
pT p  1 (1 1)
PT P  PPT  I (4  4)
Gp  Lp  0 (3  1)
U T U  UU T  I (4  4)
GGT  LLT  I (3  3)
U 0T U 0  I (2  2)
Derivatives: Derivative properties include:
pT p  0 (11)
Gp  Lp  0 (3 1)
LGT  LGT (3  3)
Ab / a 
d
dt
tilde(Gp)  GGT  GGT (3 1)
tilde( Lp)  LLT   LLT (3 1)
( LGT )  LGT  LGT  2LGT (3  3)
bb/ a   Ab/ a AbT/ a  (2LGT )(GLT )  2LLT  2tilde( Lp) (3  3)
28
bb/ a  2Lp (3 1)
(25)
Non-uniqueness: The quaternion vector p created by a right-handed rotation by  about u is the same as the
vector p created by a right-handed rotation by – about –u, in other words:
pu ( )  pu (   ) (4 1)
The quaternion vector with the opposite sign maps to the same direction cosine matrix. Changing the sign
corresponds to adding 2 to the rotation:
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
pu ( )   p(  2 ) (4 1)
Quaternion to direction cosine matrix: Given pb/a the corresponding direction cosine matrix is:
Ab/ a  Ru ( )  Lb/ aGbT/ a  I  2e(e0 I  e) (3  3)
(26)
The expanded version is:
 e2  e2  1 e e  e e e e  e e 
x 2
x y
0 z
x z
0 y
 0

2
2
1
Ab / a  2 ex e y  e0 ez e0  e y  2 e y ez  e0ex  (3  3)


2
2
ex ez  e0 e y e y ez  e0ex e0  ez  12 
(27)
Simulation: The update equation for the quaternion is:
pb / a  12 LTb / abb / a (4 1)
(28)
At each step of the simulation, the differential equation is used to update the quaternion vector and then the
direction cosine matrix is created for the new orientation. To derive Equation (28) multiply Equation (25) on the left
by LT/2.
The quaternion vector pb/a has unit norm, and the update equation does not guarantee that the unit norm property
is maintained. This problem is studied in Phillips16 who concludes that for implementations using double precision
algebra and fourth order Runge-Kutta integration the numerical errors are not significant for “short” simulations less
than a few million iterations. Nevertheless the quaternion vector can be re-balanced at each step or to save time
periodically using:
( pb / a )new 
pb / a
T
pb / a pb / a
 pb / a

3
2
 12 pbT/ a pb / a

Multiple frame transformations and converting Euler angles to quaternion: Quaternion vectors for multiple
transformations between frames can be combined in a way similar to what is done for direction cosine matrices:
Ad / a  Ad / c Ac /b Ab/ a (3  3)
pd / a  Pd / c Pc /b pb/ a (4 1)
Note how the subscripts multiply. This transformation can be used to convert from Euler angle rotations to the
equivalent quaternion vector, for example:
Ab/ e  Rx ( ) Ry ( ) Rz ( ) (3  3)
pb/ e  Px ( ) Py ( ) pz ( ) (4 1)
Direction cosine matrix to quaternion: From Equation (27) it follows that:
e0   (1  a11  a22  a33 ) / 4
29
[a23  a32 a31  a13 a12  a21 ] / (4e0 )
[e a / e a / e ]
 x 12 x 13 x
eT  
[a12 / e y e y a23 / e y ]
[a13 / ez a23 / ez ez ]

if e0  0
if e0  0 and ex   (1  a11 ) / 2  0
if e0  0 and e y   (1  a22 ) / 2  0
if e0  0 and ez   (1  a33 ) / 2  0
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
Quaternion to Euler angles (special case): Given a quaternion vector pb/e the problem is to compute the yaw,
pitch, and roll Euler angles  for respective rotations about the zyx axes. Compute the direction cosine matrix as
shown Equation (27) and compare this term-by-term with the following version:
c c
s c
 s 

Ab / e  Rx ( ) Ry ( ) Rz ( )    s c  c s s c c  s s s c s  (3  3)
 s s  c s c c s  s s c c c 


(29)
It follows that:
2
2 
  atan2  a13 , a11
 a12



  atan2(a12 , a11 ) if    2
  atan2(a23 , a33 ) if    2
The function atan2 is the four-quadrant version of the arctangent. At the gimbal lock condition  = /2:
0
0
1

Ab / e  Rx ( ) Ry ( 2 ) Rz ( )   sin(  ) cos(  ) 0  (3  3)
cos(  )  sin(  ) 0 


   atan2(a21, a22 ) if    2
The difference  –  can be determined but not the individual roll and yaw angles. There is no ambiguity with
either the direction cosine matrix or the quaternion, just the Euler angles.
References
1
Takahashi, M. D., “A Flight-Dynamic Helicopter Mathematical Model with a Single Flap-Lag-Torsion Main
Rotor,” NASA TM 102267, USAAVSCOM TM 90-A-004, 1990.
2
Padfield, Gareth D., “Helicopter Flight Dynamics: The Theory and Applications of Flying Qualities and
Simulation Modeling,” AIAA Education Series, Reston, VA, 1999.
3
Thompson, Peter M., "Rotorcraft Dynamics: Derivation And Literal Expansions " Systems Technology, Inc.,
Hawthorne, CA, Working Paper STI-WP-1380-04, June 2008.
4
Ferguson, Samuel W, “Development and Validation of a Simulation for a Generic Tilt-Rotor Aircraft,” NASA
CR-166537, April 1989.
5
Klein, Gary D., “Linear Modeling of Tiltrotor Aircraft (In Helicopter and Airplane Modes) for Stability
Analysis and Preliminary Design,” Naval Postgraduate School, Monterey, CA, Jun 1996.
6
Kleinhesselink, Kristi M., “Stability and Control Modeling of Tiltrotor Aircraft,” Masters Thesis, Univ. of
Maryland, 2007.
7
Cooper, Jared, David. G. Ward, Peter. M. Thompson, “Data-Driven Simulation-Updates of Tiltrotor Aircraft,”
American Helicopter Society 65th Annual Forum, Grapevine, TX, 2009.
8
Jonkman, J. M., “Modeling of the UAE Wind Turbine for Refinement of FAST_AD,” National Renewable
Energy Laboratory, NREL/TP-500-34755, 2003.
9
Lee, Donghoon. and Dewey. H. Hodges, “Multi-Flexible-Body Analysis for Application to Wind Turbine
Control Design,” National Renewable Energy Laboratory, NREL/SR-500-35228, 2004.
10
Jain, Abhinandan, “Unified Formulation of Dynamics for Serial Rigid Multibody Systems,” J. GCD, vol. 14,
Man-June 1991, pp. 531-542.
30
Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 26, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2011-6277
11
Stoneking, Eric, “Newton-Euler Dynamic Equations of Motion for a Multi-body Spacecraft,” AIAA (rest to be
determined).
12
Thompson, P. M. C.-Y. Liang, D. H. Klyde, and R. W. Allen, "Combined Terrain, Vehicle, and Digital Human
Models Used for Human Operator Performance Analysis,” SAE 2004-01-2152, 2004.
13
Haug, Edward J., “Computer Aided Kinematics and Dynamics of Mechanical Systems Volume I: Basic
Methods,” Allyn and Bacon, Boston, 1989.
14
Kane, Thomas R. and David A. Levinson, “Dynamics: Theory and Applications,” McGraw-Hill Book
Company, New York, 1985.
15
Britting, Kenneth R., “Inertial Navigation Systems Analysis,” Wiley-Interscience, New York, 1971.
16
Phillips, W. F., Hailey, C. E., and Gebert, G. A., “A Review of Attitude Representations Used for Aircraft
Kinematics,” Journal of Aircraft, Vol. 38, No. 4, 2001, pp. 718–737.
17
Nikravesh, Parviz E., “Computer-Aided Analysis of Mechanical Systems,” Prentice-Hall, Englewood Cliffs,
New Jersey, 1988.
31
Документ
Категория
Без категории
Просмотров
0
Размер файла
888 Кб
Теги
6277, 2011
1/--страниц
Пожаловаться на содержимое документа