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 a1Hb and H 2 H a1 ( 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 (31) 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 (31) 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 (31) 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 31) xb , yb , zb unit vectors at point b respectively along x, y, and z-axes (each 31) f xa , g xa , hxa vectors at point a (unit norm, usually but not necessarily orthogonal, each 31) 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 Jbbb / 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 (31) 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 a0 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 (31) aea /b vea /b rea /b the letter a is reserved for acclerations computed in the inertial frame (31) Aa /b Aa /bbb/ a ab/ a Aa /b derivative of a direction cosine matrix (3 3) Aa /b Aa /b (bb/ a bb/ abb/ a ) (ab/ a ab/ aab/ 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 2bb /e rbb /e bb /e rbb /e bb /ebb /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/ ab/ 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 / ab / a , where: bb / a Bbb / ab / a , where: bb / e Bbb / eb / 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 (11) 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 J111/ 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 /bbb/ 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 mbbb / 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 mbbb / ebb / e rbp / b Febt (6 1) J bp / b bb / e bb / e J bp / bbb / e nbpt (4) Where: vep / e veb / e Ae /bbb / 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 a1Hb H 2 H a1 ( 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 ps rcs / c b b b e c gb gb ps ps 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 (123) 13 q [I (3 1) 0 I b pp / b 0 s/c bc I 0 0 (129) 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 (912) (91) 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 (1112) (111) 0 0 I (1211) [ 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 (121) T hb ] Ab / c fc I 0 0 0 I 0 H1 H 2 0 0 0 B 0 Ac / b (122) 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 (124) 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 (128) Joint: Spherical + Parallel rep / s hb ]T (125) 0 0 0 q [0 Y Ac / b T T 0 B veb / e bb / e q2 vec / e c / b (1012) (101) (1210) 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 (812) (81) 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 (127) (712) (71) 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/bbb/ e rbp / b derivative needed for H 2 (3 3) bcs / c Ae/ ccc / 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 /bbb / e ) (vec / e bcs / ccc / 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 (11) 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 /bbb / e Bcc /bc /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 ] (112) 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 (11) 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 /bbb / e Bcc /bc /b relative rotational rate (3 1) c /b Euler angle (11) 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 fccc /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 a1Hb 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 / ki / 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 ps 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 (11) u d / unit vector along the line (3 1) f ( ) k ( 0 ) b spring-damper force (11) 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 (11) n( ) k ( 0 ) b spring-damper force (11) 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 ] (112) 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 Bc /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 pe 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, 31) B Bbb/ e rotational axes (B is shortcut name, 3 3) m, J mass properties (11 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 (mv) m0 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 m0 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 (31) 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 (11) Derived terms are: bb/ e Bb / 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 (11) 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 bbbb / e ) (vepi / e bik i ) 0 i kki /b i Aki / bbb / e k /e k i i Biki /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 (121) B3 (2412) 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 (2424) 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 (245) (241) t / e n /t h / n k / h (51) 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 (71) 0 k /p bk 3 Bk3 3 Bk3 (367) 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 (11) 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 41 vector called a quaternion: cos e 2 0 (4 1) pb / a u sin 2 e Where: e0 cos 2 Euler-parameter (11) 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 (11) 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 ( ) pu ( ) (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 / abb / 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

1/--страниц