MATHEMATICAL MODELLING OF FLEXIBLE MULTIBODY DYNAMICS WITH APPLICATION T O ORBITING SYSTEMS by Ahmed El-Hady M . Ibrahim M. Eng., McGill University, 1981 A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies Department of Mechanical Engineering We accept this thesis as conforming to the required standard The University of British Columbia April 1988 © Ahmed El-Hady M . Ibrahim, 1988 In presenting degree freely at this the available copying of department publication of in partial fulfilment University of British Columbia, for this or thesis reference thesis by this for his thesis and scholarly or for her Department DE-6(3/8-n Columbia I I further purposes gain the requirements agree that agree may be It is representatives. financial permission. The University of British 1956 Main Mall Vancouver, Canada V6T 1Y3 study. of shall not be that the for Library permission granted by understood allowed an advanced shall for the that without make it extensive head of copying my my or written ABSTRACT A relatively general formulation for the governing equations of motion, applicable to a large class of flexible multibody systems, is developed using a concise matrix format. The model considered consists of a number of arbitrarily connected flexible deployable members forming branched and closed loop configurations. Joints between bodies are permitted up to six degrees of freedom in translation and rotation. To be effective, the matrix-Lagrangian formulation necessitates development of the kinetic energy expression in a quadratic form in terms of the system velocities. The mass matrix associated with such a quadratic form is known for simple systems such as a collection of point masses, a group of connected rigid bodies, and a discretized flexible structure. However, for a multibody system, where the contributing forces arise from system's translation, rotation, elasticity, deployment, and their interactions, such an expression is not available. To fill this gap, multibody kinematics is developed in terms of the elements of the geometry matrix, which uniquely describes the configuration of branched systems. The characteristic dynamical quantities, i.e., elements of the mass matrix, are identified and the formulation is approached in an increasing order of complexity. The concept of specified and generalized coordinates together with established procedures of analytical dynamics lead to characteristic quantities ( Lagrangian, Hamiltonian, etc. ) and finally result in governing equations of motion which are new to the multibody dynamics. To account for flexibility in a consistent manner, a second-degree nonlinear displacement field is permitted. Alternatively, a linear displacement field can be used if the nonlinear terms up to the fourth-degree are preserved in the strain energy. An algorithm for calculating the stiffness matrix of a flexible element is developed, where terms up to the third-degree of nonlinearity in displacement are retained. Application of this versatile formulation is illustrated through a set of examples of contemporary interest. They pertain to a spacecraft comprising of a central rigid body with attached flexible appendages. The configuration corresponds to a large class of present and planned communication satellites. It can also represent the Space Shuttle based deployment of beam and plate type appendages aimed at scientific experiments or construction of the proposed Space Station. The system static equilibrium and stability are discussed. A computer code is developed and specialized to the specific cases in hand. Typical results ii iii of an extensive parametric study are presented for two particular situations : (i) the Space Shuttle based deployment of a beam or a plate type structural member; (ii) the configuration similar to the Waves In Space Plasma (WISP) experiment jointly proposed by Canada and the U.S.A. The problems are analyzed systematically, through progressive introduction of complexity, to help appreciate interactions between librational dynamics, flexibility, deployment, inertia parameters, orbit eccentricity, initial conditions, appendage orientation, etc. The information is fundamental to the missions concerned and essential to help develop appropriate control strategies. CONTENTS ABSTRACT ii LIST OF TABLES viii LIST OF FIGURES ix NOMENCLATURE xvi ACNOWLEGMENT xxiv PART I : MATHEMATICAL MODELLING 1 INTRODUCTION 1.1 Preliminary Remarks 2 1.2 A Comment on the General Dynamical Model for Multibody Systems . . . 4 1.3 Inertial Forces : Present State of the Problem 5 1.4 A Brief Review of the Literature 5 1.5 Outline of the Present Investigation 2 3 17 KINEMATICS OF A MULTIBODY SYSTEM 2.1 Modelling of a Multibody System 21 2.2 Kinematics : Single Body with Local Coordinates 25 2.3 Matrix of Geometry 29 2.4 Kinematics of a Multibody System 30 2.4.1 Position vectors 32 2.4.2 Velocity vectors 37 2.5 Kinematics of Flexible Multibody System 38 2.6 Inertial Versus Relative Coordinates 45 KINETIC PARAMETERS 3.1 Kinetics of a Single Body 48 3.2 Center of Mass 53 iv V 3.2.1 Center of mass of a body in a multiframe system 53 3.2.2 Center of mass of a multibody system 54 3.2.3 Velocity of the center of mass 56 3.3 Kinetic Energy of a'Body in a Multiframe System 59 3.3.1 General form 59 3.3.2 Translational component of kinetic energy 61 3.3.3 Coupled translational-rotational component of kinetic energy . . . . 62 3.3.4 Rotational kinetic energy 3.4 Kinetic Energy of a Multibody System 65 3.5 Effect of Hinge Flexibility 69 3.6 Holonomic Constraints 73 3.7 Potential Energy 74 3.7.1 4 63 Gravitational potential energy 75 3.7.2 Nonlinear strain 79 3.7.3 The strain energy and stiffness matrix 81 EQUATIONS OF MOTION 4.1 Preliminary Remarks 84 4.2 Functions of Interest 87 4.2.1 Kinetic energy 87 4.2.2 Lagrangian 87 4„2.3 Generalized momenta 88 4.2.4 Hamiltonian 88 4.2.5 Integrals of motion 89 4.3 Lagrange Equations For a Multibody System 89 4.3.1 Structure of the linearized equations 93 4.4 Equations of Motion in Terms of the Hamiltonian 95 P A R T II : A P P L I C A T I O N TO O R B I T I N G S Y S T E M S 5 A T T I T U D E DYNAMICS OF GRAVITY GRADIENT F L E X I B L E SPACECRAFT 5.1 Spacecraft Model 97 5.2 Kinematics 100 5.3 Kinetics 103 5.4 Potential Energy 108 5.5 Equations of Motion 110 vi 5.6 Static Equilibrium and Stability of the Model 6 Ill RESULTS AND DISCUSSION 6.1 Remarks on Computational Approach 128 6.1.1 Computer program 128 6.1.2 Computational considerations 131 6.2 The Rigid Space Shuttle Dynamics 135 6.3 The Orbiter Based Beam Response 141 6.4 Dynamics of the WISP Type Configuration 167 6.4.1 System treated as rigid 6.4.2 System treated as 171 flexible 6.5 The Orbiter Based Plate Dynamics 7 185 206 CONCLUDING REMARKS 7.1 Multibody Mathematical Model 222 7.2 Examples Illustrating Application of the Formulation 224 7.3 Recommendations for Future Study 227 BIBLIOGRAPHY 229 APPENDIX A : DERIVATION OF KINETIC PARAMETERS A . l Displacement Field 235 A.2 First Moment and the Center of Mass 237 A.3 Linear Momentum 238 A.4 Moment of Momentum 239 A.5 Kinetic Energy 241 A.6 Mass Moment of Inertia 244 APPENDIX B : RELATIONS IN NONLINEAR ELASTICITY 253 APPENDIX C : DIFFERENTIATION OF T H E MASS MATRIX Cl Partial Differentiation of the Mass Matrix 256 C.2 Components of the Mass Matrix 258 C.3 Partial Differentiation w.r.t. Mass 260 C.4 Partial Differentiation w.r.t. the Deployment Coordinate 260 C.5 Partial Differentiation w.r.t. the Translational Coordinate 262 C.6 Partial Differentiation w.r.t. the Deflection Coordinate 263 vii C. 7 Partial Differentiation w.r.t. the Rotation Parameter APPENDIX D : D. l 264 CHARACTERISTICS OF T H E B E A M MODES Fixed-Free Beam D.2 Free-Free Beam 266 267 APPENDIX E : INTEGRALS ASSOCIATED WITH FLEXIBILITY 270 APPENDIX F : BASIC KINETIC QUANTITIES 273 APPENDIX G : MODIFIED EULERIAN ROTATIONS 275 LIST O F TABLES Table D - l : The first four eigenvalues and associated constants for a fixed266 free beam Table D-2 : Table D-3 : The first four eignvalues and associated constants for a free-free beam 268 Some numerically evaluated integrals 269 viii LIST O F FIGURES 1- 1 Newton-Euler approach to rigid multibody sytems : (a)free body diagram for body Bk ; (b) augmented body Bk with its barycenter 2- 1 8 Multibody domains and subdomains. Note the typical direct path to a point mass in the system depending on its location 2-2 22 A cut introduced to the closed loop to define unique paths to points in a system 23 2-3 Multibody coordinates . . 2-4 Position vector, in local coordinates, to a point mass on a flexible body with i 24 deployment 26 2-5 The matrix of system geometry 31 2-6 Position vectors in a multibody system 33 2-7 Intermediate coordinates for hinge translational motion 40 2- 8 Illustrations of discretized system topology leading to : (a)Lagrangian formulation ; (b) Newtonian formulation (inertial reference point) ; (c) Newtonian formulation (reference point is system cm.) ; (d)Lagrange-Newton formulation 46 relative to frame L A . . . . 3- 1 Motion of mass element in generic body 3-2 Center of mass of a generic body in a multiframe system 55 3-3 Translational motion of the system coordinates relative to an intermediate frame 71 5-1 Geometry of orbiting spacecraft with beam and plate type appendages . 98 5-2 Reference coordinate systems for a spacecraft with beam and plate type appendages 5-3 49 99 Vibration modes forfixed-freeand free-free beams : (a)the first four modes of afixed-freebeam ; (b)the first two modes of a free-free beam ix 101 X 5-4 Modified Eulerian rotations showing yaw(/3), roll(7), and pitch(a) librations. 107 5- 5 Geometry of a beam in orbit for static equilibrium and stability study : (a) beam in the orbital plane ; (b)beam in the plane normal to the orbit ; (c)beam in the plane tangential to the orbit 122 6- 1 A flow chart showing an approach to system simulation 129 6-2 A schematic diagram showing input to and output of the computer program. 130 6-3 Schematic diagrams of the configurations showing similarity with the cases studied : (a) Space Shuttle ; (b) Orbiter based construction of beam as proposed by Grumman (c) WISP configuration jointly under study by the U.S.A. and Canada ; (d) Solar array experiment conducted in 1984 6-4 134 A comparison between linear and nonlinear responses to a small disturbance with the Orbiter in three different flight configurations 6-5 136 Plots showing inadequacy of the linear analysis to accurately predict instability of the Orbiter 6-6 137 Librational response of the Orbiter to an independent excitation in pitch, yaw and roll. Note the pitch and yaw disturbances lead to essentially uncoupled motions. The system appears to become unstable in yaw through its coupling with roll 6-7 139 Librational response of the Orbiter in the Lagrange configuration with the time history of roll control actually used in practice 6-8 140 Effect of smoothing the peaks and frequency of the roll control on the Orbiter's librational response. Note, the deadband limits can be relaxed appreciably with the stability assured 6-9 142 A schematic diagram showing the orbital coordinates X ,Y ,Z 0 0 0 and beam coordinates x,y,z with their origins at the instantaneous center of mass. In general the two coordinate systems are not coincident 6-10 Response of the beam, located in the orbital plane, to a tip excitation in the first-mode of 4% of its length. Note the transverse oscillations 6X and 6y are 144 xi uncoupled 145 6-11 Plots showing the deflected equilibrium position of the beam acting as a small initial disturbance, in this case in the orbital plane. Due to an initial tip disturbance in the out-of-plane, the beam response is coupled as shown by the 8 -6 plot x 147 y 6-12 Response of the beam to a tip disturbance when located out of the orbital plane. The transverse motions S and 8 are coupled and the vibration plane x y of the beam tip precesses at a uniform rate 6-13 Response plots showing the effect of a change in the beamflexuralrigidity. 149 150 6-14 Effect of the orbital eccentricity on vibratory response of the tip for two orientations of the beam. Note the frequency condensation effect is particularly pronounced near the apogee 151 6-15 Response plots showing the frequency condensation effect in the first two modes due to an orbital eccentricity of 0.2 152 6-16 Effect of the time of deployment on the tip response of the beam located in the orbital plane. Note the beam-tip drift in the orbital plane due to the Coriolis force. In general, the steady state amplitude of oscillations are higher for the deployment case 154 6-17 Effect of deployment strategies on tip response of a beam located normal to the orbital plane. Note a reduction in beam frequency during deployment. The steady state amplitude is essentially independent of the strategy for a given time of deployment 155 6-18 A comparison between responses of a deploying and deployed beam, initially aligned with the orbit normal, to a tip disturbance 6-19 Representative controlled motion of the Orbiter during a typical orbit . 156 158 6-20 Forced response of a beam to the Shuttle librational control. Note the beam frequency is close to that of the roll disturbance 159 6-21 Forced response of a beam extended along the local vertical from the Orbiter in the Lagrange configuration. Note, the beam response closely follows the yaw and pitch excitations with modulations at the natural frequencies in different modes superimposed 160 xii 6-22 Effect of an orbital eccentricity on the forced beam response of the case considered in Figure 6-21 (beam along the local horizontal). The inplane tip response 6 ^ y increases by an order of magintude even in the presence of relatively small smooth disturbances. The frequency condensation effect is still present 161 6-23 Forced response of a beam, aligned with the orbit normal. Note the effect of the dominant roll excitation at 13 cycles per orbit which is particularly clear in the second mode 162 6-24 A comparison between the forced response of a beam, located either along the local horizontal or the orbit normal, with the Shuttle in an orbit of e = 0.2. 163 6-25 A comparison between the forced response of a beam, aligned with the local horizontal, showing the effect of deployment 164 6-26 Effect of deployment on the forced response of a beam located along the orbit normal 165 6-27 A typical case of beam response showing the effect of mode truncation. It appears that the higher modes have very little effect on the time history of the beam deflection : (a) first mode ; (b) first two modes ; 168 (c) first four modes 169 6-28 Plots showing the influence of a permissible error during numerical integration. The allowable error of 1 0 -8 or less resulted in reliable data . . . . 170 6-29 A comparison showing the effect of orbital eccentricity on the libration response of the Orbiter in the Lagrange configuration with an initial pitch disturbance 172 6-30 A comparison showing the effect of orbital eccentricity on the libration response of WISP system treated as rigid in the Lagrange configuration with an initial pitch disturbance 173 6-31 Effect of deployment velocity on the librational response of the WISP configuration, treated as rigid and in the Lagrangian orientation, with an initial pitch disturbance. Duration of deployment : (a) 0.4 orbit ; 176 xiii (b) 0.8 orbit ; 177 (c) 1.5 orbit 178 6-32 Effect of the initial deployed mass on librational response of the WISP configuration with an identical deployment time history and duration : (a) t(0) = 0 ; (b) £(0) = 100 179 180 m 6-33 Influence of the initial deployed length on librational response of the WISP configuration. Note, the deployment time history and duration are identical in the two cases : (a) deployment from zero to 50 meters ; 181 (b) deployment from 100 to 150 meters 182 6-34 Effect of sinusoidal deployment history on the librational response of the WISP configuration : (a) a(0) — 12°. Note, a comparison with Figure 6-32a shows that the time history of deployment has virtually no effect on pitch response ; . . . . (b) a(0) = - 1 2 ° 183 184 6-35 Libration-vibration response of the WISP system, in the Lagrangian configuration, to a disturbance in the form of inplane symmetric deflection in the first beam mode : (a) e=0 ; 186 (b) e=0.2 187 6-36 Libration-vibration response of the WISP system to an inplane asymmetric deflection in the first beam mode : (a) libration and beam tip response ; 189 (b) beam-mode response 190 6-37 Effect of an orbit eccentricity on the libration-vibration response of the WISP system in the Lagrangian configuration. The system is subjected to an initial inplane asymmetric deflection in the beam's first mode : (a) libration and beam tip response ; 191 (b) beam modal response 192 6-38 Vibratory response of the WISP system in the Lagrangian configuration to xiv an out-of-plane symmetric excitation in the beam's first mode. The Shuttle remained unexcited in the librational mode (not shown) 193 6-39 Vibration response of the WISP system in the Lagrangian configuration to an out-of-plane symmetric deflection in the beam's second mode. The librational response was found to be zero (not shown) 194 6-40 A comparison showing the effect of orbital eccentricity on the librationvibration response of the WISP system, in the Lagrange configuration, when subjected to an initial symmetric deflection (first mode) in the out-of-plane direction : (a) librational response for the orbital eccentricities of 0.1 and 0.2 ; . . 195 (b) vibrational response for the orbital eccentricity of 0.2 ; 196 (c) vibrational response for the orbital eccentricity of 0.1 197 6-41 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an impulse along the orbit normal : (a) libration and beam tip response ; 199 (b) beam modal response 200 6-42 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an intial pitch disturbance < (a) libration and beam tip response, e=0 ; 201 (b) beam modal response, e=0 ; 202 (c) libration and beam tip response, e=0.2 203 6-43 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an initial yaw disturbance : (a) libration time history ; 204 (b) beam tip response 205 6-44 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an intial roll disturbance of 3° : (a) libration time history ; 207 (b) beam tip response 208 6-45 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an initial roll disturbance of 9° : XV (a) libration time history. Note the instability in yaw through a coupling with roll ; 209 (b) beam tip response 210 6-46 Effect of the Orbiter'sflightconfiguration on the transient librational response during deployment of a flexible plate-type appendage along the maximum moment of inertia axis 211 6-47 Typical vibrational response of the flexible solar panel during and after the deployment for three different flight configurations of the Orbiter . . . . 213 6-48 The Orbiter-appendage coupled dynamics during deployed and transient deployment conditions of the plate out of the plane of symmetry 214 6-49 The Orbiter-appendage coupled dynamics during deployed and transient deployment conditions of a plate in the plane of symmetry 215 6-50 Vibratory response of a plate to deployment or specified motion during the Orbiter's flight in the nominal configuration 217 6-51 Influence of the deployment rates on the vibratory response of a plate oriented tangent to the orbit 218 6-52 Response of the flexible plate-type appendage as affected by its inplane and out of the plane of symmetry orientations during the Orbiter's flight in the nominal configuration 219 6-53 Vibratory response of the deploying plate as affected by its orientation during the Orbiter'sflightin the Lagrange configuration 220 NOMENCLATURE Roman Symbols a~i position vector to origin O; of the body frame Li in frame L t - i , Fig. 2-6 ya,- position vector to origin O,- of the body frame Li in frame Ly; Fig. 2-6 a* position vector to a point of body j \%i projection of position vector ya,- on local frame Lk k a? projection of position vector a* on local frame Lk 6* position vector to center of mass of flexible body Bi in frame LV, Fig. 3-2 J in frame Ly; Fig. 2-6 J 3—* 6 position vector to center of mass of flexible body L\- in frame Ly; Fig. 3-2 k— . . - jbi projection of position vector y6,- on local frame Lk B{ body number i; also referred to as domain Di c position vector, in frame L,-, to the center of mass of body Bi in its unde- 1 flected configuration iCjk, iVjk ; j, k — 1,2,3 dimensionless vector (lengthrag,)integrals due to the time dependence of the shape functions during deployment and the associated convective terms, respectively. They represent the linear terms in deflection of the mass vector associated with coupled deployment-translational motion of flexible body Eqs. (A.46) and (A.4c), respectively C instantaneous center of mass of the system Cj matrix of direction cosines defining rotation of local frame Ly with respect to inertial frame L J V matrix of direction cosines defining rotation of frame Ly with respect to frame Lk iCjl > n iVjk, i n iN™, dimensionless rag; x nsi matrix integrals. The first two arise because of the dependence of the shape functions on time during deployment, while xvi xvii j,k,l,m,n the third is due to the associated convective term. They represent the = 1,2,3 linear terms, in deflection, of the mass matrix associated with coupled deployment-vibrational motion of body Bi\ Eqs.(A4/), (AAt) and (AAg), respectively di width of plate type appendage d{ unit vector in the direction of deployment of body B{ in frame Li id{nn dimensionless integral defined by dmi differential mass element in body B{ D{ domain or subdomain number t; also referred to as body B{ iD™, iEjl" 1 Eq.(AAk) ; dimensionless ngi x n$- matrix integrals due to time dependence of the t j,k,m,n shape functions during deployment and the associated convective terms, = 1,2,3 respectively. They define the linear terms, in deflection, in the mass vector associated with coupled deployment-rotational motion of body i?,-; Eqs.(AAe) and (AAh), respectively ifjkt dimensionless vector integrals, of length ng,-, defining the linear component j, k = 1,2,3 of deflection of the center of mass of a flexible body -B,-; Fl matrix defining the velocity of a point on body f?,- due to its vibration; Eq.(AAa) Eq.(2.2) Fk matrix defining the velocity of Ok, origin of frame Lk, due to vibration of body Bj, j = k - 1; Eq.(2.18) iFjk 1 ; dimensionless matrix integrals, each of dimension nsi x ns%, defining the j,k,m,n mass matrices of body B{ associated with vibration, rotation and their = 1,2,3 coupling; g*e 3 three dimensional vector defining the first mass moment associated with Eq.(AAd) coupled deployment-rotational motion. It is due to the deployment of body B{ and its rotation about hinge at point Oy; Eq.(3.19c) ig 3mn j,m,n ; dimensionless vector integral, of length n$i, defining the constant term in = 1,2,3 the first mass moment Pg, which is due to coupled vibrational-rotational motion of body B{. It also represents the linear term in the second mass XV111 moment of inertia {J^jk] Eq.(A.4j) matrix defining the angular velocity due to elastic deflection of frame Lk, which is fixed to flexible body Bj, j = k - 1; Eq.(2.22) G/ 1 n$ x 3 rectangular mass matrix associated with coupled vibration-rotation i motion. It arises because of rotation of a vibrating body Bi about point Of, Eq.(3.19d) —X h three dimensional vector defining contribution of the deployment of body e Bi to its linear momentum; Eq.(3.3) H , H} ; matrices defining the linear and nonlinear contributions of the elastic / = 1,2,3 deflection of body Bi to the position of its center of mass, respectively; 1 Eqs.(A.2&) and (A.3a) HI, H\ matrices defining the linear contributions of deflection to h , the linear e momentum of body Bi during its deployment; Eqs.(A.2c) and (A.2d) Hg matrix defining contribution of the vibration of body Bi to its linear momentum; Eq.(3.4) H' 3 k matrix of moment or mixed moment of inertia (j = k or j ^ k) associated with bodies which have two hinges, at Oy and Ok, on their direct path Hf k matrix of moment or mixed moment of inertia [j = k or j ^ k) associated with body Bi with two hinges, at Oy and Ok, along its direct path P moment of inertia matrix for body Bi in frame Li I\ monent of inertia matrix for undeflected body Bi in frame L,- I\, I\ matrices defining linear and nonlinear contributions of deflection, respectively, of body Bi to its mass moment of inertia in body frame L -; t Eqs.(A.2r) and (A.2s), respectively {Ifyjk, [I-sljk matrices defining elements of 1%, J | , respectively, where j,k = 1,2,3 ; Eqs.(A2r) and (A.2s), respectively {J2}y&, [«/|]y& matrices defining {I^jk and [i^y*:, respectively, where j, k = 1,2,3; Eqs.(A.2t) and (A.2v), respectively ii length of Bi in the direction of deployment xix speed of deployment of B{ local frame fixed to body Bi inertial frame of reference an auxiliary frame, fixed to deploying body Bj-i, in which translational motion of coordinate Lj is described; Figure 2-7 matrix projected on frame Lj, defining first mass moments (j = k) and first mixed mass moments (j ^ k) for all bodies with Oj and Ok on their direct paths total mass of the system mass of body Bi nondimensionalized with respect to system mass m rrij m k sum of masses of all bodies on the direct path to Oj nondimensionalized with respect to the system mass m rrij m k sum of masses of all bodies on the direct path to both Oj and Ok nondimensionalized with respect to system mass m number of generalized coordinates associated with deflection of body Bi vectors of direction cosines defining rotation of body coordinates X o , Vo,Zo relative to orbital coordinates xt, yt, zt- n corresponds to rotation about the orbit normal, I about the local vertical, and m about the local horizontal origin of the local coordinate L,vector defining contribution of the deployment of body Bi to its moment of momentum in frame L,-; Eq. (3.4) marices defining linear deflection terms in the moment of momentum vector p^, of body Bi, due to deployment; Eqs.(A2e) and (A.2f), respectively matrices defining nonlinear deflection terms in the moment of momentum XX vector j? , of body B,-, due to deployment; Eqs.(A2g), [A.2h) and (A.3c), e respectively Pg matrix denning the contribution due to vibration of body to its moment of momentum in frame L,-; Eq.(3.4) Pgg matrix defining constant terms in the moment of momentum matrix Pg, of body B{, due to vibration; Eq.(.A.2t) q vector of system generalized coordinates fi position vector from origin ON of the inertial frame LN to origin O,- of local frame L,-; Fig. 2-6 f* position vector, in inertial frame LN, to a point on body Bi; Fig. 2-6 yr,- projection of position vector ya,- on inertial frame LN f projection of position vector o* on inertial frame LN 3 x J r,- velocity of origin O,-, of frame L,-, due to translational motion; Eq.(2.16e) R position vector from the center of Earth to C, the instantaneous center of c mass of the system; Fig. 5-1 Rj matrix defining rotation of frame Ly with respect to Ly s specified velocities S matrix defining geometry of a branched system S k elements of matrix S; equal to unity if frame Ly is on the direct path to body Bk, otherwise zero T system kinetic energy Uj position vector to origin Oy (of frame Ly) in frame L y _ ; Fig. 2-7 U unit matrix matching dimensions of terms involved in the corresponding x equations U g gravitational potential energy Up total potential energy V strain energy , xxi w displacement field associated with body J3,-; Fig. 2-4 w displacement vector associated with Ok, origin of frame Lk, due to defor- 1 k mation of body Bj, j = k — 1; Fig. 2-7 w\ scalar defining the contribution of deployment of body Bi to its kinetic e energy; Eq.(3.4) w^ vector defining contribution due to interactions between deployment and 1 vibration of body Bi to its kinetic energy; Eq.(3.4) W* ,W* matrices defining nonlinear terms in deflection, in the expression for kinetic Wl ,YZ? energy due to deployment of body B,-; Eqs.(A.2k),(A.2l),(A.2m), c ei e and (A.3/), respectively , W* , s matrices defining linear terms in deflection, in the kinetic energy of body Bi due to its deployment-vibration interaction; Eqs.(A.2n), (A.2o) and (A.3g), respectively Wg> W?8ii Wjm matrices defining constant, linear and nonlinear terms in deflection, l,m — 1,• • •,^respectively, in the expression for kinetic energy due to vibration; Eqs.(A2p), (A.Si) and (A.3j), respectively Wg S matrix defining the contribution due to vibration of body Bi to its kinetic energy; Eq.(3.4) Xj,yj,Zj body coordinates, fixed to Bj at point Oj, referred to as frame Lj\ Fig. 2-6 X, Y, Z ,YO,ZQ inertial coordinate system; Fig. 2-6 orbital coordinate system with origin at C; Yo along local vertical, Zo along local horizontal, and XQ aligned with orbit normal; Fig. 5-1 Greek Symbols pitch, yaw, and roll librational angles, respectively; Fig. 5-4 function defining spatial dependence of displacement field w ; Section 2.2 1 0k value of P evaluated at Ok, the origin of frame Lk, j = k — 1; Section 2.5 3 XXII P spatial function defining velocity of a point on body Bi due to deployment x and vibration; Eq.(2.5) Pk value of ft evaluated at Ok, the origin of frame Lk, j = k — 1; Eq.(2.21) Si vector of generalized coordinates associated with deflection of body Bi 8jk kronecker delta; equal to unity when j = k, otherwise zero 6jk (1 - 8jk) 8^n generalized coordinate associated with < £ , 0 3 m n modes of free-free and clamped-free beams, respectively, to represent plate type appendage oscillations 8X[ ,6yi generalized coordinates associated with the /-th mode of vibration of a fixed-free beam; Section 5.2 A ; k = 1,2,3 fc vector of shape functions defining the linear displacement of body 5,- in Pic direction; Appendix A - l A] , / = 1,2,3 vector defining the fc-th column of matrix \I>}, Appendix A - l €jki permutation symbol equal to 1 when j, k, I are in sequence, -1 when not fc in sequence, and zero when any of them repeated f, 77, £ local coordinates associated with appendages; Fig. 5-5 position vector to a point on Bi in the body coordinates x-,y-,2:,- (frame t Li); f + v?; t Fig. 2-4 fj* vector defining velocity of a point on body Bi due to its deployment fj vector defining velocity of Ok, the origin of frame Lk, due to deployment k of body Bj, j = k - 1; Eq.(2.18) 0 true anomaly; 0 = Q {iti}jk vector defining the (j, k) submatrix in $*; Appendix A - l p, universal gravitational constant p\ coordinate associated with body Bi in the dirction of deployment p n-th coordinate of body Bi; n = 1,2,3 corresponding to l n Xi,yi,Zi, respec- xxiii tively position vector to a point on undeformed body Bi in body coordinates p 1 Xi,yi,Zi\ Fig. 2-4 modes of free-free beam <p m $* matrix of admissible or shape functions associated with body Bi <&k $ J evaluated at Ok, j — k — 1 tp modes of fixed-free beam ty], I = 1,2, 3 matrix of admissible functions associated with the nonlinear displacement n field of body B{ uij = W y - 1 fi angular velocity of Lj relative to L y _ i orbital angular velocity Miscellaneous c, s cosine and sine of a function, respectively a small letter (lower case) usually defines a scalar quantity a = {a} bar defines a column vector or an array a skew symmetric matrix formed from vector a A = [A] capital (upper case) letter usually defines a matrix. Exceptions are normally clear such as kinetic energy, T; potential energy, C7t; Hamiltonian, H; Lagrangian, L; etc. { a y } , {a Bja} arrays with j-th element equal to ay and a Bja, (a B),{Ba} array in row and column form, respectively [ a y * ] , [a Bjka] matrices with (j,k) element equal to a y * and a Bjka, [(a)y], [{a}y] matrices with j-th row and column equal to ( a ) y and { a } y , respectively T T T 1 T respectively respectively xxiv A C K N O W L E G M E N T I would like to thank Dr. V. J. Modi for the opportunity he has given me here at UBC and for the enthusiasm and support he has shown. I am also grateful to my colleagues for the exciting academic atmosphere they provided during the course of this work. Without the support of my family back home this work woud not have been possible. Thanks to my wife Sabah and to my daughters Zainub and Roukaya for their encouragement. This research was carried out under a grant from the Natural Sciences and Engineering Research Council of Canada. Additional support was recieved from the University Graduate Fellowship. As an international student, I am much indebted to Canada for her hospitality. XXV To my mother and the memory of my father PART I : MATHEMATICAL 1 MODELLING 1. 1.1 INTRODUCTION Preliminary Remarks The classical dynamicists of the past laid down the principles necessary to describe dynamic behaviour of mechanical systems. However, before the advent of computers, their applications were limited to relatively simple systems. This limitation was accepted with a sense of abandon since, in the absence of analytical tools leading to closed-form solutions, useful information concerning the motion of more complex systems was considered beyond reach. Hence, it is not surprising that, for a long time, the attention remained focused on analytical intricacies of the mathematics of dynamics, while the process of arriving at the equations of motion came to be regarded as a rather routine matter. However, the computer's ability to integrate a system of equations efficiently and with an acceptable degree of accuracy is now established beyond doubt. Increasing attention is, therefore, directed towards analysis of complex systems using realistic mathematical models for machine solution instead of seeking exact or approximate analytical answers for their simplified representations. At the outset it was recognized that geometric complexities, nonautonomous character, coupling, constraints, nonlinearities and hybrid character inherent to the systems made the task enormously challenging. However, resolution of the new generation of problems demanded that the challenge be met. It was apparent that inability to adequately model the problems of contemporary interest can be as great a hinderance now as their analytical solutions was before. There was a renewed interest in classical dynamics with the goal of effectively developing satisfactory mathematical models. Although a great deal of effort was devoted to merely reviving forgotten techniques and classifying existing approaches [1-3], a fresh look at approaching dynamical problems was also taken [2] together 2 3 with development of inductive and symbolic computer aided methodologies [4-9]. These endeavors have been fruitfully applied to nontraditional fields as well as designs of complex systems which perform at higher speeds and demand more stringent limitations on the system weight and working accuracy. Typical examples in dynamics are large space structures, high speed transportation systems, biomechanics, robotics, mechanisms, etc. [10-15]. The design of a complex system is a process of repeated modifications to achieve compromise between competing factors. The initial design is generally based on previous experience, simple analytical models, experimental tests, and numerical analyses. Only seldom is the information available from all the desired sources. Inevitable design changes, imposed by the evolutionary process, lead to iteration cycles aimed at accommodation, checks and balance. Furthermore, new approaches need to be verified and parametric assessment is a must to arrive at a safe regime of operation. Obviously, availability of a working mathematical model for the problem is essential for effectively performing such a comprehensive analysis involving aspects of continuum mechanics, dynamics, stability, control, system identification, etc. It is apparent that availability of a relatively general formulation applicable to a large class of systems would be a distinct advantage here. This has been recognized by the scientific community at large together with industrial corporations and government agencies. Unfortunately, the time-bound character of the industrial contracts and sponsored research have generally limited efforts in that direction. By and large the attitude has been to investigate the specific configuration in hand to an extent considered possible within the constraints of allocated time and money. In the long range, the procedure turns out to be inefficient, lacks perspective and, of course, is devoid of elegance associated with a general approach. University based research, particularly in Canada, may have a useful role to play here. 4 1.2 A Comment on the General Dynamical Model for Multibody Systems A mathematical model for system dynamics is represented by a set of equations viewed by d'Alembert as a balance among forces involved. The forces may be classified into two categories : those which are functions of the system state, its dynamical parameters and possibly time; while the others, in addition, depend on the system geometry, boundary conditions, etc. Forces belonging to the first category are derivable from a single function called the Lagrangian, while those of the second category are obtained through experiments, simplified analytical models or integral forms over boundaries. Examples of the first category are forces due to gravity, elastic deformations and inertia. Forces of the second kind are encountered in the problems of fluid dynamics, machine tool mechanics, etc. This categorization demonstrates our ability or inability to mathematically model forces in nature. Forces of the first type can easily be treated on a general basis and are reducible to various special cases, while forces of the second kind may be accounted for on an individual basis through appropriate input information. The main goal of a multibody dynamical formulation is to determine a general expression for inertial forces (inertia, centrifugal, coriolis, etc.) of the type encountered in complex mechanical systems. These forces are associated with translation, rotation, elastic deformation and deployment, as well as their interactions. A relatively general configuration may be represented by a number of contiguous, flexible structural members (beams, plates, shells, membrane, tether) forming branched or closed loops. The members may be permitted to undergo large relative excursions in length, location, orientation and deflection as would be the case during deployment, slewing maneuvers, etc. Furthermore, the entire system, or a part of it, may be translating and/or rotating in a specified manner. Obviously, a large class of marine, ground-based and aerospace systems would be represented by such a multibody model. 5 1.3 Inertial Forces: Present State of the Problem Historically, efforts aimed at mathematical modelling of dynamical systems have taken two distinct directions : structural vibration against rigid multibody dynamics. More efforts were devoted to vibrational analysis which has become an established discipline as evident from numerous textbooks on the subject with courses in structural dynamics, matrix methods and finite element procedures. Exacting performance requirements, particularly for the aerospace systems, have dictated the necessity of considering elastic deformations in multibody models. Such studies have resulted in computer codes such as DISCOS [16], NBOD2 [17], A L L F L E X [18], TREETOPS [19], and several others, although there is still an on-going debate about many aspects of the problem and reliability of the results. Partly in response to these concerns, a NASA program called the Spacecraft COntrol Laboratory Experiment (SCOLE) is proposed. It will serve as the focus of a design challenge for the purpose of comparing different approaches to modelling, control synthesis, order reduction, state estimation and system identification in conjunction with large flexible structures. It is of interest to recognize that the only available book on the subject [20] covers a single approach to rigid multibody dynamics. In a sense, this suggests a considerable scope for advancement of the present state of the art. 1.4 A Brief Review of the Literature Apparently the advent of the space age together with the computer revolution were required to motivate refinement of this classical discipline. Analyses so far were mostly confined to specific configurations as pointed out before. Numerous contributions in the area of flexible spacecraft dynamics aimed at vibration/libration interaction problems have been reviewed by Modi [21]. A review of modelling techniques for the open and closed loop dynamics of large space systems is presented by Bainum [22]. In the present study, the focus is on the second category of problems classified as multibody dynamics. In spite of the fact that the underlying physics is classical and undisputed, the development of such general models is not straightforward. The difficulity is contributed, in part, by the multiplicity of options open to dynamicists in formulating the problem : (i) The principles of mechanics used in the formulation : Newton-Euler, d'Alembert, Lagranian, Hamiltonian or their variations. For example Kane's approach represents one such variation [23]. (ii) The system geometry considered: chain, cluster, tree or arbitrary configurations with closed loops. (iii) The strategy followed in development of the equations of motion : progressive intoduction of complexity in terms, i.e., number of bodies, flexibility, deployment, etc. (iv) The choice of system variables: absolute, relative, true or quasi coordinates; generalized momenta or velocities; etc. (v) Hinge modelling and permitted motion of the individual body : rigid body translation and rotation, flexible deformation and axial deployment. (vi) Choice of the reference point : center of mass, origin of an arbitrarily located frame, or some point on one of the bodies. (vii) Degree of nonlinearity preserved in the formulation. Obviously, depending on the approach adopted, the resulting equations would vary in form, generality of the model, scope of applicability and computational efficiency. The following review attempts to clarify these points. The two pioneering contributions in the area of multibody dynamics, by Hooker and Margulis [24] and Roberson and Wittenburg [25], were motivated by problems in spacecraft dynamics. They provided the first spark for the subsequent research that has continued to date. These original researchers used the Newton-Euler formulation procedure and limited 7 their scope of work to a system of rigid bodies having a branched configuration. The hinges between the bodies were allowed up to three degrees of freedom in rotation. They used the center of mass of individual bodies and the center of mass of the whole system as reference points. Denoting the mass of body Bk by m , and the associated external and k constraint forces as F be m (p k k k and F , respectively, the translational equation of motion would k + x) = F + F . k k Here, the inertial position of the center of mass of each body is decomposed into two components: one is relative to the overall system center of mass, p , and the second is the position of the system center of mass, x, with respect to the k inertial frame ( Figure 1-1). With m and F representing the system mass and the totality of the external forces, respectively, the sum of the translational equations over all the bodies results in the relation for the system center of mass as, rrix — F. The acceleration of the system center of mass ,x = F/m, can now be eliminated from the translational equation for body B , i.e., k m (p k + F/m) — F k k —F. k While summing the preceeding equations over bodies on either side of a hinge, the constraint forces cancel in pairs at all other hinges leaving the resultant equation with one constraint force as a function of p . k The accelerations p , and in sequence each constraint force, are obtained in terms of the k system geometry and bodies' angular velocities, U . k The rotational equation for the individual body would read as I oJk + w X I - U = k T k + ^2i(rki x F ) + T , where I kl k k k k is the mass moment of inertia of body B k k and the summation is over the number of attached bodies. The external and constraint torques on B k are denoted by T , and T , respectively. The vector r i, from the center of mass k k k of Bk, is normal to the line of action of F . kl The expressions for constraint forces obtained before are now introduced in the rotational equations. F kl Therefore, while the constraint torques appear explicitly in the rotation equations, the forces due to constraints do not. Simplifications of the equations revealed the authors' most celebrated concept of the augmented body and the connection barycenter. The terms in the rotational equations now reflect significance of the barycenter as the mass center of the augmented body, which 8 Figure 1-1 Newton-Euler approach to rigid multibody sytems : (a) free body diagram for body Bk ; (b) augmented body Bk with its barycenter . 9 consists of an individual body augmented at each connection point by a particle having the mass of the corresponding branch set. In most applications, joints between bodies have less than three degrees of freedom. The locked gimbals produce unknown constraint torques in the rotational equations. Hence, the equations are augmented by constraint relations. The latter expresses the fact that the constraint torque at each joint is orthogonal to the difference between the angular velocity vectors associated with the two bodies adjacent to each joint. The well known Lagrange multiplier method was used to obtain a set of governing equations [24,25]. Interestingly, Roberson and Wittenburg [25] used the graph theory to describe the configuration of the bodies. The incidence matrix used in this method can uniquely define the specific topological tree configuration. It also systematically governs the transformation of variables and simplifies the algebra involved. Roberson [26] has further extended this approach to allow for relative translational motion between bodies at the joints. As before, the unknown constraint torques are treated using the Lagrange multiplier method. Computationally there was still a hurdle to overcome. Numerical experiments showed a problem of accumulated error if the number of equations exceeded the number of degrees of freedom. In response, Hooker [27] developed a procedure to reduce the Euler equations to its minimum set. This was achieved by expressing the inertial angular velocities in terms of the relative angular velocities between the system joints. The rotational equations associated with unlocked gimbals were then obtained by summing the rotational equations of all bodies forming a branch-set on one side of the gimbal. Since this process is algebraically involved, particularly when flexibility is included, Ness and Farrenkopf [4] as well as Ho and Gluck [5] used an inductive approach applicable to a flexible n-body system. Ness and Farrenkopf developed a nonlinear rigid body model with flexible terminal connected bodies while Ho and Gluck used a linearized all-flexible configuration. In both the cases, the dynamical equations for each body were first derived analytically. Then they were combined with the equations of the adjacent body through 10 a computerized inductive method. The interacting forces and torques were eliminated at each combining step. Modifications for the case of the locked gimbal axes were performed simultaneously. Ho [28] has pointed out the drawback of this inductive approach. The loading and reloading of the coefficient matrices in the governing equations for the subsystems required considerable computer time. Moreover, the inductive method does not provide overall dynamical equations for the entire system. Hooker [29], using the augmented body approach, developed a mathematical model for a system forming a topological tree configuration with flexible terminal bodies. Likins [30] extended this approach to a set of hinged substructures forming a tree configuration. Each substructure consisted of a central rigid body with attached rigid axisymmetric rotors and nonrigid appendages. The common feature in the above mentioned approaches is the application of the augmented body approach. This involves application of the system center of mass equations to the individual member's translational equations as explained before. Hence, the choice of the system center of mass as the reference point is essential. The procedure is motivated by a desire to simplify translational equations for the whole system which, in many cases, can be uncoupled from the rotational equations. Unfortunately, the simplification is not without cost since the rotational equations are more complicated now (in their coefficients) despite the elegant augmented-body interpretation [31]. Ho [28,32] has used a different approach which he refers to as the direct path method. Here the system reference point is not necessarily its center of mass. A direct path, originating from an assigned home-body, uniquely defined a point on a given body. The reference points were the centers of mass of the individual bodies or hinge locations along the direct path. The translational and rotational equations are now coupled, but the calculation and interpretation of the inertial coefficients are straightforward: they are conventional and mixed inertias of the bodies. Hooker [33] gave up his original idea of the augmented body approach in favour of the direct path method. Once again he considered a tree configuration of rigid bodies 11 with the terminal bodies treated as flexible. Relative translation between the bodies was not allowed. He applied the concept of nested bodies as suggested by Russel [34] together with the approaches presented by Ness and Farrenkopf [4] as well as Ho and Gluck [5], and obtained rotational equations about the attachment points rather than the center of mass. The equations are now free of the constraint forces. The constraint torques were eliminated by projecting the rotational equations along the gimbal axis. The translational equations of the system were obtained using the Newtonian approach in conjunction with a reference point fixed to the main body. The equations associated with the elastic deformations were obtained through the Lagrangian procedure. Ho [28] considered a tree configuration of rigid bodies with the terminal bodies treated as flexible. No relative translational motions were allowed at the hinges. Three different approaches to the problem were compared : the Newtonian, Lagrangian-Newtonian and Lagrangian. In the Lagrangian-Newtonian approach the kinetic energy expression for the individual body was developed and its equation of motion derived using Lagrange's equation with the reactions between bodies treated as external forces. Using a summation and projection procedure the reaction forces were eliminated from the equations of motion. In the complete Lagrangian approach, relative motions at the hinges were used as generalized coordinates. Accordingly, forces due to holonomic constraints did not appear in Lagrange's equation as they do not contribute to any virtual work. Ho developed the kinetic energy expression for individual bodies, using the inertial coordinates, with the body center of mass as a reference point. This was followed by a transformation to new coordinates which described the relative motions at the hinges. The system kinetic energy is the sum of the kinetic energies associated with the individual bodies. Ho indicated that with the Lagrangian procedure the resulting equations are the same as the ones obtained using the Newtonian-Lagrangian formulation. Hughes [35] used the direct path approach to formulate the dynamical equations for a chain of bodies connected by joints possessing one, two, or three degrees of freedom. 12 In contrast to the previous contributions, here the two bodies at the extremities are rigid while the interior bodies structurally flexible. This study was motivated by the Canadarm with the Space Shuttle and a payload at opposite ends. For the purpose of real time computation associated with the manipulator control, the equations were linearized in the angular rate and elastic deflections. However, formulation permits unlimited, slowly evolving angles at the joints. The Newton-Euler procedure was used in conjunction with the center of mass of each body. The sum of all the translational and rotational equations provided the system equations free from any constraint as explained earlier. The sum of translational equations and the sum of the rotational equations, for bodies on one side of a hinge in the chain, gave the constraint force and torque at that hinge, respectively. The constraint force thus obtained was substituted into the summed rotational equations. Similar relations were obtained at other hinges. Projection of these equations along the rotational directions resulted in a minimum set of the governing equations. An important configuration concept should be noted here : the path between any two bodies is unique. By analogy with a botanical tree, which has smaller branches proceeding from larger ones but no closed loops, such a system is said to have a tree topology. In a system with the configuration of a toplogical tree the constraints are holonomic. The degrees of freedom for a rigid system are equal to six plus the number of degrees of freedom at each hinge. In the literature discussed so far, the constraints are dealt with using the Lagrange multiplier method or the summation and projection process as explained. In the following, other approaches are reviewed where the system constraints are treated differently. Bodly, et al. [16] considered a general configuration of flexible bodies with closed loops. Six degrees of freedom were allowed between adjacent bodies. The dynamic equations were retained in primitive or free body form and supplemented by the constraint equations. The Lagrange multiplier method was used for both the holonomic and nonholonomic constraints. Instead of simultaneously solving the motion and constraint equa- 13 tions, a minimum set of equations was obtained as follows. The constraint equations are differentiated and together with the motion equations solved for the Lagrange multipliers. Expressions for the latter are, in turn, substituted back into the equations of motion. The resulting set of equations completely define the time derivative of a nonholonomic velocity vector. The latter is numerically integrated. Jerkovsky [36] considered a system of rigid bodies in the tree configuration. He started with the free body equations which were linearly transformed by an operator. In the original equations the velocities were inertial, whereas in the new equations the velocities are relative. In the case of closed loops, the Lagrange multiplier method has to be used. Jerkovsky describes his approach as a mixed formulation procedure since the final equations were expressed in terms of either intermediate variables, which were algebraic functions of generalized coordinates and momenta, or the time derivatives of the generalized coordinates. Such intermediate variables are kept for computational simplicity of the resulting equations. Keat [37], in his doctoral dissertation on multibody dynamics, developed a model accounting forflexibilityand closed loops. Following Boland [38,39] and Wittenburg [40,41], he used d'Alembert's principle of virtual work to obtain primitive equations, i.e., in terms of absolute velocities. Jerkovsky's transformation operator was generalized to account for flexibility. This transformation was used to introduce relative velocities at the joints, hence the holonomic constraints were satisfied. Yet another transformation was used to deal with the nonholonomic constraints. It was obtained by differentiating the constraint equations with the second time derivatives of the generalized coordinates expressed in terms of a reduced set. The transformation is essentially a form of the orthogonal complement of the matrix of nonholonomic constraints and, when substituted in the equations of motion, eliminates the constraint forces. Ohkami, et al. [42] are in a process of developing a general multibody program for the National Aerospace Laboratory, in Japan. The model, at present, accounts only for 14 rigid bodies in a general closed loop configuration. Each hinge is modeled to allow up to six degrees of freedom. An extension is underway [43] to include nonlinear hinge forces, system flexibility, etc. The Newton-Euler equations are used for the translational and rotational motions for individual bodies. The equations are supplemented by constraint conditions which express the relation between bodies' absolute velocities and hinge relative velocities. Additional kinetic constraint equations are used to model spacecraft docking. The equation of motion as well as the constraint relations are expressed in terms of time derivatives of nonholonomic angular and linear momenta with the constraint forces and torques expressed as vector x. The system of equations are arranged in an algebraic form Ax = B and solved for the time derivative of the momenta vector, constraint forces and torques. This is followed by a numerical integration routine for the time derivative of the momenta vector. The authors have discussed the inversion of matrix A. They referred to their method as the unified matrix approach since it accounts for both the kinetic and kinematic constraints. In connection with the human-body models, Huston and Passerello [44] developed dynamical equations of motion which were applied to simulate a crash-victim. The model consists of rigid bodies in a tree configuration, where each body is allowed up to six degrees of freedom. Interestingly, the approach proposed by Kane [23], which is essentially the Lagrange form of d'Alembert's principle, is used. This procedure automatically eliminates nonworking internal constraint forces. The overall system coordinate is fixed to a reference point in the main body, and local coordinates are fixed to the center of mass of each component body. The inertial linear and angular velocities and accelerations for each body are expressed in terms of its velocities and accelerations relative to an adjacent body, thus associating them with the generalized coordinates. Partial angular velocities Uki and partial velocities Vki are obtained for each body A; by partial differentiation of the corresponding inertial velocities with respect to the generalized relative velocity qi. All external forces and torques on body k are replaced by a force Rk through its center of 15 mass and a torque J * . Similarly, the body inertial forces are replaced by a force Rk and a torque T£. The generalized active forces associated with the coordinate / are calculated as Fi = S£ (vjfc/.i?fc + Uki-Tk), where n is the number of bodies in the sytem. Similarly, =1 generalized inertia forces are obtained as Ff = Ylk-i(^kl-Rk + ^ki'Tk ). Lagrange's form of d'Alembert's principle states that sum of the generalized active force and inertia force for each generalized coordinate of the system must be zero. Accordingly, the equations of motion are given by Fi + Ff = 0, where / ranges over the generalized coordinates. To summarize, a review of the literaure suggests that the problem of multibody dynamics has been approached using a variety of methodologies : (i) Primitive equations for individual bodies are obtained using Newton-Euler's, Lagrange's, or d'Alembert's method. The system coordinates are dependent, hence the motion equations are supplemented by the holonomic and nonholonomic constraint relations. An augmented body or a direct path approach may be pursued depending on the choice of the systems' coordinates and the management of the equations. (ii) Using Lagrange's approach, d'Alembert's principle of virtual work, or their manipulation as suggested by Kane, the equations of motion for the systems' coordinates are obtained free from any holonomic constraints; however they must be supplemented by appropriate equations describing nonholonomic constraints. (iii) The holonomic constraints are dealt with as follows : a) Use of the Lagrange multiplier method to solve the motion and constraint equations simultaneously. b) Based on physical understanding of the constraint forces, the equations of motion for certain individual bodies are summed and projected along the motion direction. c) Differentiating the constraints equations and solving them, together with the 16 motion equations, for Lagrange multipliers, which are substituted back into the equations of motion. d) Premultiplying the equations of motion by a matrix which is the orthogonal complement to the constraint matrix. The new reduced set of equations together with the constraint relations are then solved simultaneously. e) The equations of motion are linearly transformed using an operator. In the original differential equations the velocities are inertial, while in the transformed equations they are relative. The holonomic constraints are trivially satisfied. f) The equations of motion and constraint relations are expressed in terms of the time derivative of nonholonomic momenta. These expressions are arranged in algebraic matrix form, the solution of which gives the constraint forces and the time derivative of the momenta. The latter is integrated to obtain the required solution. (iv) The equations of motion for a multibody system have one of the following forms : p=F(q,u) Jii + Q, = G{q,u)+Q, q= WJ~ p; 1 q= W u • where : q generalized coordinates ; u time derivative of quasi-coordinate ; p generalized momenta ; Q generalized forces ; J generalized mass matrix ; W linear map of u to q ; F, G nonlinear inertia forces . 17 1.5 Outline of the Present Investigation With the vast body of literature as background, the present thesis attempts to take an important step forward. It develops a rather general formulation for studying dynamics of a large class of systems, consisting of arbitrary flexible deployable structural members, forming topological tree-type configurations, open or closed. Such a broadbased approach accounting for rotations, translations, deformations, deployment and their interactions has received virtually no attention, not to mention the development of associated algorithms for their evaluation. Furthermore, it achieves this entirely through the Lagrangian approach which has never been attempted before. Lagrange's approach is not popular in multibody dynamics because the kinetic energy expression can become extremely large and perhaps unmanageable, as indicated by Hooker [27] and others. On the other hand, its effectiveness has been attested by a variety of problems in analytical dynamics for more than 200 years. More specifically, the approach automatically satisfies holonomic constraints without any additional relations. It provides expressions for useful functions such as Lagrangian, Hamiltonian, conjugate momenta, etc., and the form of the governing equations displays clear physical meaning in terms of contributing forces. Equally important is the fact that the equations are readily amenable to the stability study and well suited for control design. A key to the use of Lagrange's approach in multibody dynamics is the development of the kinetic energy expression in a concise form, which can be differentiated as required. Obviously, the favoured form for the kinetic energy is | y M y , where y and M are the system velocity vector and mass matrix, respectively. Also the mass matrix should clearly display the system's dynamic character in a simple and meaningful form. Such a form for the kinetic energy is known for configurations such as a system of point masses, discretized vibrating structures, rigid bodies connected in a chain form [31], etc. To arrive at the form for a complex system under study would be the first challenge. 18 An emphasis on generality is one of the distinctive features of this thesis. It is achieved by considering an arbitrary number and type of mass elements connected through general joints to form any desired configuration. To be more specific, the model considers a number of arbitrarily connected members such as membranes, beams, plates, shells, etc., which form a configuration with branches and closed loops. Each member is modelled as an axially deployable flexible body. Joints between the bodies permit up to six degrees of freedom in translation and rotation. The system is divided into domains, each containing one or more members whose motion may be holonomically constrained. Necessary cuts are introduced such that the system has a branched configuration with constrained domains. A direct path is defined which uniquely describes the position of the origin of the inertial frame with respect to the system domains and subdomains. The configuration is described using the matrix of geometry thus extending the application of kinematics to such an involved multibody system. The distinction between local coordinates associated with a constituent body on the one hand and the multicoordinate overall system on the other hand allowed the progress of matrix formulation in a logical sequence. Furthermore, it helped identify, with clarity, the characteristic dynamical parameters associated with the system resulting in further appreciation of the contributing terms. Two sources of constraints are distinguished in the model. They are : (i) constraints among a set of subdomains forming a given domain; (ii) constraints among the system domains due to cuts introduced in the closed loops. In general, these constraints are holonomic, hence they are incorporated in the kinetic energy expression by using a transformation matrix, which maps the dependent coordinates into the independent ones. Any nonholonomic constraint in the system must be considered together with the equations of motion. The well-known Lagrange multipliers can be used, however, it is preferred to obtain a reduced set of differential equations for a numerical solution. This can be achieved through premultiplying the equations of motion by a matrix which is the orthogonal complement to the matrix of constraints. Note that the motion of domains forming 19 a branched configuration are free of any constraint relations, since the generalized coordinates are used which represent orientation of the individual body relative to its neighbour on the direct path. The comprehensive character of the formulation can be appreciated by the following observations : (i) It is believed that a formulation of comparable generality has not been recorded in the open literature. For example, it permits modelling of deploying appendages using quasi-modal as well as discrete-coordinate representations. The system reference can be its center of mass or an arbitrary point. The matrix of geometry provides considerable flexibility in the choice of the system coordinates thus synthesizing the Newton-Euler and Lagrange procedures to any desired degree. (ii) The formulation allows a varying degree of sophistication in modelling. For example, an element in the system can be modeled as a point mass, a rigid body or a flexible member. Flexibility can be discretized using lumped parameters, admissible functions or the finite element method. The system nonlinearities are preserved in the formulation, thus permitting a degree of simplification as warranted. (iii) The formulation accounts for ignorable coordinates and specified motions. (iv) An expression for the gravitational potential energy is derived, which is new to the multibody system. It accounts for the nonlinear contributions to the order 1/TQ, due to the inverse square gravitationalfield,thus making the formulation applicable to very large spacecraft, such as the proposed Space Station. Here ro is the distance of the system from the inertial frame. (v) Structural flexibility may require special treatment for adequate modelling. This is achieved by permitting a second-degree nonlinear displacement field. Alternatively, a linear displacement field can be used together with a fourth degree nonlinear strain energy. Thus the stiffness matrix contains terms up to the third-degree of 20 nonlinearity. (vi) It is now possible to give the governing equations of motion a variety of different desirable forms, a concept quite new to the multibody dynamics. Also, characteristic dynamical functions are obtained quite readily. As pointed out before, the holonomic constraints are satisfied during the process of formulation. (vii) Management of contributions from various sources, during development of the formulation, using concise matrix format provides better physical appreciation of the constituent terms. Equally important would be to demonstrate application of this versatile formulation to the solution of problems of practical importance. To that end four configurations of contemporary interest are considered for dynamical and stability studies : (i) The rigid Space Shuttle, in an arbitrary orientation, negotiating a circular trajectory. (ii) The orbiter based deployment of a beam-type appendage. This has some similarity to the experiment carried out by astronauts Jerry Ross and Sherwood Spring in November 1985. Obviously this is one of the fundamental steps in construction of the proposed Space Station. (iii) The Space Shuttle based deployment of a plate-type appendage similar to the Solar Array Flight Experiment (SAFE) conducted in September 1984. (iv) The Space Shuttle based deployment of twin beam-type appendages (each 150 m in length) similar to the Waves In Space Plasma (WISP) experiment jointly under investigation by Canada and the USA. The configurations suggest progressive introduction of complexity to help appreciate parametric interactions involved. The study establishes a formulation with a potential in tackling a large class of dynamical problems. 2. K I N E M A T I C S 2.1 OFA MULTIBODY SYSTEM Modelling of a Multibody System Modelling of a dynamical system is, of course, not unique. A mathematical model for a system may depend on the physical problem, purpose of the analysis, available computational facilities, and similar other factors. In the present development, several alternatives are allowed in the model construction which are subject to engineering judgement. Consider a system consisting of a number of arbitrarily connected members, each of which modelled either as a point mass, a rigid body, a flexible body, or a flexible deployable member (Figure 2-1). Let the system be divided into domains, each containing one or more members, whose motion may be holonomically constrained. In general, the system domains form a topology of branched and closed loop configuration. Through appropriate introduction of cuts at joints (Figure 2-2) it is possible to arrive at a branched system with a unique path, referred to as the direct path, to any point in the system starting from the origin of the inertial frame. To account for the system continuity at the cuts, a set of constraint relations must be satisfied. Some domains in the system may be further divided into holonomically constrained subdomains. In that case, motion of the latter is described relative to coordinates which are denoted as global (Figure 2-3). In particular, a flexible structure may be divided into members or finite elements where continuity should be satisfied at their boundaries. Let one of the domains be taken as the main or home domain (Figure 2-1) with an arbitrary point on it denoted as the main point. In general, the main point can be approached from the origin of the inertial reference frame through an intermediate reference point. Thus, the main point may be looked upon as a connection between the main domain and the reference system. Domains at the far end of each branch in the system are referred 21 22 sub domains terminal domains terminal joints inertial coordinates Figure 2-1 Multibody domains and subdomains. Note the typical direct path to a point mass in the system depending on its location. 23 Figure 2-2 A cut is introduced to the closed loop to define unique paths to points in a system. 24 s u b d o m a i n coordinate coordinate - - ^ ( x 4 , V4, Z4) ( z s , VsiZs) inertial t e r m i n a l c o o r d i n a t e Figure 2-3 Multibody coordinates. 25 to as terminal domains with the corresponding joints denoted as terminal joints. Starting from the main point a path can be defined to any of the terminal joints. One can achieve this through a minimum number of intermediate connections. On each of the subdomains a reference point with a corresponding coordinate system is defined which is connected to the origin of its global coordinate system (Figure 2-3). The totality of the connecting lines, initiated from the main point, form a topological tree configuration. Each line representing a unique path from the origin of the inertial frame to any of the system domains and subdomains is called the direct path. Let X, Y, Z (LN) and xo,yo,^o (Lo) he the inertial and reference coordinate systems, respectively. Fixed to the main domain D\ (body B\), at the main point, is the coordinate system xi,y\,z\ (Li). Similarly,fixedto each domain and subdomain D{ (body Bi) at joints on the direct path are local coordinate systems The rigid body motion of a domain Z),- is defined by the degrees of freedom of its joint. This is represented by the motion of the coordinate system L,- relative to the adjacent domain on its direct path. An exception is the motion of the main domain D\ which is defined by the motion of the reference coordinate Lo relative to the inertial frame LN, hence the reference point works as the joint for the main domain. The motion of a point in the domain, due to deployment and flexibility, is described relative to the domain coordinate Li. The motion of a group of subdomains forming a domain is described relative to its global frame. Note that the motion of the system domains in a branched configuration is free of any constraint relations, while the motion of a group of subdomains is not. 2.2 Kinematics : Single Body with Local Coordinates Here kinematical relations are developed which describe, relative to local coordi- nates, the position and velocity vectors for a generic point on aflexible,deployable member. As an illustration, consider a single flexible body Bi, deploying as shown in Figure 2-4. Fixed to it, at an arbitrary refrence point Oi, is the local coordinate system X{,yi,Zi. Let 26 Figure 2-4 Position vector, in local coordinates, to a point mass on a flexible body with deployment. 27 the position of a generic point on the body, in its initial undefiected state, be described by vector p with projections p\,p ,p % l 2 3 along the body coordinates x -,yi,2,-, respectively. Now t let the generic point on the body undergo an elastic displacement w*, which is a function of time varying parameters S{ and the position vector of the point p*. The spatial functions, regarded as preassigned, have continuous derivatives and satisfy at least geometric boundary conditions. The time varying parameters Si are unknown and represent the generalized coordinates. This representation of the deflection is a special case of a more general form where the deflection vector w may be expressed as a power series in the time 1 varying parameters Si, i.e., the first term is linear in Si, the second is quadratic in Si, and so on. This assumed form of the elastic deflection w is in agreement with the engineering 1 practice and has its roots in the theory of elasticity. As discussed later, under certain restrictions, quadratic and higher degree contributions can be expressed as functions of the leading term. The functional dependence of the deflection vector may thus be written as w = 1 tZJ*(/3\ Si), where /3 l is a function of the spatial vector p and the body dimensions. Note, 1 with deployment the body dimensions are changing with time. Let the body coordinates (i.e., Xi, yi, Zi) be so selected that one of the axes is aligned with the direction of deployment. Let the component of p along the direction of deployment be p\. The deployment velocity 1 is represented as ti. Denoting the r-th element in the array of generalized coordinates <5t as 81, the position and velocity vectors of the generic point , in the local frame, can be written as : £ =p l + w ; (2.1) l dp\ dt ' 28 = li rf + F iT (2.2) Si ; where t _ dp*' t duT ~ ~dp~\ + dp* / dw 1 ~dp\ + dp W ~dU ' A,- x duT 1 = di + ( 1 + j? ) UT"' ; (2.4) P= (2.5) Here i i / = dw /dp^ and d,- denotes the vector of direction cosines of the deployment direc1 1 tion in the body coordinates. Note that, because of deployment, the total time derivative —t of the position vector £ contains a convective term. The special cases follow quite readily with Bi modelled as a : • point mass: t = 0; £ == •t 0; • rigid body: f = t i i -- • rigid deploying body: f = t e == U di ; • flexible body: f = t+ w ; • flexible deploying body: f = f + v? ; 1 = 0; £ == F •i iT ti ; = Urf + F iT ti In absence of rotational and translational motions, a linear displacement field is usually used in conjunction with admissible functions or afiniteelement analysis. However, an adequate modelling of the problem may require the quadratic term in the displacement field. For example retaining quadratic terms in the transverse displacement field accounts for the axial foreshortening due to the transverse vibration as discussed by Vigneron [45,46]. 29 Concerning the problem of deployment analysis, there are two distinct approaches suggested in the literature : (a) discrete-coordinate formulation; (b)quasi-modal formulation. Note, the model presented here accounts for both the procedures. The first approach has been used in the early rigid-multibody-models for applications in space dynamics. Here an appendage is modelled as a collection of point masses or rigid bodies. The approach can be useful when large appendage deformations are expected. In the quasi-modal approach the appendage flexibility is discretized using vibration modes with time varying length, often referred to as quasi-modes (with associated quasi-frequencies) to distinguish them from the natural ones. Obviously, this is an approximation as for a system with time dependent configuration, as in the case of deployment, natural frequencies and modes are not defined. Instead one visualizes a succession of natural modes and frequencies for instantaneous lengths during deployment. The resulting equations of motion contain time-dependent coefficients representing instantaneous member length, mass, extension rate, and stiffness. A t each instant, the equations give correct modal representation, thus promising the solution to converge to the true value when the step size of integration is sufficiently small compared to the extension rate. A review of the rather limited literature on the deployment dynamics suggests that the approach has been used by Cherchas [47], Lips and Modi [48], Hughes [49], Ibrahim and Misra [50], and others. 2.3 Matrix of Geometry Equations of motion are usually derived for a specific system with known configu- ration. In the present task of a general formulation applicable to a large class of multibody systems the configuration is unknown. This added difficulty dictates development of equations using parameters which can describe the system geometry. In the present development, this is achieved by introducing a matrix of geometry, S. A n element of the matrix is denoted as Sj , where the superscript refers to domain or coordinate number i while the subscript denotes coordinate j. Recognizing that attached to each domain is a 30 coordinate system with its origin uniquely defined on the direct path, it is apparent that for each pair of domains there is a uniquely defined element in the matrix S . Also note that the geometry of a branched system is completely defined once the presence or a lack of direct path between domains i and j is established. Thus the geometry matrix element Sj responds as true or false depending on whether j is on the direct path to % or not. Here, number 1 is asigned for 'true' and 0 for false : 3 _ f 1, \ 0, if j is on the direct path to i ; otherwise . , » \ •) 5 is a square matrix of length equal to the number of coordinates in the system. Each of the diagonal elements has value 1, since each domain lies on its own direct path. Since all direct paths originate from point 0\ on the main body B\, when j equals 0 or 1, 5j takes the value 1 for all *. The general rule to remember when using the matrix S is that quantities multiplied by Sj are to be regarded or disregarded depending on whether Sj has the value 1 or 0, respectively. It may be pointed out that when Sj appears inside a summation operation over i or j , again the algebraic value defined by Eq.(2.6) is used. On the other hand, when Sj is not operated on by a summation symbol, it sets a condition on the validity of the expression. The expression is valid only if Sj equals 1. An illustration of the geometric matrix describing a multibody configuration is presented in Figure 2-5. 2.4 Kinematics of a Multibody System From the mapping point of view, which is adopted here, physical quantities of interest, such as position and velocity vectors, are expressed in one frame or the other through a transformation matrix. These transformations for a multibody system can become rather involved. Here the process is systematized and operations carried out step by step to avoid confusion and possible errors. Bi 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 1 1 0 0 0 0 0 1 0 1 1 0 0 0 1 0 0 1 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 1 0 1 0 0 0 0 1 0 0 0 1 1 0 0 1 Figure 2-5 The matrix of system geometry. 32 2.4.1 Position vectors As before, ^ is a position vector for a generic point on B{ in the local coordinate system Xi,yi,Zi, which is denoted as frame L,- in Figure 2-6. With respect to the body frame Lj, with coordinates Xj,yj,Zj, the position of a point on B{ is denoted as J a\ Consider a transformation matrix T- (Li (->• Lj), which maps the position vector the local coordinate system Li into vector a 1 in in the local system Lj. The mapping is an 1 orientation and metric (magnitude) preserving transformation, therefore it can be uniquely decomposed into a translation and a rotation, ya"i : Lj »-» Lj and C\ : L ; i-> Lj, respectively, hence S) >tf = S} T> f = Sj {jcii + c{ f} (2.7a) . The transformation matrix T / can be decomposed as : T( = Tj Tj^....T^i +l The TJ' . 1 plus or minus sign corresponds to changes along the direct path and not to an algebraic operation. Operating on £* using the transformation matrix in the decomposed form leads to : S}ia i = S) Ti_ x = S}T>_ T\~ f x (ai + C r 1 = S}Ti_ Till(ai 2 {ai_ 1 f) , + & r \ ai + C\Z\ C ^ x 1 f) (o-_ 2 + C\zl o,--! + C\Z\ ai + = S}Ti_ 3 = S} ( a f), 1 + C t 2 = S]Ti , j + 1 + C j j + l a j + 2 + + Ci_ , C i - z f ) , o-_ 2 + C | L 3 + 2 a ) + S} C{ f { • Here, a; is the position vector to origin Oi of the body frame Li in frame L»_i. Note, the above decomposition is performed along the direct path. The expression can be rewritten in the form of a summation over all the system coordinates and, using the geometry matrix S, over the appropriate coordinates on the direct path as follows: n Sj ' V = S} ] T St S j - C\_ 1=0 1 x ai + S*j C\ f , (2.7b) 33 Frame L i - i Figure 2-6 Position vectors in a multibody system. 34 = SJ £ t S +1 5 J _ 1 L «i; C =.-+1 = t • (2.7c) Comparing Eq.(2.7a) and Eq.(2.76), n s i- &= 5 J E s t s i~ l c i - i • (- ) 2 7d 1=0 The first term on the r.h.s. of Eqs.(2.76) and (2.7d), contains geometry matrix elements Sf and 5 J . The former sets the condition to include only those coordinates which are on - 1 the direct path to i. The latter, S J , requires exclusion of the coordinates on the direct - 1 path to j. By projecting of on Lk and using Eq.(2.7) we introduce the vectors ^a and ya,3 1 which are defined by : SJ{a = SJC a i 3 (2.8a) = SJ a +C?f, i k j i n - Sj Y, ! S 5 J _ 1 1 c C i - i * i + <tf r , (2.86) 1=0 n = t s E s 1=0 i + = ' - i « ' » 5 ( - ) 2 8C n z=o where 5 j Ja,- = 5 ; C) -aZi . (2.8c) As a special case of Eq.(2.8), the projections of a? and ya- on the inertial frame L J V are 3 t denoted by r* and yr,, respectively, and can be written as follows : J 5; r = Sj Cj a = SJ jfi + d f , 3 i 3 (2.9a) 1 n = Sj Y, Sj Sj' C / _ a , + C - I*" , (2.96) 1 x t 1=0 n = S J X) 1=0 5 / + 1 c ' - i ' ^+1 = ^ : ( 2 9 c ) 35 n where yr- = 5} Cy ya- . t (2.9e) t Note that : J V= jo ' ; 1 ya- = Jo,- ; r — a , t N - _ TV- . j't — y <*i , - _ t-l- . o.i = i-i^t ; - _ ~i ~ NNi -t' _ r = ) a N-i N a The convention used here is explained below : • The general form for expressing a position vector is j.0 or *Oi, where the former 1 defines the position of a generic point on S - relative to Lj and as projected on Lk . t The latter represents the position of the joint at C\, the origin of frame Li, relative to Lj and as projected on Lk• When k = j the general forms described above are denoted by •'a* and ya"i, respectively, for simplicity. • When k = N, where N corresponds to the inertial frame, the general forms are denoted by f* and yr,-, respectively. J • When k = j = i; the vector ^a is denoted as £*, and when k = j = N it is 1 designated as f*. • When k = j = i — 1 the position vector ya,- is simply written as a -, and when t k = j = N it is denoted as r,-. Furthermore, it should be noted that Eq.(2.9d) takes special forms in the following 36 four cases : (i) 3 > 1 (frame Ly is a local body frame), ij S S 3*i = Y, l Sj ' S • 1 1 (2-10a) l=i (ii) j = 1 (frame Ly is the main body frame), n iri = E S} C , _ i ai . (2.106) 1=2 (iii) j = 0 (frame Ly is the system reference frame, Figure 2-6), n or; = Y SJ C , _ i o, . (2.10c) l=i (iv) j = — 1 (frame Ly is the inertial frame) : r ' = T.-f = r + Cif , (2.11a) 1 { n = X)s/C,_ia, + C -?, z=o (2-116) t n = £ C , _ ! a,, o, +1 =? ; (2.11c) n rt = J ^ 5 / C j _ i oi . (2.11d) From Figure 2-6, it is apparent that Vi = Sf{ri + ,r<) ; l i fi ^ S f{¥i+ r ) . (2.12a) (2.126) For frame L i , the only coordinates which affects its motion are those on its direct path. This condition is enforced through the use of Sj on the r.h.s. of the above two equations. 37 2.4.2 Velocity vectors The time rate of change of the rotation matrix Cy is obtained as ^ = X:s/c,_ c;-'cj. 1 In general, if matrix C represents rotation of frame L A relative to frame LB, C = CU, 6J = C C, T or where u is a skew symmetric marix formed from the angular velocity vector of L A relative to LB as projected on the former frame. Introducing this relation into the above equation dt 1=0 n Citif = Y,S{ z=o C) . 1 (2.13) The inertial velocity of a generic point on Bi is obtained by differentiating Eq.(2.11c), re ? =E n s i+1 3 c ;-i +E j 5 y=o +1 ^'-i > = T• (- ) 2 14 y=o The second term on the r.h.s. of the above equation can be developed further by using Eqs.(2.13) and (2.9c) as follows : £ y=o C y - iffy= £ S} j=o £ +1 « r 1 [ C, wj" C J ^ ] ay , 1 +1 = f , i=o = j r c ul- Cf [ J2 1 1 t n = £ ay 5/ Cy-l ay ] , ay = f , +1 c> u j - c r V , 1 n = £ 5/ ( Ci ujj- ) xr . 1 (2.15) l 1=0 Substituting from above into Eq.(2.14) results in re F* = F ' + £ 1 5/ ( Ci oJJ- ) x ¥ , 1 (2.16a) 38 where = £ j S +1 i-i C *i ' = ? > ( 2- 166) y=o re Sj"' Cy_i Oy + C - £ , t (2.16c) y=o •t r,- + CiC ; (2.16d) and rt- = J2 Sj c 3-i *i • ( 2- 16e) 3=0 In Eq.(2.16a), r* represents the inertial velocity of a generic point on Thefirstterm on the r.h.s. in Eq.(2.16a) represents the contribution of the translational velocity while the second term is due to coordinate rotations on the direct path to body i. In Eq.(2.16d) the translational velocity, r*, has two components: the first, f,-, is the result of the translational velocities of body frames on the direct path to i, while the second term is due to the motion of the point on B{ relative to the body frame L,-. 2.5 Kinematics of Flexible Multibody System Two methods are considered here for discretization of structures. In thefirst,the flexible body is discretized as a number of interconnected elements (subdomains) and the deformation within each element is described in terms of generalized coordinates with respect to a common global coordinate system. The continuity of the displacement field between the elements can be assured by imposing constraint relations. Here, this approach may be used when the flexible elements are not connected through hinges and are assumed to have small elastic rotations. Actually, with motions described relative to a global coordinate, large elastic rotation may be accounted for only on the basis of the full nonlinear strain-displacement relation, which is discussed in a later section. If the elastic rotations are large or the bodies are connected through hinges with 39 rigid degrees of freedom, another approach is used. Here each body is represented as a separate domain, where its rigid motion is defined relative to the adjacent body along the direct path. The rigid mode is described by the motion of the body coordinates, which are located at a hinge point. This motion will have contributions from the hinge degree of freedom as well as the elastic motion of the adjacent body, at the point of connection. Also, there can be relative translational and rotational motions of the body at other joints, which are caused by the interacting forces and torques from other adjacent bodies. Here the motions are relative and, hence, free of constraint relations. However, the displacement field of a typical flexible interconnected body, with at least two or more hinge points, should satisfy time-varying boundary conditions at the joints. For example, a special case of this procedure is applied in the finite element method to account for large elastic rotations of a beam (in the absence of rigid body motion). Here the beam is subdivided into a sufficient number of elements, hence the elastic rotation of the element can be restricted within each element relative to its own reference frame to any desired extent. In this case, member elastic properties may be obtained on the basis of small rotations and strains. Consider body Bk of the multibody system shown in Figure 2-7. Fixed to it at point Ok is the local coordinate Lk. The body adjacent to Bk on its direct path is denoted as Bj with local frame Lj fixed at Oj. Body Bk is connected to body Bj at point Ok through a joint hki which may have rigid degrees of freedom in translation and rotation. The terms defining location and orientation of Lk relative to Lj will have contributions from the deformation and/or deployment of Bj, as well as from the joint degrees of freedom. Thus, the translational motion of Lk realtive to Lj will have contributions from hinge h-k's translational degree of freedom and from By's deformation and/or deployment. Similarly, the orientation of Lk relative to Lj will have contributions from hinge h^s rotational degree of freedom as well as the slope (elastic rotation) at hk due to the deformation of Bj. Assuming that a general expression for the displacement field of a flexible body Bj with time-varying boundary conditions at the joints Ok is available, expressions for 40 Figure 2-7 Intermediate coordinates for hinge translational motion. 41 the translational and rotational motions of a frame Lk fixed to the body are developed first. This is followed by a discussion on simplification in the displacement field leading to particular cases. The position of point Ok on body Bj in its initial undeflected state is described by vector p = {p }ok J k a n o ^ its elastic displacement is defined by vector Wk = {w }o - For 3 k each frame Lj, we define a conjugate frame L*j (Figure 2-7), fixed to body Bj with its origin denoted as Oj. Its matrix of direction cosines with refrence to frame Lj is denoted as Rj. The translational motion of frame Lk due to hinge hk's degrees of freedom, is measured relative to the conjugate frame Lj, and is denoted as Uk- Note that frame L*- will translate relative to frame Lj according to the deployment of Bj. In light of Eqs. 2.1 to 2.5, the position and velocity vectors of point Ok, in the local frame Lj, can be written as : ak = {I }o 3 = {p + w }o » } k 3 k - {p }o* + Rj ^k + {w }o 3 3 k , = Pj*+ Rj Uk + Wk , (2.17) = ij W + k (2.18) F k T Jy + Rj *k ; The first and second terms on the r.h.s. represent translation of body B due to deployment k and vibration of body Bj, respectively. The last term describes the translational motion of body Bk relative to body Bj, where : (2.19) = dj + {i + $ ) k K ; (2.20) 42 Here, wJ^. = •! dw /dp , \ 3 , and dy denotes the vector of direction cosines of the deploy- 3 } °k ment, in the body coordinates. Thus, for clarity, relevant parameters are : Py = position of O* in frame Ly with body Bj undeformed ; pj = Pj + Rj Uk, location of Ok in frame Ly with body Bj undeformed ; Rj — matrix defining rotation from L*- to Ly ; Uk = location of Ok with respect to frame Ly with body Bj undeformed ; u>k = {w }p , value of w evaluated at p ; 3 3 k k Fk = [ F ]p , value of F evaluated at p ', 3 3 k k rf — {rf }p J h i e of rj evaluated at p • 3 3 va k k k Now the rotational motion of frame Ljt relative to frame Ly is to be developed. With point Ok as origin, consider another coordinate system x i, y i, z i (Figure 2-7) fixed k k k to Bj. Now C^, is due to the elastic rotation of Bj at point Ok, and C k is because of the rotational degrees of freedom at the hinge. Hence the matrix of direction cosines between the coordinates Xk,Vk,Zk and be written as C{ = C{, C f . (2.22) C{ = C , C f + C{, C f Also 3 k C{ "I = C> , SJ, C f + C{, S# ; k u{ = C , cb{, C with k nl = c£,u , k + u + k> / k k , i.e. . (2.23) Here U is the angular velocity of Lk relative to Ly and as projected on the former. k Note that U k is the angular velocity of Lk relative to Ly due to hinge degrees of freedom, and Vpy is the angular velocity of Ly relative to Ly because of the body By's elasticity. Recognizing that k> = k> "i> ' C C 43 Thus it follows that i-th. element of u , can be written as 3 k £ = ligh C ,(m,h) C ,(m,g) , 3 3 k k m=l where rji h is equal to unity when the subscripts appear in a cyclic order; otherwise it is g zero. Now consider the expression for C , for a particular case of relevance. When elastic 3 k elongation and shear contributions are small compared to the elastic rotation in Bj, C , 3 k can be written as [51] (2.24) where d/dp 3 3 is defined as d/dp n , n = 1,2,3 . Also 3 d d w '3 d p 3 Note, d w jd t — F 3T 3 3 8 . Here F , y represents the partial differentiation of F with Ph 3 3 respect to p n (therc-thelement of p ). denoted as f ] Since F 3T 3 is a 3 x n s matrix, its ra-th row is . Hence, T rn C ,(m,g) = 6 P a ' k> 3 k J p Thus the i-th component of u> , is, k ^i'W = { £ ^igh i>( > ) C m h {/my h, k } • m=l The expression for the angular velocity vector due to elastic rotation can be put in the form &[, = G tj , (2.25) k where G is a 3 x n £ matrix. The i-th row of G is equal to {52m=i Vigh C f(m, h){f , y }- } 3 k k J m t Pg k > f From Eqs.(2.23) and (2.25) U = W 3 k k k +C k y Gk 8j , j = k-l (2.26) 44 The two terms on the r.h.s. of the above equation represent the angular velocity vector of frame Lk relative to frame Lk-i due to hinge (h-k) degrees of freedom and body's (Bk-i) elasticity, respectively. The problem now is to find a proper expression for the displacement field w which 3 satisfies the boundary conditions, including those of the moving boundaries, consistent with the simplification of small shear and elongations, and large elastic rotations. Instead of looking for an expression for such a complex three dimensional displacement field, the elastic deflection w is expanded in a power series with the leading term describing, for 3 example, the deformation of the elastic axis of a beam or the middle surface of a plate. With the shear and elongation assumed negligible, the higher terms in the series may be expressed as functions of the leading term, which is assumed known. Such expansions are known for plates and shells and, though more complicated, are also reported with additional assumptions for beams [51]. Note, these reported displacement fields, expressed as series, do not account for the moving joints. Hurty [52] developed a method for describing the kinematics of aflexiblebody, where the displacement field is composed of the rigid body's constrained and normal modes. It is clear that, in general, six independent rigid body displacements exist (motion of Lj). Hurty refers to a constraint as 'redundant' when it limits the body motion beyond the rigid body displacements (motion of Lk). Added to the rigid body contributions is a set of displacements produced by giving each redundant constraint, in turn, an arbitrary displacement keeping all other constraints fixed. The resulting deformations are called the constrained modes, and there are as many of them as there are redundant constraints. The rigid and constrained modes permit description of any arbitrary displacement of all movable joints. However, it is still necessary to provide for the displacement of other points on the structure relative to the constraints (joints Ok). Such relative displacements are provided by the introduction of a set of independent modes in which all constraints are fixed [Oj and Ok). An infinity of them exist and are called natural or normal modes. 45 In connection with the multibody dynamics Ho and Herber [18] have further divided the constrained modes into, what they refers to as, quasi-static translations and quasistatic rotations. To obtain a quasi-static translation mode, each redundant constraint in turn is given an arbitrary displacement but no rotation. Similarly, a quasi-static mode in rotation is realized by introducing pure rotation at each redundant constraint in turn. With this type of discretization, the displacement field w , in the above development, is J now represented as a sum of displacements due to the natural and the redundant (the quasi-static translations and the quasi-static rotations) modes. Only the latter contributes to the translation and rotation of an attached body, as described by Eqs. (2.17), (2.18), (2.24) and (2.25). 2.6 Inertial Versus Relative Coordinates Consider a multibody system which, upon introducing some cuts, has a tree config- uration as shown in Figure (2-8a). Alternatively, one may introduce cuts such that each body in the system is completely separate as indicated in Figures 2-8b, c. The two alternatives correspond to two approaches discussed here. In a tree configuration, translational and rotational displacements are described in terms of the relative motion between two adjacent body coordinates. The motion of separated bodies, on the other hand, is with respect to an inertial reference, and hence constrained. If the Lagrange procedure is used to derive governing equations for the system in terms of the inertial coordinates, they will contain constraint forces as described by the Lagrange multipliers. They correspond to reaction forces in the Newtonian procedure due to construction of a free-body diagram for each member. Hence, one may refer to the two approaches (Figures 2-8a, b, c) as the Lagrangian and Newton-Lagrangian, respectively. One may approach a problem here through either of the procedures, or even their combination depending on the character of discretization (Figure 2-8d). Irrespective of the discretization procedure used, one can find direct paths to all 46 Figure 2-8 Illustrations of discretized system topology leading to : (a) Lagrangian formulation ; (b) Newtonian formulation (inertial reference point) ; (c) Newtonian formulation (reference point is system cm.) ; (d) Newton-Lagrangian formulation . 47 constituent bodies, and hence the matrix of geometry can be determined as discussed before. The multibody kinematics may describe inertial or relative motions depending on the method of discretization used, and hence on the matrix of geometry. For example, if each body in the system is considered separately, as in Figure 2-8b, elements of the matrix of geometry S*- are all zero except for j = i, when they are unity. Using these values, for example, in Eq.(2.13) or (2.16e) the system velocities are inertial. Actually, a matrix of geometry which is unity suggests motion with respect to inertial coordinates. Note that equations like (2.11), (2.13) and (2.16) represent transformations from inertial to relative coordinates and, for unit geometry matrix, lead to a trivial identity. The Lagrangian approach, in conjunction with a system of branched configurations, involve generalized coordinates free of any constraints, which is favoured here. So the approach to discretization of the sysem, as explained in Section 2.1, mostly aims at representation in the form of a branch configuration. The kinetic energy expression developed in the next chapter can be a mixture of the two approaches. The matrix of geometry together with the system constraints define the configuration of the system. The matrix of geometry, system constraints and kinetic energy thus uniquely define the dynamics of a given problem. 3. K I N E T I C 3.1 PARAMETERS Kinetics of a Single B o d y Previous researchers in the field of multibody dynamics did not find the Lagrangian formulation attractive. As pointed out before, this is because of the wildly unmanageable form of the kinetic energy expression in most problems of practical importance. Hence, the goal here is to express the kinetic energy of an arbitrary multibody system in a quadratic form of the system velocities. The associated mass matrix should have, as far as possible, a clear and concise form suitable for differentiation as required in the Lagrangian procedure. This section identifies basic kinetic quantities, associated with the mass matrix, which characterize the dynamics of a single deployableflexiblemember. These dynamic quantities later serve as building blocks during the development of kinetic energy expression for a multibody system. In the Lagrangian mechanics, all the dynamic characteristics of a body are contained in its kinetic energy expression. Consider a flexible deployable body i?,-. Let a coordinate system, fixed to it at a given point, undergo arbitrary translational and rotational motions relative to frame L A (Fig. 3-1). The velocity r*, in the coordinate system L A , of a differential mass element dm,- in body J5,- is given by 1 r = o\- + C, f + Ci 57,- x | * , where a,, C,- and w,- define the location, orientation, and angular velocity of frame L,—i relative to frame L A , respectively. £ is the position vector of dm,- in frame L,-. Note, the above expression can be obtained from Eq.(2.16) by taking / = i. The kinetic energy of 48 49 Figure 3-1 Motion of mass element in generic body B{ relative to frame L A 50 body Bi with, respect to frame L A can be written as Ti = \ f ?T 1 J ™i f drrn , = i | m - of + 2 aj Ci J t +Ef/ ? J mi £ drrii + x j f dm,- + s f £2 / J mi drrii + aj [r Ci aJ* / «7 " f ? T f* drrii ] 57.-) , J (3-1) where U is a 3 x 3 unit matrix. The following expressions can be readily identified to have designated meaning : • £* drrii '• first mass moment of Bi in Lt- ; fm fm- £ drrii '• momentum of Bi in the local frame L - ; t —t — 1 f'mm. . t £ • | fm £ drrii '• x 2 m t £ 2 drrii '• kinetic energy of Bi in Li ; im-[£* ^ ' moment of momentum of Bi in L - ; ^ »" — m i : moment of inertia of Bi with reference to L t By analogy between the potential energy and its complementary form on one hand, and the kinetic energy and its complementary form on the other hand, Tabarrok [53] showed that | m i ; is strictly a complementary kinetic energy. Accordingly, 2 as defined in Eq.(3.1) represents a complementary kinetic energy. However, the magnitude of the kinetic energy and its complementary form are equal in classical mechanics by virtue of the linear relation between momentum and velocity. Hence, as it is customary, Ti is referred to here as the kinetic energy. Using the definition of the position and velocity vectors of a generic point in body Bi with respect to L - as given by Eqs.(2.1) and (2.2), the above integrals can be rewritten t so as to isolate the velocity terms as follows : 51 / —1 / £ drrii = m,- b ; Jm{ —x I f / f x f drm = mi {lipi f f.f dm = rm { 1] < / [ ? ' drrn = mi{iiK + # P + Pi r ^ - f f ti } ; T + 2 U wi { r h }5 t + of W T s { ss t { }; (3-2) [r]; T l where 6* = — i m Ti ~ e H l s I Wi dmi ; ™t Jmi — I m,I — mi Pi = — < = m Jm. T Jm. x dmi i dmi • Vi ^F dmi; iT I i *fi Wi dm ; { Jmi — mi WU = — i m P= { Jmi — wis = F f m Ke f* dm, ; Jm{ — = Pe = / J J m; [ f Jm. Vi p l dmi ; T [ F F dm ; i iT i Jm{ T f U - f f T } d m i . (3.3) 52 It may be pointed out that, in general, f is a three dimensional position vector and S{ is an array of length n , which is equal to the number of generalized coordinates associated —t — i l s with the deflection of B{. Note also that iu* e is a scalar; b , h , and are three dimensional e vectors; and w ^ is an array of length n\. The two matrices H\ and Pg are of dimension 1 ra£ x 3, while P and W$s are square matrices of dimensions 3 and n , respectively. Further % s details and physical meaning associated with the terms are elaborated upon in Appendix A. Substituting from Eq.(3.2) into Eq.(3.1), the kinetic energy of f?,- can be written in the matrix format as / — ai Ti 1 2 T -\ U Ch { CiH e h <. Si W db T iT 6 ai T eS k > pi Sym. Si ™i Note, the conventional 'mass matrix' associated with the kinetic energy expression is now dependent on the system coordinates as well as the amount of matter contained in the body. The mass matrix is divided into submatrices, which are associated with system velocities and their interactions. Since any general motion can be described as a combination of translation and rotation, the mass matrix can accordingly be divided into mass, first mass moment and second mass moment. The translational motion may arise from a variety of sources including conventional rigid body motion, elastic deformation and deployment. Thus it is apparent that there are several forms for the mass and first mass moment; and the structure of the mass matrix, as shown in the above expression for the kinetic energy, is readily identified. Note, m - U, ra - w\ and m - Wg are masses associated t e t t S with a translational motion in cartesian coordinates, deployment and elastic deformation, respectively. The coupling of deployment with elastic deformation and rotational motion are represented by the mass vectors m; w\ s T and m - p^ , respectively. Similarly one can T t give physical meaning to other elements of the mass matrix. 53 Using a displacement field w with the second-degree nonlinearity in the generalized 1 coordinate 8{ associated with deflection w and retaining terms up to the same degree of 1 nonlinearity in deflection, Eq.(3.3) takes the form (Appendix A) : b =bl -T + Bl T he = h1 Si + iSi W2 T S{ + {S{ + Hi + [ B\X T fc = •>iT Pi+ir pi = pi+[ pii T w\ + wi Wig T = ®1 + ^2 = 1,2,3 B3lS{}, = 1,2,3 H^Si}, = 1,2,3 bi 6i + {6i = 1,2,3 P3i6i}, *.•] + [ ? f Psmi h - Si —T + 6} = 1,2,3 , m = 1,2,-.-.nj ; W lz Si ; —T i + i 6i S i, = U Si } , W Wis = Wl + [ W}t T Si} + [ S{ s W3lm Si } P = I\ + [ {P2}fm Si} + { S* [J ], 8 m 1> 2, • • •, n ; Si ] , ,m = l , 2 , - - . , n j ; ,m = 1,2,3. (3.4) It is apparent that now all the kinetic parameters are expressed in terms of functions that are independent of deflection, linear in deflection, and quadratic in deflection. Thus, for example, center of mass of body Bi with respect to frame Li, which is represented by T? has contributions from three sources: center of mass of the rigid body Bi ( V ), linear contribution due to deflection <$•, ( B\ T Si), and the quadratic term { Si B3 Si }. Similarly, the total moment of inertia of the body is composed of its rigid term I{ , the linear term [ {l2}Tm ^i 3.2 3.2.1 ] a n < ^* n e nonlinear contribution [ [i^/m Si ]• Center of Mass Center of mass of a body in a multiframe system In the previous section the kinetics of a flexible deployable body employed relative translational and rotational motions of the body coordinates L - with respect to a reference t frame L A - Of course, in a multibody system, the orientation of a body is affected by the 54 motion of all other bodies on its direct path. Hence, the relative positions of the bodies (along the direct path), as specified by the coordinate systems used, are of significance. Locations of the local body coordinate systems L% were identified earlier in Chapter 2 in the context of multibody kinematics. Now, any arbitrary motion of a body can be expressed as linear motion of any point on the body and rotation about that point. The point in particular can be the center of mass of the body. Thus, location of the center of mass of body Bi relative to Lj is an important dynamic parameter and it appears in the formulation process frequently. It is, therefore, appropriate to develop the necessary expression for this parameter before analyzing the kinetics of a multibody system. Integrating Eqs.(2.7a), (2.8a) and (2.9a) over the mass of Bi : mi • mi • — S) m i j I 3 a? drm -- S} 3t - S) jai + C{ t ; {a* dm{ = S} {b* = S lj Jo,- + C} t ; (3.5a) (3.56) J mi I >Na* dmi = S) jNV = S) jfi + d 6* ; (3.5c) J TO.- where (Figure 3-2) bl center of mass of body Bi with respect to frame L,- projected on frame Lk ; bl center of mass of body Bi relative to frame Lj (j = k) ; bl projection of 3b l on the inertial frame, LN [k = —1) ; V center of mass of body Bi with respect to L - (i = j = k) . k ] 3 N 3.2.2 t Center of mass of a multibody system As for a single body, the reference point for a multibody system can either be a specified point or the system center of mass. Here the expression for the position vector of the center of mass of a multibody system is obtained. Figure 3-2 Center of mass of a generic body in a multiframe system 56 Taking the mass moment of the whole system about the origin of the inertial frame, / f^drrii = ma0 . (3.6) From Eqs.(2.1l6) and (2.106), n r = a~o + Co a>\ + l Sy Cj-i a~j + Ct- £ , y=2 —% = a-o + Co ai + ir,- + C{ £ . Substituting from above into Eq.(3.6) and using Eqs.(3.2) and (3.5), ft ft ft m C0 ai + £ mi {if, + C-6*} = m C0 a"i + £ m - ^6* = ^ m, ^6* = 0 . t=i t=i t=i Thus the expression for the cm. has the form t t ai = 1 n Cl £ — m{ Nb { , (3.7) (3.8a) i=\ l = n Yl ™i 0? • i=i (3-86) Note, a"i represents position of the origin Oi of frame L\ on the main body B\ with respect to frame Lo attached to the cm. (Figure 2-6). 3.2.3 Velocity of the center of mass Using Eqs.(3.7) and (2.106), n n m Co Si = { ~ £ M t » N * > 6 n Cy_i ay + Cf &*' } = ~ X > - { £ 5 j y=2 t=i = n t=i = • + 1 Cy_x Oy } , y=2 (3.9) t=i The arrow over ay is used to distinguish this expression from equations like (2.96) and 57 (2.11c). Differentiating Eq.(3.9) and noting that n m- m = and t m= 0, i=0 gives t=i [ y=i J y=2 di+i = b . (3.10) Now from Eq.(2.13), £ Cy_ ay = £ y=i £ x j=i 5/- 1 [ C, a,;- C } . , ] ay , 1 /=o = £ Cr { £ C, a;!" 1 5/'- 1 Cy_! dy } , y=i i=o Also from Eqs.(3.5) and (2.10a), { £ SJ 5/'" Cy_ dy } = { £ SJ Sf- 1 +1 X 3=1 Cy_ dy + C V } , 1 x t j=l = S[ Sf ,r< + Ci t = S[ Sj b { N . Hence n £ 3=1 Sj+* Cy_, dy = J2Cl 1=0 n = E 1=1 s i Z'r 1 C l Cf 1 c { Sf i r S l t f } , » f=i Substituting from above into Eq.(3.10), »=1 [ m = ~ m °o £ I i=i l 3=1 ^ JTS* + mi £ y=i 3=2 Sj ( Cy_ ST*" )* 1 x d t + 1 = V . 58 +» E i m (3.11) o-j + mi Ci b S j=2 However, with the system center of mass as the reference point, (3.12a) ui = 0 , as now the system rotation is described by UJQ. Also, from Eq.(3.2), 6 =lihe + (3.126) H tsSi Substituting from Eq.(3.12) into Eq.(3.11) gives Co" E *i = ~ \ +< m i=i Eho * / i Sj ( C _ ! W r j y l )** f N y=o E ^° hi j +m» S Cj-! tj + mi d { l t { e + H\ h } \ , which can be written in a matrix form as f rnP \ e w s 7*1 m J V ± C T b N 3 mjCjhg m?8JQ8JICJ—\ mjCjHg — 6JQ8J\C\L^C^ j< 6 (3.13) Here : re m = j ]P i=i Sj mi ; (3.14a) 59 n C ^ m ^ ^ C 1 5 t'=l n (3.146) | 0, \ 1, for j = k ; for j ^k . Note, the velocity of the center of mass of the system is a function of the translation, rotation, deformation and deployment velocities of the constituent bodies. 3.3 3.3.1 Kinetic Energy of a Body in a Multiframe System General form The kinetic energy, Ti, of body Bi is defined as *t *t where r is the inertial velocity of a generic mass element dmi. Substituting for r from Eq.(2.16) in the above expression gives +( n n ( c uj)x V ) • ( J^si (c U ) X V ) ) y j=0 k k dmi•i > k=0 (3.15) where : 60 y=o y=o = 3 ,/m t / 'f* x r dm,- ; J m; J c f EX>y T * / * *k ££ = y 5{ / 5 ( ( Cy V ) • ( (C u7y ) X k U k )X V ) dm, , j=0 k=Q ft ft = EE H{ k = S} S Cf[f l k [ >V Jm- = 5; S C f \ f l k y ^ 5 5 • V [ >Y f k y T c 17 iT k f [ / [,y • f c f i u r >r { iT } drm ] dmi\C k 1 «] *^; dm k f i 3 y T c , Cfc ; U — is a 3 x 3 unit matrix. (3.16) In arriving at the expressions in Eq. (3.16), the following indentities between vectors a, b, c, d are used : a • (b x c) = b • (c x a) ; (a x 6) • (c x d) = (a • c)(b • d) — (a • d)(6 • c) ; a • b U — a b = —b a = b a = b a T . T Note, the three terms on the r.h.s. of the kinetic energy expression (Eq. 3.15) are associated with the translational, coupled translational-rotational, and rotational motions, respectively. Recognizing that the kinetic energy expression is composed of velocity-square terms, Eq.(3.15) can be written as EE^ ^^ ) > T Ti = \ (vfEiVi + 2 £ c J \ y r C ^ + y=o y=ofc=o J or in the matrix form 1 I Vi 2 I ufi Ei kT G1 Vi U k (3.17) 61 Here, v~i is an array of the translational-type velocities ( rigid body translation, deployment, and elastic deflection ) corresponding to body f?,-. The mass matrix, associated with the kinetic energy in the above equation, is basically composed of three submatrices : E{, G k and H jk. Ei represents the mass matrix associated with the translational motion. The second matrix, G k, is the first mass moment of body 2?,- about the body-fixed reference frame Lfc. It is associated with the coupling between the translational and rotational motions. The third matrix, H jk, is the second mass moment of Bi with reference to frames Lj and Lk- It is apparent that when j = k, H? k represents the conventional mass moment of inertia diadic. For j ^ k, it corresponds to the 'mixed' mass moments of inertia matrix. 3.3.2 Translational component of kinetic energy From Eqs.(3.15), (3.16) and (3.17), it follows that e,- = vf Ei Vi = j F* . f 1 drrii . J 171; Substituting for f* from Eq.(2.16d) and using Eq.(3.2) « „ « f —i e,- = m,r - . r,- + 2r,- . C / t = rm fj t J m; ft + 2 rm rj d f f drrii + / f J m,- _i» . £ dm,- , {I Ve + Hi T ti) + rm (I 2 w\e + 2t T vTeS t{ + o\ W lss ti) , which has the matrix form < _Sym. <s T W lss | f | • (- ) 3 18 J Ui J It should be recognized that Eq.(3.18) shows components of the translational velocity vector, Vi, and the translational mass matrix, Ei. The diagonal submatrices ( U, uJ^ , WgS ) e represent masses associated with the velocity components, while the off-diagonal submatrices correspond to the coupling terms and may be designated as 'mixed-masses'. 62 3.3.3 Coupled translational-rotational component of kinetic energy From Eqs.(3.15), (3.16) and (3.17), the moment of momentum of body Bi about the hinge at Oy is given by, f. = G{ Vi = SJ Cf I > V x f drrn • J rt%i Substituting for Sj 'Y and r* from Eqs.(2.9a) and (2.16d), g{ = Sj Cjl I J {jfi + d f) x Vi drrn + jfi x C - I f t mi J drrn + C f t m^ J T *t drm \ . m,- In light of Eq.(3.2), g{ = rm Sj Cf {{yr- + C t} t { x f,- + yf- x d{l ti + H e t l T s tt-} + d{l % + , which can be written as g\ = rm Sj Cl [ &> CT ?J G/ } | % J , { (3.19a) where : G'J = Cj [ yfi + cTb* } d , = [jo - + 6 ']=j6 '; 1 (3.196) t t r y e = C f { yf,- Ci V + d pi } , e = { )ai V + pi } ; (3.19c) e G /= { CT [ jfi Ci Hl = [y Hi T ai T + dPl + P }. iT s T }, (3.19a-) Thus components of the first mass moments of body Bi about Oy take the form as shown in Eq.(3.19). The expressions in Eqs.(3.196), (3.19c), and (3.19d) represent the body as a point mass, a deploying member, and a flexible body, respectively. Thus, in the absence of deployment, the second term vanishes, while for a rigid body the third term has no 63 contribution. Note that the r.h.s. of Eqs.(3.196, c, d) consists of two terms: the second term represents the mass moment of the body about its own coordinates L,- while the first term is the first mass moment of body Bi about origin Oj of frame Lj as discussed in Section 3.1 . 3.3.4 Rotational kinetic energy Earlier Eqs.(3.15), (3.16), and (3.17) gave the expression for the moment of inertia matrix as Hi" = SJ Si Cj[ f [ f f j {k ] drm C , iT k = Sj Si Cf [ f [ V ' . V U - V >V ] drm } C J T k . (3.20) J m( The second term on the right hand side is referred to here as the pseudo-moment of inertia. In the following, an expession for the pseudo-moment of inertia is first developed, followed by the expression for the matrix of inertia. Substituting from Eq.(2.9a), the pseudo-moment of inertia of body Bi can be expressed as f S} S V >V dmi = Jmi Jmi l k T j {Si ¥i + d t} {S} jfi + d f } T k dmi , = mi Sj Si fi jfi + Si fi {d j t dmi} J mi T T k k + S) {d I t dmi) jTi + d f f f J mi J mi T T dm Cf , { which with Eq.(3.2) takes the form Sj S I l k V >V drm = mi [ Si Sj r r , + S* r'i t r r k { y kk T C?' + SJ d t r - + d J Cf } T y t { J mi where J - = / . £* t m dm,-. Using the relation obtained earlier, 5*-N V = Sj jfi + Ci V ( Eq.3.5c), the pseudo-inertia for body Bi with refrence to coordinate frames Lj and L k 64 can be written as SJ Si V I >r iTdmi SJ S{ [c [J - - Sj b?\Cj = rm t + t k N P j N V Thus the expression for the inertia matrix H? , as given in Eq.(3.17), takes the form k Hj k = Sj Si Cf [ d = Sj Si C{ [U - [ /,• - m t mi b { b iT V b iT } Cf + + rm i~b { m i N V k N b {T] C , k fb {T Cl (3.21) With the individual contribution due to translation (Eq.3.20), rotation (Eq.3.21) and their coupling (Eq.3.19) in hand, one can now proceed to combine their effects using equation (3.17), T Ti < 1 Ci Tie Ci Ke . si Hg T d si Ks fb iT f kT ( — c{ Ti ci U < > si Si G\ t -±-Hi Sym. \ ci • Si k u k A word concerning the notation used here would be appropriate : u; with a superscript, as in U , indicates that the quantity forms an array of all the angular velocity 3 vectors. For example, ZJ 3T = (UQ ZjJ .... tJ^). Subscript appearing in ri,£,-, and 6i implies that the parameter is associated with body Bi. From Eq.(2.16e), r,- = X2y=o S)Cj-i<Lj , hence the above equation takes the form a Ti = sisici'] j k k —1 ii ^c? j x _ 1 £ e sicfm J t w eS Wis Sym. 7 , o sjsicr w Tci x a Si(K Tik~ai T + p /)Ci i Si[Hii~ai T + PJ]Ci ±H? k ii Si (3.22) 65 To recapitulate, this is an expression for kinetic energy of a flexible deployable body in a multiframe system, where the velocities are associated with the motion of each coordinate frame relative to its adjacent one on the direct path to body Bi. 3.4 Kinetic Energy of a Multibody System The kinetic energy of a multibody system is the sum of the kinetic energies of the constituent bodies. Hence, summing over the entire system consisting of n bodies Bi (i = 1,2, ...,n), the total kinetic energy expression can be written as ( zj \ a m C Z\ 3k 3 k m S*Ci- xh) m S^C{- H l V mj6 w{ jk rnj6j w' e k ( •k \ a C^V^C{ kT k k mjS g C 3 eS 3kT k '& 3 k > < 6 ™j6j W m,*iGi C' 3 k kT s ik Sym. k •k > 5 —k where : 1=1 Lf = J2m S}Sici b C}; k i i »=i n W k = E Sj Si (3.23) Hl . k t=i Here, m , LJ and H ' k jfc 3 k are the masses, first mass moments and the second mass mo- ments, respectively, associated with the reference frames Lj and L . When j ^ k, these k quantities may be recognized as mixed masses, first mass moments and second mass moments, respectively. Other terms in Eq.(3.23) were defined earlier. As before, a superscript associated with a velocity symbol indicates an array ranging over all relevant bodies. For example, I is an array of deployment velocities associated with all deploying bodies in the 3 system. Thus, each term in Eq.(3.23) is a matrix of size j,k ( j, rows; k, columns ). Each 66 element in turn may be a scalar, a vector or a matrix. The state of the system, as described so far, is defined by relative coordinates emanating from the main point, located on a predefined main body, together with the relative position of frame L 0 with respect to the inertial frame, represented by (^OJJ/OJ^O) a and C . In practice, the coordinate system x ,yo,z 0 0 0 0 is either fixed to the system center of mass or to a specified point on a body, which then becomes the main body. When the system reference point is specified, i.e., when the refrence frame xo,yo,zo is located at a definite point, a~i and C® are specified and a~x = cJi = 0. In this case the kinetic energy expression, as given in Eq.(3.23), can be rewritten as ( lj a 1 . T r hihiM kClz{ mJ^Ct 1^ mJ^Ci-'H^ hAiC^ 1 L? T C>k -] > ( a-k \ < u 6 SikWit m,-$klSiG> kTCi mi Sym. &ji&kiH' k s• k —fc (3.24) Here, T refers to the kinetic energy of a multibody system when its reference point is B specified. In above 6j is equal to zero or one depending on j = fc or j ^ k, respectively. k In the case where the coordinates xo,yo, ZQ are located at the system center of mass, a"i and a"i become functions of the generalized coordiates. Moreover, C° is a specified quantity and TJJi = 0. In this case the kinetic energy, T, takes the form T = T +T = c s (M + M ) x-, e (3.25) a where T is the kinetic energy of the system excluding the motion of the center of mass 3 oi, as defined by Eq.(3.24), and T is the kinetic energy associated with the motion a of c x the system center of mass. The mass matrices associated with T and T are denoted as s M a c and M , respectively. The kinetic energy contributed by the motion of the center of c mass can be obtained using the the components associated with a"i in Eq.(3.23) , 67 a T = lmii Z +i%CZ[ C h T e 1 1 mk mCH m 4i^_i k k A e k 6 C L\ Cl\< kT > kT k kl 1 6 (2.26) The first term represents the kinetic energy of the system treated as a point mass and located at the center of mass. The second term denotes the coupling of the center of mass translational velocity with other coordinates in the system. Substituting in Eq.(3.26) for a"i from Eq.(3.13) and noting that the first mass moment of the system about its center of mass equals zero, L\° = 0, gives i r>r N° m,-hi Cj ' rn \ 1 ti ( m' £fc T m^-o^-iCj-i [],? m C t k k e m^fco^iCfc-! mCH k k k h»h±C L\ kT 1 Cl)< a* s u' 3 w J 1 ( m> pT (m ti 1 2m i m ^ 5yiCT_ J 0 1 [0 2m C t k k e 2mH C -x kl k 2m C H k k k 28 6 iC L Cl]< 1 kT k0 k 1 1 6 The mass matrix M associated with T is obtained from the above equation as c c k a* 68 Mr — [-jft mCt k k m hi(2 k e - 4o)C _! mCH koki^L ^ k k fc C^ 1 k 1 m,H 3/Cf The mass matix M of the system, M = Ms + Mc, can now be developed as 0 0 0 mAiSih'J Ci_ M 1 hikxM Ci~_\ m^Stf™ m^Ci-'H^ k = 0 m^wis r 3 i ^ i C i - m -6 W> ik 3 C mj ^ 1 8S G C 3 s kl 3kT k 3 k k Sym. (cont.) 0 rn 6 tf Ci m]khe TC 3khe mfatf/CU T k ko 1 m fc-i 6 h C LilfcT m h\ CiH T kT 3 T ik kimj m'"*$ *,Ai(l + S o)C Z\ miS C - H 3 yo k 3 k jl 1 kT k 3 3 k ^i^i^om'Cj'- ^* 1 m HC H ]k 3 e m -6 6 H C L kT 3T k 3 3 kl k0 k r kT k -hohi5™kxL) C L . 3 Sym. 3 k kT k With each mass my in the system nondimensionalized with respect to the system total mass m, a rather simple form of the kinetic energy expression can be developed as follows: 1 •T T = - m x" 2 where * - 1 Mx~=-mx~ T[Mi + 2 = (m k l k U k 6 k 57*) A TBA}x~, ; (3.27) 69 0 0 S mjwi jk Mi e 0 0 0 6j mjw{ 0 0 k m,-6,-kWis 0 T s m,-6klslgi kTC>k 0 m,-6klS>kG>s kTC>k .sym. HI and the matrix B is given by 0 -mikC 3k -m kSkoC Jk_1 m,4(Sj-m*)q(-i 0 -mjkCl —rrijkCl 0 -mj-SnSkxC' —mjSnSkoSkiCl SiiSjoSkoSjiSkiCl _ sym. The mass matrix thus obtained is applicable to both the cases, i.e., when the reference point is the center of mass or a specified point. In the latter case, only the underlined terms in the B matrix are to be included. 3.5 Effect of Hinge Flexibility It is necessary to develop a transformation matrix that can account for two factors : (a) contribution of the elasticity and deployment of the system members to the translation and rotation of other bodies; (b) translational and rotational motion of the system reference 70 point relative to some intermediate coordinate system. The array of system coordinates can be written as m" fk £fc ao . - . • fc 6 . or* . S* • fc 6 w 0 57* where the translational and rotational coordinates of the overall system reference point ( ao, aJo ) are indicated separately. Hence, for a and U , superscript on the left hand side varies from 0,1,2, • • •, n, while on the right hand side it ranges over 1,2, • • •, n. For other coordinates, fc = 1,2, • • •, n on both the sides. The translational velocity of the system reference point, i.e., origin Oo (Figure 3-3) can be written as a —Eu 0 0 0 + E r u , r where : Eo UQ = translational velocity of the system reference point relative to the intermediate curvilinear coordinates ; E u r r = translational velocity of the origin of intermediate curvilinear coordinate system in an inertial curvilinear frame. Also the rotational motion of the system reference frame is expressed as aJ + C° a7 , where 0 r uJo is now the system angular velocity relative to the intermediate frame. The angular velocity of the intermediate frame relative to the inertial frame is denoted by uJ . Note r that the angular velocities UQ and u7 are expressed in terms of the system coordinates. r 71 Figure 3-3 Translational motion of system coordinates relative to an intermediate frame. 72 Using these transformations for the translation and rotation of frame LQ together with Eqs.(2.18) and (2-26), the array of the system velocities can be written as / m • it m l l K k k a~o U k E 0 >=< S jk ii + E 0 R - u ° + Sj 3 k 8 r ir r V Wk + Sjk F T k f 6 U k + Sjk B U k k G 6 k hence • lk a u Sjk a0 Er fe Sjk rfk & T «o EQ Sjk Sjk Rj • fc Fk T (3.28) Sjk S WT w0 Sjk Bk Gk Wo Sjk Wk. '0 •J I 9 ) where rf , Rj, Gk and Bk = C £ , are as defined in section 2.5, and : k j,k = 1,2,... ,n, n = number of bodies in the system ; W 0 = system angular velocity relative to the intermediate coordinate system ; W r 0 r = angular velocity of the intermediate coordinate relative to the inertial frame; Wk Ok = angular velocity of body B k relative to the preceeding body on its direct 73 path due to rotational degrees of freedom of hinge Ok ; Sjk = unity if Bj is adjacent to B on its direct path, otherwise it is zero. k With this transformation the old state parameters, as defined by Eq.(3.27), are related to the new parameters as x = Ty v , where T y represents the r.h.s. of Eq.(3.28). Substituting from above into Eq.(3.27), the v system kinetic energy takes the form T = - f T ? M T l v y - , = i i f T„ [Mi + A T T 3.6 B A}T y-. v (3.29) Holonomic Constraints The division of the system into domains and subdomains, as discussed in Section 2.1, result in constraint relations which are generally holonomic. Nonholonomic constraints are treated later during the development of the equations of motion. Here the attention is focused on the holonomic constraints, which have the form C q — 0, where C is a matrix of constraints and q is a vector of the dependent coordinates. From such constraint relations one can develop a transformation matrix, T , which maps the dependent coordinates , q, c into independent ones, q, as q — T q. In the following, three methods of obtaining the c transformation matrix T are briefly discussed. c (i) The transformation matrix necessary to identify independent coordinates can be obtained by solving the constraint relations as simultaneous equations. The solution leads to the selection of certain of the coordinates as independent and expressions for the remaining coordinates in terms of those which have been selected to be independent. To clarify the method, the constraint equation, C q = 0, is written in the form Aq + Bq = 0 , where q are the selected independent coordinates and q are the remaining dependent 74 coordinates of the set q, which are not included in q. The dimension of q is such that A is a square nonsingular matrix with distinct columns of C corresponding to q. Similarly matrix B comprises of those columns of C not included in A but corresponding to q. The constraint equation can now have the form q = —A~ Bq. 1 •A~ B U Therefore, the required 1 transformation matrix T = c Note, the transformation matrix T is not unique c in form. (ii) The alternate approach is based on the the zero eigenvalues theorm [54]. With the same notations as above, let the symmetric matrix E be defined by E = C C. As seen T before, the solution of the constraint equations may be expressed in the form q = T q. c According to the theorm, T is a matrix whose columns are those of any modal matrix of c E corresponding to its zero eigenvalues, and q is a column matrix compatible with T . c (iii) In the finite element method of structural analysis the coordinates of a free body element may be chosen as displacements and rotations at node points of structural elements. In such analyses, the equations of constraints require matching of the displacements and rotations at the nodes. A set of independent coordinates is determined simply through the use of a single symbol for each set of displacements and rotations in conjunction with the equality constraints. This is the basis of the now widely used procedure of superposing stiffness matrices or mass matrices of structural elements to determine a corresponding effective representation for the entire structure. 3.7 Potential Energy The system's potential energy comprises of two contributions : position in the inverse square gravitational field (U ), and strain energy stored by the structure due to its g flexibility (V). Each of the contributions are evaluated in some detail here. 75 3.7.1 Gravitational potential energy The gravitational potential energy, U , of a multibody system, located in an inverse g square gravitational field, can be expressed as t =i j=i Jm i Jrn j i T where : G universal gravitational constant, 6.67 x 1 0 n number of bodies in the system ; p number of primary bodies ; m- mass of body Bi ; rrij mass of the primary body j ; t rj N.m /kg 2 - 1 1 2 ; distance from drrii to drrij . l In the present development, only one celestial primary body, the Earth, is considered, hence j — 1. Also, the Earth is taken to be spherical and homogeneous, i.e., it can be represented as a point mass. Therefore, the gravitational potential energy expression reduces to u * = E / ' i=i i Jm ( - °) 3 3 where : fi = G m = 3.986 x 10 Nm /kg 14 2 e m = 5.97 x 10 e 24 ; kg ; r = distance from center of the Earth to drrii • % From Eq.(2.12b), r =r + 0 r =a + 0 r . With the inertial frame located at the center of the Earth, the position of the system 76 reference point, a~o, can be written as a~o = ao I , where ao is the magnitude of the position vector a"o and £ is a unit vector in its direction. The position of a differential mass relative to the reference point as projected on the inertial frame, ° r * , can be expressed using Eqs.(2.9a) and (2.10c) as : n ° r * = fi + C° £ ; o¥i = £ 0 5/ C / _ i ai . l=i Now („1 i 0°-*' — i 0-t 0-i<i — i (a + 2 r . ao+ r . r ) ?, 0 - ( 1 + x) ao 7, where a; ~2 r •(2ao + a0 r ) . Hence 1 1 , — = —(1 r* a0V _ - a ^ " a l 1 1.3 2 x -\ x2 2 2.4 1.2.3 o x 3 + ...) 2.4.6 ' ^ (rr °=* _i_ ^ = " 0 ( G 0 - " + 2 _i_ l r r ) + - 3 2^V I r — °=i\2 , t— 0-^'wO-t O-t'-v , ^ / 0 - t ( a °- 0 + i oa r )( r . r ) +-( 0-i\2 \ , r • r ) J+ Substituting from above into Eq.(3.30) , £2. ± j Vom, + 3/i / 2a; ^ ao where + ii(ZI ) + ap 0 f ( V . V - 3(1 V ) ) 2 <V)«m.- - ^ ± j (V.V) dm . + ... 2 t * ( - t r J E T 0 + 3£.i7°.£) - ^Zh° 2ap 2aQ + ..., (3.31) 77 £0 = J27=i fm °r*drrii Y17=i * %^ > * ' m * n e n r s m a s s moment of the system (as a point mass) about the origin Oo of frame LQ as projected on the inertial frame ; N l? = center of mass of body Bi in frame Lo as projected on the inertial frame ; H° = second mass moment of the system about the reference point as projected on the inertial frame ; 7i° = X^iLi ~ Er=i drrii ; third mass moment of the system about Im . the reference point as projected on the inertial frame. The first term in Eq.(3.3l) represents the familiar inverse square force field contribution with the system represented by a point mass at the reference point. The second term is due to the torque as a result of the reference point not coinciding with the system center of mass. The third term is proportional to the gravity gradient evaluated at the reference point while the final term accounts for a higher degree of nonlinearity in the gravitational field. Usually only the first three terms in Eq.(3.31) are included in the analysis. This is adequate for the case where the spacecraft is small compared to its distance from the center of the Earth, i.e., ^f / a | 1 0 << 1 . In the above derivation an additional higher order term, purposely retained, required consideration of the third mass moment of the spacecraft. This term may be important when the spacecraft is very large or its moments of inertia are nearly equal, as pointed out by Meirovitch [55]. The expression of the third mass moment of body Bi can be further developed as follows. From Eq.(2-9a), °f* can be written as V = i.e., 7i + Cit , 0 h°{= I J V °r iT V drrn , m,- = mi 0ri ofi T Ci l Jmi on + [2 0¥i {2ff T 0 r, T + r, 0 + f Tfu}dmiori T or~i ] Ci | * dm,+ + Ci f Jmi f fT f drrn Cj , 78 = rrii r,- r , 0 0 T 0 r,- + [ 3 tr. R % - d 2 R l [ 3 tr. f - 2 P } C f 0r< + ??r?dm.-. (3.32) —\ Here definition of the center of mass b and mass moment of inertia P from Eq.(3.3) have been introduced. R l is the second mass moment of inertia of body Bi represented as a point mass located at the origin O,- of the body frame L,-. In Eq.(3.32), the last term on the r.h.s. represents the third mass moment of body Bi about its own body coordinates, while the first three terms transform it to the reference point. The last term in Eq.(3.32) can be developed further as follows. From Section A . l with a displacement field uT , which 1 is linear in deflection, the position of a point in body B,- with respect to the local frame Li can be written as, T = r + $ i T it. Substituting from above, the expression for the third mass moment about the body coordinates takes the form m = 1,2,3 ; / = 1,2, • • •, nsi . Here the two terms, the first and the last on the r.h.s., represent contributions due to the body treated as rigid and flexible, respectively. The remaining terms are due to the interactions between them. The above equation together with Eq.(3.32) define the dynamical quantities required to account for the additional nonlinear gravity gradient term. 79 The derivation presented here is new in that it extends previous formulation of the gravitational potential energy for a singel body to a multibody system consisting of flexible deployable bodies that may be translating and rotating. 3.7.2 Nonlinear strain Here we consider problems with geometric nonlinearities only, i.e., the angle of rotation can be large with strains not exceeding the limit of proportionality. In the case of a thin beam or a thin shell we assume that every point remains after the deformation on the same perpendicular to the beam axis or the shell middle surface as before and that the distance of every point from the beam axis or the shell middle surface remains unchanged by the deformation. By expanding of the displacement field as a power series about the beam axis or the shell middle surface and invoking the above stated assumptions one can obtain a displacement field which is fully defined in terms of the deformation of the beam axis or the shell middle surface [51]. Use of a displacement field which is formed from the linear terms only in this expansion, together with the complete nonlinear strain energy, results in equations of motion which are adequate for large elastic rotations and account for axial foreshortening due to transverse vibrations [56]. Although this method is not new, its use is limited due to the involved algebra [57]. Actually, the explicit use of the complete nonlinear strain-displacement relations result in terms up to the fourth degree in the strain energy expression and terms up to the third degree in the stiffness matrix. In this and the next section a procedure is introduced for developing an algorithm to calculate the stiffness matrix from the complete nonlinear strain energy. The expressions for the nonlinear strains in terms of the linear strains and elastic rotations are given in Appendix B. These relations can be written in a matrix form as U) = { « } + { <? E{ e } + { e E{ w } + { c T e T e E' u } , 3 e j = 1,2,..., 6 , (3.33) where : {e} = ( en ^22 £33 ^23 £31 £12 ) ; {e} = ( en e 2 e e3 e i wf = ( W i W W ) . T e 2 33 e2 e3 2 3 e 12 ); Here w , ejk, eyjt, j,k = 1,2,3 are elastic rotations, linear and nonlinear strains, reej 80 spectively. E^E^E^ are described in detail in Appendix B. Note, the nonlinear strain expression in Eq.(3.33) is composed of four terms : the first is linear while the other three are nonlinear. The second and the fourth terms are quadratic in the linear strain and elastic rotation, respectively, while the third term represents a coupling between the two. The linear strain-displacement and elastic rotation-displacement relations can be written using operator matrices L and L e w (Appendix B) as: {e} = [L ] {w} ; {u7} = [L ] {w} . e (3.34) u e From Section A . l , a linear elastic dispalcement vector, w, can be expressed in terms of admissible functions as {w} = $ 6. (3.35) T Therefore, the linear strain and rotation vectors can be written in terms of the shape functions and generalized coordinates as : {e} = L w = L $ 6 = 5 ° 6 , B° = L $ ; (3.36a) {w} = L„ w = L„ $o" = G°6 , G° = L„ $ . (3.366) e e e On substituting from Eq.(3.36) in Eq.(3.33), the nonlinear strain vector is obtained in terms of the generalized coordinates as {e} = { B° 6 } + { 6 T B{6 = { B° 6 } + { S B T }+{ B{6 f }+{ 6} , 3 ! B{ 6 } , T (3.37) where : B = B{ + B + B ; 1 3 B{ = B 0T E G°; 3 2 B{ = B 3 2 3 B =G 3 0T 3 E\ B° ; 0T E G° . 3 3 Note, contribution to the nonlinear strain comes from two sources : linear as well as quadratic generalized deflection. It is apparant that the nonlinear term contains contributions from the linear strain, elastic rotation and their interactions. 81 3.7.3 The strain energy and stiffness matrix The strain energy can be written as V = i J^a T edr = ^ J e C e dr , T (3.38) where : T = body volume ; a — stress vector ; C = matrix of material constants . Introducing the expression of the nonlinear strain in terms of the generalized coordinates from Eq.(3.37), the integrand in Eq.(3.38) can be put in the following form , e T C e = {B°8 = 6 T + 8~ B T [B 0 T 8} j [C] { B° 8 + 8~ B 3 8 } T T C B° + B 0 T C \6 B j] + \B iT 8} C B° + [B T 6} \C] [6~ B ] ] 6 . T jT k Substituting from above in Eq.(3.38), the strain energy expression can be rewritten as 1 -T Vi + 6 V + V T V = -6 l 2 2 6+6 V lT T 2 lm z 6] 6, l,m = 1,2,... ,n ; (3.39) 6 where ng is the number of generalized coordinates associated with the deflection of the body, i.e., the number of elements in the vector of deflection 8. Note that only the symmetric parts of V\ and V3 lm are to be considered. Denoting the s-element of 8 as 6 : 8 = j^B 0T V l 6 i.e., V^k) = J2 C B° dr , 6 * X (' ) / C / m Z=lm=l [8 V }= J B T 2 (f; l C [ ) = (X X s v2 i(s,r) /r 8 s=l m=l 6 i.e., 0 T J j=l C °( l>fi dr ; B°(m,k) B t 8 B T ] j dr , / * ° K O (X^ KJ") J 8=1 t V2 (s,r) = X X m = l j=l C ( '-7) J/ t ° ( > ) m B m 1 BJ ( > ) s r '( > ) s r hr d T ) . r 6 l BJ d r 5 . 82 [6 S n T V S] = J [ B lm 3 S 6 n £EW («.r)) B 6 T , k S *8 - n (EE (^) /(5> * (*.0)i (X> *('.»0)*m*), ' = l B 6}[C}[ 6 B ]dr j T C y fc=i \y=i = lfc=l r=l.= l 6 fl y J 8=1 r r=l ' m 6 i.e., y=i fc=i Note the notations used in expressing an element in a matrix. For example: • is the (/, r) element in the matrix [6 ( X)s=i S V (s,r) ) l 3 2 lr V ] ; 2 • ( E r = i E . = i M r V 3 ' m ( a , r ) ) , m is the (/,m) element in the matrix [ ~f V lm 3 Also note that Vi,V are square symmetric matrices, while V lm 3 l 2 6 ]. is, in general, a square asymmetric matrix. The stiffness matrix can now be obtained by differentiating the strain energy expression (Eq.3.39), dV (3.40) [K +K +K ] 6 =K 6, x d6 2 3 where K (= K\ + K + K3) is a nonlinear stiffness matrix obtained as follows. 2 From Eqs.(3.39) and (3.40) : KiS= --^=(6 V 6) 1 = \[V T 1 1 V }6, T 1 (3.41a) Ki = V ; i.e., x K 2 6=\ A( f [t V ]6 + f l 2 = \{2[6 V }6 T i.e., lT 2 ]o], [i K ^ 3 [V} 6}6) + 2[V 6]6 l 2 1 5 T + , 2[V 6}6}, l 2 / = 1,2,-.., n s ; (3.416) * = jj= ( * [ * Vi 6 } 6 ) , T -T S i.e. + K* = 2 [V ± 6 lm 3 T T m +V lmT 3 Vl 6 m , } 6 6, l,m = 1,2,-- -,n s (3.41c) 83 Note that K\ and K are symmetric square matrices while K 3 2 is an asymmetric square matrix. To summarize the formulation of the stiffness matrix : • choose shape functions $ to describe elastic deformations of a body ; • using Eqs.(3.36a) and (3.366) calculate B° and G° ; • with Eq.(3.37), calculate any or all of B{, B and B , depending on the problem, 2 3 and hence B ; J • Using Eq.(3.39) calculate V V and V ; l lt • K\,K 2 2 lm 3 and K , are now determined from Eq.(3.41) once 8 is known at a given 3 instant. The attractive feature of the algorithm lies in the fact that, with the shape functions identified, it requires evaluation of integrals in Eq.(3.39) only once. 4. E Q U A T I O N S O F M O T I O N 4.1 Preliminary Remarks With the kinetic and potential energy expressions in hand one can approach the dynamical problem using the celebrated fundamental equation, due to Lagrange, / d,dT, dT „ £ Ww)~ir .- °) Q \ , c % q = 0 - , (41) where 8q are virtual displacements associated with the generalized coordinates q {s = s s 1,2,... ,n). The kinetic energy T and the generalized forces Q are, in general, functions s of all the system coordinates (q ), velocities (q ) and time (i). g 3 The generalized velocities, q , satisfy the linearly independent system of constraint s equations B q-+b T = 0, (4.2) where B is an / x n matrix of constraints and b is an /—dimensional vector. Therefore, among all the q , q which satisfy the constraint relations we seek those which also satisfy a a the fundamental equations of motion. The virtual displacements 6q satisfy the set of equations B 6q = 0 . T (4.3) From the linear independence of the constraint equation (4.2), it follows that the set of equations in (4.3) are also linearly independent. As shown by Lagrange, if q and q satisfy the fundamental equation (4.1) as well a 3 as the constraint relations (4.2), then one may write [58] 84 85 where Q and A are the vectors of generalized forces and Lagrange multipliers, respectively. Here, 8q are completely arbitrary, i.e., they are no longer subject to the constraints, hence the necessary and sufficient condition for the above equation to be satisfied is The system of equations given by (4.2) and (4.4) are to be solved simultaneously for the generalized coordinates q and the Lagrange multipliers A. The constraint relations in terms of the virtual displacements, Eq.(4.3), can be rewritten as B i SQi = ~B 8q , 2 2 where 8q is an /—dimensional vector from 8q, and B\ is the associated / x / matrix from x B. The rest of 8q and B form the (n — /)—dimensional vector 8q , and the I x (n — I) 2 matrix B , respectively. The square matrix B\ is nonsingular; otherwise the equations of 2 constraint would not be independent. It follows that the constraint equations can be solved for / virtual displacements 8q as linear combinations of the remaining n — I ones, 1 1 1 Sq = -IB; } x llT [ lT D c- [B Y 6q 2 2 . Denoting the fundamental equation (4.1) can be written as Rf 6q + R% 8q = { -B x 2 B^ R + R 1 2 x As 8q are arbitrary, 2 1 —B -Bj R\ + R = 0 , 2 2 } T 2 8q = 0 2 86 i.e., n n . ! ( d i r 1 dT , d,dT, 8T „ , , The equations in (4.5) together with the constraint relations (4.2) constitute the formulation by embedding of nonholonomic constraints. They form a system of n equations with n unknowns, g, and involve no Lagrange multipliers. Hemami [59] proposed a method to account for nonholonomic contraints, where one premultiplies the equations of motions by a matrix, which is an orthogonal complement to the matrix of constraints. The new, reduced system of equations and the constraint relations form the set of equations to be solved. Since the orthogonal complement to a matrix is not unique, the reduced equations of the system is not unique in form. Consequently, for the reduced system, a variety of representations are possible. The method of embedding of constraints, discussed above, is essentially similar. Note, different representations for the reduced system may be obtained by using different 6gl . Thus the method of embedding of constraints provides the orthogonal complement of the matrix of constraints as [ -B 2 Bf1 U ] , where U is the unit matrix. The equations of motion as given in (4.5) can be obtained by premultiplying Eq.(4.4) with its orthogonal complement eliminating the Lagrange multipliers as follows, [ -B 2 sr 1 u \ In the above equation, the matrix of constraint is decomposed into two submatrices By and JE?2 a s before. 87 4.2 Functions of Interest 4.2.1 Kinetic energy In many dynamical applications, one is required to specify some of the coordinates as functions of time. Since in the present work relative coordinates are used they can be rearranged to express the kinetic energy as T=-x 1 •T M* x , 2 T = -q TMq~ A 1•T = -Q~ = M NT r + J N M3 N q~+ -j A Ma j , Mq-+n Ti[+to, (4.6) $2 + ti + to j where s are the specified velocities. M, Ma and N are masses associated with generalized coordinates, specified coordinates and their coupling, respectively. M* is the mass matrix associated with the system coordinates (generalized and specified) denoted as x. The kinetic energy is composed of quadratic, linear, and independent terms in the system velocities, represented here by t2,ti and to, respectively. 4.2.2 Lagrangian Denoting the system potential energy as Up, the Lagrangian of the system can be written as L = T — Up , 1 •T = - x 2 = ^f M* x — Up , Mii+n T il+{to-Up) . (4.7) 88 When the kinetic energy is only a quadratic function in the system velocities, n and to vanish and the Lagrangian is called natural. 4.2.3 Gerneralized momenta The generalized momenta are defined by (4.8) p= — = M q + n . dq Note that when the kinetic energy is quadratic in the system velocities the generalized momenta have the form p = M q . 4.2.4 Hamiltonian From Eq.(4.8), the generalized velocities can be expressed in terms of the system momenta as q = M~ x {p — n} . Substituting from above into Eq.(4.7), the Lagrangian can be written in terms of the generalized momenta as L = -{p - n} T M " 1 {p-n} + n T M~ {p-n} l + t -U 0 p . Using the above two expressions of the generalized velocities and the Lagrangian the Hamiltonian, H, can be expressed using its definition as H = = p p q--L, T T M - { p - n } - L , l = \{P ~ n} T M' = ti ~ h + U . p 1 {p-n}-t 0 + U , p (4.9) Note that t\ does not appear in Hamilton's function. When the Lagrangian is natural 89 (*i = to — 0), the Hamiltonian is just the system's total energy. 4.2.5 Integrals of motion When the Hamiltonian is not a function of time, it represents a constant of the motion. In addition, the case when the kinetic energy is quadratic in the velocity represents conservation of the total energy. Hamiltonian, when constant, can serve to check the accuracy during the numerical integration. It can also act as the Liapunov function in the stability analysis. If the external forces associated with some generalized coordinates g, known as ignorable coordinates, are zero, then the momenta associated with g, denoted as p , are g constants. Hence p = M g where M g g i} + n = c , (4.10) g is the mass matrix associated with the ignorable coordinates. Note, n = N y g with y denoting the system generalized and specified velocities; while N g g represents the mass matrix, which couples the ignorable coordinates with the system. Here c is an array of constants defined by the system initial conditions. From above, t=M~ 1 {c-fig} (4.11) . Thus with the ignorable coordinates g, the order of the differential equations to be solved is reduced. 4.3 Lagrange Equations for a Multibody System Substituting for Lagrangian in Eq.(4.4), the governing equations of motion can be written as M*x + M* x — Ux T M* x) + U , p x =Q + BX . 90 •T Separating x into generalized and specified velocities q and s, respectively, x •T = (q •T s ), Lagrange's equations can be written as M N N M T q + s a M N N M T l(x s 8 + [M\ ] T x) q 2{x [M*, ]xj T ( U, p q {U , + a p a (4.12) B a The first set of equations describe the system dynamics, which must be solved for the generalized coordinates q, while the second set of equations are for the specified coordinates 5 and, when solved, gives applied forces associated with s. An alternative form of the equations describing the systems dynamics can be developed as follows. Noting that (Eq.4.6) : N s + N s = fT ; 1 1 .T j_T -x M*, x=-q M, qk q + n, qk q+ t , qk Q qk , k = 1,2, • • •, n ; q Eq.(4.12) takes the form M ^ + [ M - n J j ^ - i { t M ^ } + ^ - r - ( L 7 - t ) ^ = Q BX , T g p g+ J 0 g k = 1,2,... ,n ; (4.13) q where n is the number of the generalized coordinates. Noting that : q n q G= E *>ij J + q )1 *~ + >t 5 and W = 3=1 M = Y, «3 M V + Mt; 3=1 [ >*3 M ?J]?={£E > q q m % ),qj)k M }= {t m [M{k,m) } tqj k f} , •y=i = i(EE^' M{]\k), ) qm k } = {q T [M{j,k) ] !qm k q} 91 Note, M(k, m) . and M(j, k) i9j are the (j, m) elements in the matrices [M(k, m)^^ and tqm [M(j, k) ]k, respectively. Eq.(4.13) can be rewritten in the form jqm Mf+[M t + G] q-+^{f T q:} + n t + {U -to),<i = Q + B A , k i p k = 1,2,... ,n , (4.14) q where G = [ ,9fc n 1 - [ ,q 1 > n k r*(m,j) = 2M(fc,m) . i? = M(m,k) . jq * = l,2,---,n, ; - M(j,m) , + M(k,j) - M{j,m) tqk >qm , )qk j,k,m = 1,2, • • • , n q . Here G is a skew symmetric matrix associated with the gyroscopic forces while T , k = k 1,2, • • •, n , are symmetric matrices representing Coriolis type forces. q Lagrange's equation can be expressed in yet another form. Let n denote the number a of specified coordinates and note that : i=i n a n = JV s = i / n/ ; l=i na n,t = Yi s n n XX a i + n a im ^' m s ; 1=1 m=l 1=1 n n a G 1 =X ' I 5 - 1=1 Gi = [ nj, to = s M gfc s a [ "Fwfc ] ] = X G l 1=1 ] - [ nf, qk s= Y ' A; = 1,2, • • •, n ; ], q X *'^ m Mail,™) 5 (4.15) 1=1 m = l where ni is the /-th column of the matrix JV. Substituting from above into Eq.(4.14) gives 92 Mq-+ [E*/[ M +Gi ]] tai T -{ f £} + E ' ^ S k / J=i +E E 1=1 -SiS {n/,5 m M — M (/,m),g} 8 + U ,q = Q + B X . (4.16) p m=l As before, G/ is a skew symmetric matrix (gyroscopic force) associated with the specified coordinate sj. It should be emphasized that Eqs. (4.13), (4.14) and (4.16) governing the generalized coordinates q account for the coupling with the specified coordinates. The equations of motion, or more appropriately force equations, for the specified coordinates are similar in form with q and s interchanged and, of course, the corresponding adjustment in the associated mass matrix, external forces and constraints. Ignorable coordinates can be accounted for as follows. By subdividing the system •T •T •T g velocities a s x = ( g •T s ), accordingly identifying the coefficient submatrices in Eq.(4.12) and following the same procedure as before, the resulting equations (similar to 4.13, 4.14, 4.15) will now have additional terms to account for the separated ignorable coordinates. As for the ignorable coordinates, Eq.(4.13) reduces to M g+ M g where M g g g+ n = g — { M g g+ n } = 0, g is the mass matrix associated with the ignorable coordinates g, and n is a g coupling term associated with the coordinates q and 5. Hence V= ~ g~ M M g i.e., l { g>t V+n , } ; M g t g + rig = c , ^= Mg' 1 { c-fTg } ; where c is a constant to be calculated from initial conditions. The nonconservative external force, Q , associated with the generalized coordinate r 93 q can be derived as r n d r*- t 3 «r = E E / t=l y=i r*i J a Q r where : n = number of bodies in a system ; = j-th component of the external force per unit mass acting on the differential element dm,- ; r*. = j-th component of the position vector f to dm,- i n the inertial frame . 1 The position vector f is defined by Eq.(2.11). Its partial derivative with respect to the 1 generalized coordinate q is given by r n ^tqr { Cl-l = ^E ,q + r 0 / , g +C{, r £ + £,> Ci £ , qr q r j, 1=0 where the partial differentiations involved are listed in Appendix C . 4.3.1 Structure of the linearized equations In some applications, such as control permitting small deviations from the equi- librium state, it may be adequate to use the linearized form of the governing equations of motion. Such a form can be obtained quite readily by neglecting the quadratic terms in the velocity and expanding the terms i n the equation about some nominal equilibrium point : h fii = s { rti + [ ni, t = h{ni q + ^[G l ] q} , + E ]q}; l and sis {ni, m Sm -M {l,m), } s q -- sis {ni, where terms on the r.h.s. m 3m -M (l,m),,+i[Gj, 8 Sm +Ei, Sm -2M (l,m),^] 8 q} ; of the above two expressions are evaluated at the nominal 94 equilibrium point. Here : Gi — [ ™hq } ~ [ ™Tiq ] k Ei = [ ni, +nf , qk Gi , Sm E — [ ni,q k i>a =[ni,q m k ], qk },s >?fc Lam ' [l ],s -[nf,q k skew symmetric matrix ; a symmetric matrix ; n m a J m ], k am a , skew symmetric matrix ; a symmetric matrix ; Using the above equalities i n Eq.(4.16), M f + ^ B j [ M, +Gi ] f + Sl lYh[Gi i=i + Ei] i=i 1 ns ns 1 m=l res res res 1=1 1=1 m=l which can be rewritten as M f + [ M , +G + D ] ^+ t [ i [G , + £ , t t ] + (LT - t ),qq-Q,q p 0 ] f = (t - Up),g~n,t + Q , 0 where I? is a matrix due to damping. Here, again all the coefficient matrices are evaluated at the equilibrium point. Note that M , D, E, U qq and to,qq are symmetric matrices and p> G is skew symmetric. The forcing terms U q, to,q, n,t and Q vanish if the linearization is Pi about an equilibrium point. Note that both the specified coordinates and external forces can lead to asymmetric geometric stiffness matrix. The linearized equations shown above have a standard form which is well studied in the context of stability analysis. Specialized variations of this general form are used in the analyses of various systems (conservative, psedo-conservative, circulatory, gyroscopic conservative, etc.). It is also used i n the development of modal control design, where the set of equations is transformed into canonical form [60-61]. This procedure leads to computational advantages as well as results i n a better physical implementation of 95 decoupled control through a repeated pattern of control units. 4.4 Equations of Motion in Terms of the Hamiltonian Hamilton's canonical equations are dH - +Q +BX, dq ' ^ ' " ' a q . q= q * 3H — dp Using the Hamiltonian (Eq.4.9) and noting that : MM' 1 = 1; q= M M — oq k — = —M dq k k —— M oq 1 k {p - n} ; 1 p = -\{{p - n} [M-\ T qk + { {to-U ), +Q p = \{{V-n} q M' T 0 T p q [M, qk q } {p - n}} + [n , ] M T k +B \ q + { (t -U ), +Q = i{t = 0; + — M oq 1 [M, } qk q q 1 {p-n}] T q + [fi , } M T _ {p-n} 1 q }, } ^} + [n , ] q-+{(t k {p-n} 1 }, M' k + BX - q 0 U ), +Q p q g +BX }. q Note that Hamilton's equations, derived above, can easily be obtained from Eq.(4.13). P A R T II : A P P L I C A T I O N TO O R B I T I N G S Y S T E M S 96 5. A T T I T U D E 5.1 DYNAMICS OF GRAVITY FLEXIBLE SPACECRAFT GRADIENT Spacecraft M o d e l Consider a satellite consisting of a central rigid body with an arbitrary number of plate and/or beam-type flexible, deployable appendages fixed to it in desired directions. The satellite is free to negotiate a specified trajectory. Let the position vector Rc and the true anomaly 0 define the location of the instantaneous center of mass O of the spacecraft with respect to the inertial coordinate sytem X, Y, Z having its origin at S, the center of the Earth (Figure 5-1). A n orthogonal orbiting reference frame X0,Y0,Z0 with its origin at C is so oriented that YQ and Z0 are along the local vertical and horizontal, respectively, while X0 is aligned with the orbit normal. The body coordinates xo,yo,zo with origin at C coincide with the orbital coordinates in absence of any librations. The orientation of the body axes £0,2/0,20 a ^ a n Y instant relative to the orbital coordinate frame X0,Y0,Z0 can be described by a set of Euler rotations. The angular velocity in the body coordinate system can be written in terms of the librational angles and velocities. Here we choose the rigid body as the main body By, and the system's instantaneous center of mass as the reference point. Fixed to the main body, at its center of mass, are the body coordinates £1 > 2/1 > 21 • At any instant, the orientation of a given appendage with respect to its nominal undeflected configuration is denoted by a body coordinate system xt-, 2/1,2;; i > 1. For a plate type appendage, yt- is taken normal to the nominal plate plane, Z{ along the direction of deployment, and £ ; normal to Z{ in the plane of the plate. For a beam type appendage, Zi is along the nominal beam direction, while x,- and y,- complete the orthogonal set (Figure 5-2). 97 Figure 5-1 Geometry of orbiting spacecraft with beam and plate type appendages. Figure 5-2 Reference coordinate systems for a spacecraft with beam and plate type appendages. 100 5.2 Kinematics As described earlier, it is convenient to represent the vibrational displacements in terms of a set of admissible functions which are somewhat arbitrary as long as they satisfy at least the geometric boundary conditions. The exact modes of vibration of a rectangular plate clamped at one side and free at the remaining three sides are not known. Hence, the modes of a fixed-free and free-free beams are used to describe the plate elastic deformations. The lateral displacement v for a given plate and the transverse oscillations Ub, v&, for a p given beam, in orthogonal directions x and y, respectively, are assumed to be of the form : Vp {x,Z,t) =m X X *mn(*) <Pm{x) tjj {z) ; n n u {z,t) = Y n b S {t) ijj {z) ; Xn n v {z,t) = ] T S {t) ip (z) ; n b yn n where : f>mn — 6 ,8 Xn yn unknown generalized coordinates associated with the plate vibration ; = unknown generalized coordinates associated with the beam vibration in the x and y directions, respectively ; tp n = characteristic modes of a fixed-free beam (Figure 5-3a) ; tp m = characteristic modes of a free-free beam (Fiugure 5-3b) ; <Pm = shape functions for a cantilever plate ; n, m = number of modes considered in the analysis . The characteristic functions ip and tpm, their properties, and associated integrals n of interest are presented in Appendix D. With the modes and their properties identified, the first step is to present elastic deformations in the matrix format : normalized undeformed axial coordinate Figure 5-3 Vibration modes for fixed-free and free-free beams (a) the first four modes of a fixed-free beam ; (b) the first two modes of a free-free beam . 102 • For a beam-type appendage, the matrix of admissible functions has the form <f>u _o 0 0 <f> 0 0 0 0 {«i}ii 2i = { o*} ; = { } h 5 = re ^ ^ • For a plate-type appendage, the matrix of admissible functions has the form 0 0 0 0 0 <?i 0 0 0 {/C,} 2 2 = { <t>i } ; fa = { <Pm }i • • The generalized coordinates associated with the beam and plate-type appendages, respectively, are * = <*•><• • The unit vector in the direction of deployment, relative to the body coordinates, for both the beam and plate appendages is given by For the system configuration under consideration, the relevant kinematical parameters with reference to inertial coordinates can readily be identified : a) The matrix of geometry S ; 6 J k is such that : i *L°k f ? Xj-i,yj—i,Zj-i Xj,yj,Zj, j > |5 ^ T, \' . for all A; if j > 1 . \ S!- = 0, J b) The directions of the coordinates dinates *' 1 1, are related to the body coor- through the transformation {x } 3 = Cj - 1 {xy_i}. As the coordinates are time invariant, it follows that uJy =0, c) The location of the body frames Xj-i, Xj,yj,Zj, j ? 0 . j > 1, in the coordinate system yy_i, Z j - i is given by the constant vectors ay, therefore 103 5.3 Kinetics The integrals required to account for flexibility of the appendages (Eq.A.4) are shown in Appendix E while the basic kinetic quantities arising from the flexibility, as defined by Eq.(A.2), are listed in Appendix F. The dynamic quantities, as defined by Eq.(3-4), for beam and plate-type appendages are listed below. Note that I, d and m represent the instantaneous length of a beam or a plate in the deployment direction, the width of a plate, and the instantaneous mass of the beam or the plate, respectively. Terms on the r.h.s. of the expressions are the integrals due to the flexibility as presented in Appendix E . Beam : —T — flit —T — /22 6 1/2 i { c i H 6 = [/u / v* y - 2 2 2 s 0] ; 2 {-tjf 2 + d t 1 \6 [Ell-E\l T L ee f [ FU = 1+ 1 6 T = + / 2 2 } T s Mi-ti-fnVs I w 2 + Dl\ -Dir]6 T J -^22 [ C{\ + CH + F i Ell + EM - D\\ l Y +V - D\\ ]; 2 2 2 2 - 2 Nil - 2 Nil ] * ; 104 Wss = [ {Fil + iF2 22 2 ] ; [Jx] = | m 1 0 0 0 1 0 0 0 0 e 0 m r 0 -t{gh} 6 o -nghYI T Sym 0 6 [F%}6 -6 [F**]6 T 0 T m 0 6 IFR)6 T Sym. 6 [FR + Plate —T b= 1*22 6 I 2 h = P H s = [0 ( 22 ^22} ^ c J / T 0] ; 2 2 I {-v% + d 2 + / 22 2 2 } 6 T 0 d {v£ 2 d\ } 6 T 2 0 -d c " w ee = 1+ y w eS = d { v\ T 22 {^ - <4} «+ ^ T 2 - d\ } + j 8 T 2 2 T [ C 2 2 2 2 [ E ^ - D\l } ; 2 + V% - 2 N 22 ] i 105 W = [ iF% ] ; ss 4£ 2 0 -3£d m 12 0 4(£ + d 2 ) 0 0 -d 0 {g\ } S T 2 0 [/2] = m o 6 [F%]6 0 ['s] = m -3£d 0 4d 2 2 0 Sym. 0 0 6 {F%)6 T In the present case there are no holonomic constraints inside or among any of the system domains and the kinetic energy, Eq.(3.29), has the form T =- 2 The x~ T Tvj 1 M Tv X~ V mass matrix of the system, Eq.(3.27), reduces to -W cT o T 3 i i f S 3k^3 ie \ W 3 0 \-mjkK kK' C T - V Cj 0r,C0 0 U 0 CgH°C r<o pfcr 0 T C0 Sym. T •}W c{ b k k 0 -)WcJ 0 0 fc V D' J k 1 < a 2 0 u • h 0 7 wo ( 6 0 J Sym. £ JkT 0 a0 < w jpjk I fcT orkCkHs 6 -m}kH 3C 3kH kT\] The kinetic energy can now be written as , - SjkmjW 1^ 6ji.Sk! ' m? 3 -m]k he C k H« 0 • A; & J T 106 where : I = CZH°C D jk ; 0 = Sjihi (Sjkmj-w E = 8 m C^[-C P^ k kl F jk 9 k 3 - m V C\Ji^ T ee jk + f C Hl ] T 0 k h jkT ; ; T k k = SjAtlSjkmjW's - m H CiH \ ; j jk jT e k T 6 s 8nmAt CJ'-K Cjor-j]Co; = T T = hAi^km^J - m h c{H ] J T jk e kT . In attitude dynamics of spacecraft the motion is described relative to an orbiting frame, hence : tUo = uJ + fin ; oo : Ret + £lR u ; c where : u = W0~; f = (sin0 0 cos0) ; u T = (cos0 0 - sinfl) . Here, 6 is the orbital angle and VI the orbital angular velocity. The unit vector n is along the orbit normal. The spacecraft angular velocity with respect to the orbital frame is denoted as U. It is expressed in terms of the modified Eulerian rotations, where 6 is a vector formed from the Eulerian angular velocities and W is as defined in Appendix G using the /?, 7, a rotation sequence (Figure 5-4). 107 7 roll Figure 5-4 Modified Eulerian rotations showing yaw(/3), roll(7), and pitch(a) librations . 108 The system coordinates can be written as, rn? b "1 0 V a R t + QR ti 0 c U0 c 57 _u' . 6 - + fin I 0 0 0" 0 1 0 0 0 0 0 0 0 0 0 0 0 n W 0 0 0 0 .0 J & 0 t R ti c 0 1. Introducing the transformation matrix T as defined above, the kinetic energy takes the v form, T rn? I \ \ 2 -W l R T l 0 k -W 1 0 T D' 3 k R t + UR u c 0 h 9 jT u c 0 + VLn We I T J r JS^C&S* V 2 D' 3 1 R 1 fi > -}b Cp 3T R 5.4 c r fin 0 u •k 6 f •k\ m 0 l ^ T k W o J2 tf + n / n T c n IW nE W IW WE T T T k T pjk k R o < 2 Sym. { gJT 3 T c Wt+ 3T g'n c k c 0 R t + QR u < -R )b C]u c UJ 0 pjk _ J7 c < 0 j k T E Sym. m o - 3 k fi > t •k [6 J Potential Energy From Eq.(3.31), the gravitational potential energy considering the system reference point as its center of mass can be written as 109 The higher order terms in Eq.(3.31) are omitted as their effect on dynamics of the system under consideration is expected to be negligible. On the r.h.s. of the above equation, the first term represents potential energy due to the satellite treated as a point at its mass center. The other two terms are due to the finite size of the system. The strain energy stored in the appendages during their vibration is given by v = J2 f+Il * h> v v j=i k=i where n and rib are the number of plates and beams, respectively. The strain energy p expressions associated with a plate and a beam are : Here G p and G b are the flexural rigidities of the plate and beam, respectively, and v is Poisson's ratio. Substituting the expressions for the plate and beam displacements in terms of the shape functions and the generalized coordinates, the above equations take the form: E E E E m where J^s , J%3 , Jmr , J m r n r +2d - » ) J i J i j ] ; 8 are constants given by : T2 rp f d ^ n d l j j a 110 dz . The values of the above integrals are given in Appendix D. Note that the system strain energy can be written as where Here mn and s define the number of shape functions associated with plate number j and beam number k, respectively. The diagonal and off-diagonal terms in the V matrix can be identified from the strain energy expression. Note that for a beam the strain energies for different modes are decoupled, while in case of a plate they are coupled for theflexiblefreefree modes (6j j — 3,4, k A; = 1,2, •••). The strain energy expression developed here implies small elastic deflection and rotation with the shear deformations ignored. Thus they correspond to the first term on the r.h.s. of Eq.(3.39). 5.5 Equations of Motion The equation of motion can be written as where the system velocities are Ill The system mass matrix is w ci\b T k o pit M* = •)V Cjt T 0 -R )V Cju T c 0 gn gW u Rt u o jT iT T r R\U + n In T n IW T W IW T Sym. T jk T 0 n' E T k WE T k F> k As explained in Chapter 4, the array of system velocities x contains both the system generalized and specified velocities, q and s, respectively. Therefore, the system generalized coordinates q, mass M , and coupling momenta n depend on the case under consideration. Here the stiffness matrix K is full of zeros except for the submatrix, associated with the appendage's elastic deflection, which is equal to V as defined in Section 5.4. The e differentials appearing in the above equation were treated in detail in Chapter 4 and the associated Appendix B. The set of n second order governing equations, defined above, are transformed into a set of 2n first order equations, as required by the numerical routine solver : with 5.6 I £]{!}+[£ ^{fM/l-; Static Equilibrium and Stability of the Model For a triaxial gravity gradient rigid satellite in a circular orbit, Likins and Roberson [62] have shown that there are 24 equilibrium orientations. An equilibrium state is obtained when any one of the three principle axes, each with two senses, is aligned with the local vertical. For each of these six possibilities, either of the two transverse axes, each with two senses, may point along the orbit normal thus giving the 24 equilibrium positions. The torque on the satellite, due to gravity and centrifugal forces, is identically zero only 112 for these 24 equilibria. It is shown that for the nondeploying flexible gravity gradient spacecraft under study, these 24 equilibria for the attitude motion are still valid although the flexible appendages are, in general, deflected. Furthermore, an analysis of the deformed configurations of the flexible satellite in equilibrium and their stability are presented. Stability of the equilibrium orientations is studied here using the system Hamiltonian H as the Liapunov function. The attention is directed towards a class of problems where the Lagrangian (and hence the Hamiltonian) is free of explicit time dependence. From Eq.(4.9), H = t - t + U = t + P , 2 where to and t 2 0 p 2 are functions independent of and quadratic in the system generalized velocities, respectively. U is the system potential energy and P = U — to is called the p dynamic potential. As t 2 p is positive definite in the system generalized velocities, and the free constant in P can be so selected as to make the Hamiltonian H zero at the origin, H is also positive definite when P is positive definite. By requiring the dynamic potential P to be a local minimum at the origin, sufficient conditions for the system stability are obtained. The procedure, due to Poincare, is applicable to conservative, holonomic systems with the Lagrangian free of time dependence, and was first applied by Auelmann [63] in the satellite stability study. Note, the method represents an extension of the case when the Lagrangian is natural (ty = to = 0) and the dynamic potential P is merely the system potential energy. A n earlier theorem attributed to Lagrange states : if in the static equilibrium state the potential energy of a conservative, holonomic, scleronomic, mechanical system is a local minimum, the motion is Liapunov stable. Following Auelmann, Pringle [64,65] also applied Poincare's approach to study extensively the stability of spinning satellites. He extended the method to encompass a wider range of systems by the following theorem : a holonomically constrained mechanical system with: • dL/dt = 0; • only dissipative nonconservative forces ; 113 • damping both solely generalized and complete in the noncyclic generalized coordinates ; is: (a) asymptotically stable in these coordinates if the dynamic potential energy P is positive definite; (b) unstable in these coordinates if the dynamic potential energy P can take on negative values arbitrarily close to the origin. For the system under consideration, the specified coordinates are defined by {s} (ra I J R 3 O) and the generalized coordinates by c {q} T = = (U S ) , hence the kinetic T 3 energy of the system can be written as T = t + h + t 2 , 0 (5.1) where - U= -q 1 U = - M q 2 \ 6 \ T 3 J T T nE - T }W ci \b k V = 3k k Uf 2 \ h 0 D -)W 3 j k cj t 0 + nE Js z,-ijl? k kT -R I6 J )V Cf u l T c 3 3 k ' m ' k 3 £ k < Rc u Sym. Note <2> * i k e n IW 0 WE F T \v T t W IW Sym. a n R t T c R 2 C u > Rc U+ n I n. T d to are components of the kinetic energy and correspond to quadratic, linear and independent terms in the system generalized velocities (Eq. 4.6), respectively. Let the system be so oriented that its principal axes coincide with the orbital coordinates. The dynamic potential is defined as P = -t 0 + U , p 114 = -t -±gtrI+^gf U + ± 6 0 V6 T , 9 = ^ . (5-2) Here the gravitational and strain energies from Section 5.4 are used to define the potential energy, U , of the system. The first term in the gravitational potential energy, —fj, /R , is p e c omitted as it is independent of the generalized coordinates and hence does not contribute to this analysis. The terms in the kinetic energy, which are independent of the generalized velocities are denoted as to, which is defined by Eq.(5.1c). An examination of Eq.(5.1c) suggests that the terms which render the dynamic potential time-invariant are those associated with the specified coordinate represented by the orbital angular velocity $1, since it is constant in the case of a circular orbit. From Eqs. (5.2) and (5.1c), the dynamic potential can be written as : P = -— n 2 I n - -g trl + -g f T 2 2 y y I 1+ - 6~ V 6 , T 2 3 = X o o I(j,j) j = - i n nj 2 3 = - i ( fi 2 = hUJ) Here ray, tj, and my, + - 6 VS; (5.3a) T \g ( n) + t) + m) ) + \g l) , + g ) n) - ±g m) + g l) ; + [ {hyjj 6) + {6 T [Js]yy 6 ] . (5.36) (5.3c) j = 1,2,3 are the components of unit vectors n, I, fn, which are normal to the orbital plane, along the local vertical, and in the direction of motion, respectively. Thus, for example, n\, n , n represent components of n along the system 2 3 coordinates x , yo, z , respectively. 0 0 The attitude equilibria are obtained by equating to zero the partial differentiation of the dynamic potential with respect to the attitude generalized coordinates. Using Eq.(5.3): 3 P,e r Oj,er = Yl °J>°r IU>J) = . 0 = - ( n + o ) ny nj,0 2 r r = 1,2,3 ; -g rrij rrij,er +2g tj tj,er , (5.4a) j,r = 1,2,3 . (5.46) 115 Here 0 , r = 1,2,3, correspond to a, f3, 7 describing the pitch, yaw and roll librations, r respectively. Substituting forray,£y, my,ray,^ , ij,e , my,^ from Appendix G into Eq.(5.46) r r r gives Oj,e = 0 at 0 — 0, hence Eq.(5.4o) is satisfied and the reference orientation (a = r r /3 = 7 = 0) is an equilibrium position. Thus, the 24 equilibria for the rigid satellite are still applicable to the flexible system under study. Similarily the equilibrium positions of theflexibleappendages are obtained by equating to zero the partial differentiation of dynamic potential with respect to the deflection generalized coordinates. Using Eq(5.3), 3 P °ji >T = E j=i Uth, + 6} + V6, 2 (5.5a) 3 °iibh* = 0 • = [ E °i [H-y + 1 +E v 2 y ( 5- 56) y y=i The equilibrium configurations for the deflected appendages are given from the solution of Eq.(5.56) as * = -I 2 E °i \h\a + v r 3=1 1 oAhhi } > {J2 3 = |[-ff[(n2 + l)[/3]ii-2[7 3 ]22 + [ / 3 ] 3 3 ] + ^ ] " 1 { (n 2 + 1){/ 2 }H - 2{/2}22 + {/2}33 } , (5.6) 2 where 0 = Cl /g. Note, the above uses the values of oy, j = 1,2,3, as given by Eq.(5.36) 2 at the attitude equilibrium points, which are : 0 1 02 = - ^ ( f i 2 + o) ; = 9; 03 = -\g • The sufficient conditions for stability are obtained by requiring the dynamic potential to be a minimum at the equilibrium points, hence its Hessian is required. Using 116 Eqs.(5.4) and (5.5) : 3 Pi6 9s = ^2 °j,6 e r r (5.7a) 5 I{j\J) a j=l Oj,e e r = - ( ^ + y)i j,e 2 a + 2o(£y,^ nj,e n a tj,e s r +j n r +£j Lj,e 6 ) r a nj,e e ) r a - g{mj,e a mj,e +j m r "*y»* *a) r (5.76) ; 3 (5.7c) y=i 3 (5.7a") y=i where j,r,s = 1,2,3. Substituting from Appendix G in above gives, at the equilibrium attitude points ( 0 = 0, r = 1,2,3 ) : r Oj,e P,s e r = 0; a P,e e r — 0; a r ^ s ; = 0, a hence the Hessian of the dynamic potential has the form 0 0 (o) • 0 p,pp 0 (0) 0 0 {0} {0} P where : (0) row of 0 ; {0} column of 0 ; P,gj matrix. iaa p {0} (o) P>66 - The sufficient conditions for stability are obtained by requiring that the determinent of the Hessian at the equilibrium points to be positive. Using the relations of Appendix G, the conditions for stability are : dP 2 do? = 3a (J -I ) 3 2 > 0 ; (5.8a) 117 | ^ = - n 2 (I, - h) > 0 ; (5.86) dP 2 — = 3o (A - I ) - fl 2 (J - h) = (3</ + n 2 )(/! - h) > 0 ; 2 2 (5.8c) 3 2 P >55= I E °i[h]a + V\ > 0. (5.8d) Note that, since the dynamic potential is quadratic in the deflection coordinates 6, the equations of motion for the vibration degrees of freedom are linear and the condition of static stability reduces to the requirement that the stiffness matrix be positive definite. Note, Eq.(5.8d) can be obtained from Eq.(5.6), which defines unique equilibria for the elastic deflection. Also the attitude instability, which is governed by the nonlinear equation, implies the system acquiring another equilibrium orientation, while instability for a flexibe appendage means unbounded deflection, as shown by Eq.(5.6). From Eqs.(5.8a), (5.86) and (5.8c), respectively, for stability : h > h ; h > h; h > h • (5.9) Hence, the sufficient condition for attitude stability is (5.10) h > h > h • The condition requires the satellite principal axis associated with the maximum moment of inertia to be aligned with the orbit normal, and the principle axes associated with the minimum moment of inertia to be aligned with the local vertical. The orientation is referred to as the Lagrangian configuration. Of course, i i , J , Is are the principal moments of 2 inertia for a satellite in its deflected equilibrium configuration. The condition defining stability of a flexible appendage, Eq.(5.8d), involves information on the appendage rigidity, orientation, location, and the orbital parameters. In general, this condition can be satisfied if the appendage rigidity is large enough. However, limitations on satellite weight, may impose restrictions on the alternatives available to a 118 design engineer. Fortunately, in practice, the rigidity is usually adequate to satisfy this sufficiency condition for static stability. To account for system dynamics would require inclusion of additional forces and one may have to resort to a numerical simulation. The principal moments of inertia about the system coordinates for an appendage, with its body coordinates arbitrarily oriented and located at the system center of mass, can be obtained from Eq.(3.21) : —T T — n 1 n . EE i n y=i fc=i kU'fy ' Hk 3 3 E 3 + E E ( ~ i' ) 1 y n J") 3= 1 s k n J n *)' k 3 = 1 k=l c 2 f t c 21{ 2 c <Xi c 2Pi s 2ii + a a,2 \ T , a2ft - 2 soti ca, aft eft s 7i + c a,- a f t + 2 aa< cat aft eft - 2 ( ca.i c /3i sn C7< - aa,- aft eft cu ) s 2cti 2 c 2Pi 2 2 s^i s^j 2 2 ( s a ; c 2 f t 37,- c~a + ca,- aft eft C7< ) - J?(2,2) = / r * ( M ) ) J< (2,2) 2i(3,3) 4(1,2) 4(1,3) I 4(2,3) J (5.11a) / £, 3 3 =EE ^ ^ > 3 = 1 k=l 3 =E 3 £ y=i y + 3 EE ( -^ )4* c 32 *) • £ y=i fc=i * 7f a i c 7< a * c 7i 1 2 2 1 2 4(2,2) 1,(3,3) 4(1,2) 4(1,3) I 4(2,3) J 2 2 ca^ s^j C 7 i —2 saj S7i 07; —2 s a » ca,- c 7i . 2 (5.116) 7° (3,3) = m I m T 3 3 E E m 3™ k 40',*) ' j = l fc=l 3 3 E ? m 3 +E E ( ~ 1 ) J m mfc fc ) > 119 T a" 2fa s n + 2 sa, s 2Pi s27,- — 2 sa,- ca,- S/3,- eft 37,- + c2c*i 8 2cti c 2Pi c2a,- c 2f)i + s 2ai 2 cai sft r A(I,I) 1 ^ (2,2) > i /<(3,3) 4(1,2) 4(1,3) I 7*2,3) J eft 37,- -2 ( ca,- s 2 ft sn 07^ + sai a/9,- eft 07,- ) 2 ( -ecu,- sft eft 07^ + sai s 2 ft m m ) (5.11c) Euler angels a,-, ft, 7,- in Eq.(5.11) define the orientation of frame L,- relative to the system reference frame LQ, which now coincides with the orbital frame. From Eq.(5.8d), which defines the appendage stability, and Eq.(5.6), it is apparent that only the terms in the moments of inertia, which are functions of deflection, i.e. Ii,Iz, contribute to this analysis. Introducing the values of g , F$ } mn n (m,n,j,k = 1,2,3) from Appendix E in the expressions for i^jLs (Section 5.3), and substituting them into Eq.(5.11), lead to linear and nonlinear terms in deflection, in the moment of inertia matrices, as follows : Beam : ( ca, sPi cfii C7i + sa,- c pi 57,- C7,- ) / A 2 {/ }n = - 4 £ m ( sPi cpi 57,- ( s a,- - c a , ) + sa,- ca, ( s ft - c ft s 7,- ) ) / A / ' 2 2 {J } 2 2 2 IT\ -Urn j 33 H / I m ii > 2 c --Af 4 l m j m A ( ^ ^ 2 2 CCLi ( c S P< a . _ a 5 C 2 a . ) l + i + S C i i ^ rn ca S2 ^» _ ( S7 » e 2 f f i C7 » )/ « \ • A _ p. s2 5 2 . 7 ) I, ) / A 2 (sa sft cpi C7,- - ca< c p 57,- cn)U ( sa spi cpi C7i - cai c Pi 57,- 07,- ) U ((eft C 7 , ) + (ca,- eft + sa, eft S7,) )t/ 2 2 { 2 c27,- U -cai <s7i cyt *7 { 2 2 —cai sn C7,- U (s 7,- + s aj C 7 J ) Z7 2 2 ( c Pi + s pi 2 [73)33 - 2 2 ( s ft + c pi s27,- ) U { [73)22 = rn 2 j , 2 ( ~ 2 [/3]u = 2 2 s27,- ) U 2 ( sai spi cpi C7,- + ca,- s ft 57,- cn ) J7 2 (sa sft cPi C7,- + ca,- s ft 57,- cn)U 2 { ((5ft C 7 , ) + (ca,- eft - sa,- sft 2 sn) )U 2 (5.12) 120 where A are the roots of the characteristic equation for a fixed-free beam, and U is the n ns x ns unit matrix. Plate : {h}n - - 2 m | ( sai spi cPi c^i - cai c 2Pi sn cy,- ) d (6 ( sfit cpi an ( s 2a{ - c 2a{ s^i m d(6 - ml ) + sa{ cat ( s pi - c 2pi s 21{ 2 {J2/22 = 4m { -cai {/2>33 = ~4m | ( sai spi cPi cn + cai s Pi sn cy- ) d (<5 i - [J ]u = m ( (cPi c ) 1{ 3 3 22 - -^=6m2)pn t = m ( 8*ii + s 2ai c [/alas = m ( {afit c ) 2 1{ 2 7 i t m l m ) [U] p{ s1{) 2 C 2 }5 2 } ; m 2 - c 2a{ ) - sa{ ca{ ( c 2Pi - s 2p{ s 2^ + (cai sft + sa{ 2 ) ) &5 i/A + sa{ ca,- c -y- £ o / A 2 ( spi cpi S7,- ( s 2a{ l7 ] ml -j=6m2)Pn+ 2 }; -^=6m2)Pn- ) ) £5 i/A m ; ) [U] ; + (ca,- eft + sa sp{ s-y ) ) [U] . (5.13) 2 t t Here, the subscripts m and ra correspond to the free-free and fixed-free beam modes, respectively. Pn = an/\n, where on are constants dependent on A (Appendix D). n In summary, the equilibrium position of a cantilever plate in orbit is defined by Eqs.(5.6) and (5.13). Similarity, the equilibrium of afixed-freebeam in orbit is defined by Eqs.(5.6) and (5.12). Their stability is expressed by Eq.(5.8d) and Eqs.(5.13), and (5.12), respectively. Of course, in general, the appendage may have any arbitrary spatial orientation. These relations governing the appendage equilibrium and stability are functions of 121 the attitude parameters (a, (3,7). The situation can be tackled quite readily by locating the appendage in a specific plane and allowing it to take various positions in that plane. This is illustrated in Figure 5-5 through three examples. Beam in the orbital plane Setting A> = 7t = 0 (Eq. 5.12, Figure 5-5a) gives {/.}.. = { ° } ; [J ]u 3 - m 0 0 0 u {/.}» = « m { [h]22 s a b c a b K {hhs = -Urn { ^ } ; 0 U 0 = rn ° sa U 0 U 2 b ^ } ; 0 ca 2 h U Substituting from above in Eq.(5.6), using the expression for the matrix of strain energy V from Section 5.4, and noting that Q = fi /o = 1 for a circular orbit, 2 e = - 6 £ f]*2 0 -U 0 6 = — 6 £ gm -gm u *2 0 {-2s a 2 + ca 2 b 0 -3c a U + f] + 1) U " ( 2 { sa l /X £3 ^ ] + 2 h ° 2 b b n } ' \ \ sa c a / A J 2 6 b 6 Here 0* = £ m n / G . Note that U is anragx ns unit matrix, and A is a 2ns X 2ns diagonal 3 2 2 b matrix with elements defined by A(j, j) = A(j + ns, j + ns) = Ay j = 1,2,..., ng. From the above expression for deflection, the following observations for a beam type appendage in a circular orbit can be made : • The deflection associated with the n-th mode of the beam in the y-direction is given by -6 £ f T sa ca 2 8. y n b X ( -3 Q*2 c a 2 b + A£ ) 2 b 6 £ fi* sa ca 2 b A2 ( -3f2* c a 2 2 b b + 1) ' where 0*, a dimensionless quantity associated with the n-th mode, is defined as • At a = 0, 90°, 180° and 270°, b 6 yn = 0. 122 Figure 5-5 Geometry of a beam in orbit for static equilibrium and stability study: (a)beam in the orbital plane; (b)beam in the plane normal to the orbit; (c)beam in the plane tangential to the orbit. 123 • The maximum deflection in the y-direction is for the beam inclination (a&) defined by tan a 2 b = 1 - 3H; . 2 When f2* << 1, the maximum deflection occurs approximately at a « ±45°, ±135°, b and is given by [<$ J ax = T3 I fi^/A m n /X = ^3 I 4 2 y m 2 6 G . b n • The stiffness in the x-direction, which is associated with the n-th beam mode, is equal to m fi + A^ G / £ m . The second term represents the structural stiffness 2 6 3 of the beam, while the first one, of the geometrical type, is due to centrifugal and gravity forces. The geometric stiffness element is a positive constant for all modes in the x-direction . • The stiffness in the y-direction, associated with the n-th beam mode, is equal to —3m c a f2 + \„ G /£. m. 2 b 2 b The second term represents the structural stiffness 3 of the beam, while the first one is of the geometrical type. The geometric stiffness element in the y-direction varies with the beam inclination (a ) in the orbital plane. b It has no contribution to the beam stiffness at cc = ± 9 0 ° , and a maximum negative b contribution at a = 0, 180° of - 3 m U + A* G /£ m. 2 b 3 b St rn 3 For static stability, G > b fi /A£. 2 Beam in the plane normal to the orbit (out-of-plane) Setting fi = 90° and = 0 (Eq. 5.12, Figure 5-5b) gives : b {/,}„ = -«m [I:3 11 m { ^ U 0 ° c a b / x l 0 ca 2 b U } ; { h h = 4£m { ^ 2 [J ] 2 3 2 -m U 0 } ; ^ 0 sa 2 b U {/ } 3 = 2 3 { ° } ; 0 0" 124 S = -2lgm {T? -l)U -gm (-2s 2ab 0 0 = -8£ n* 2 + [l + n )c 2ab + + l)U £3 1(3 + n 2)sab J 2 cab /X n j ' 0 (s 2ab 0 - 3c 2ab) A4 4 + U ] {4sabcab/\ 2n} ' Based on the expression for deflection of a beam in orbit, the following remarks can be made: • The deflection associated with the n-th mode of the beam in the y-direction is given by 8 I U* sab cab 8 I Q* sab cab 2 2 X { n* (s a y n 2 2 - 3c a ) + X* ) 2 2 b b • For a = 0, 90°, 180° and 270°, 6 b • When O* yn X ( n* (s*a 2 2 n b - ^a ) b +1 ) = 0. 1, the maximum deflection occurs approximately at ab — << ± 4 5 ° , ±135°, and is given by [<5 Jmax = T4 £ m f i / n ° 4 2 A h y • The stiffness in the x-direction, which is associated with the n-th beam mode, is equal to A^ G /£ m. b It is solely due to the structural stiffness and the orbital forces 3 have no contribution. • The stiffness in the y-direction, associated with the n-th beam mode, is equal to m(s Q!6 — 3c ct6) f l + A G It?m. The geometric stiffness element has no contribu2 2 2 4 b tion at ab =• ± 6 0 ° , ±150°, and a maximum positive contribution at ab = 90°, 270°, of m Cl + X* G /C m. The maximum negative contribution in the y-direction of —3m Cl + X* G /£ m occurs at a& = 0, 180°. The static stability condition is 2 b 2 G b b > 3£ mH /A 3 2 3 3 125 Beam in the plane tangential to the orbit Setting 76 = ct = 0 (Eq. 5.12, Figure 5-5c) : b {/.}» = - « m { * } ; { s ft U 0 0 u h h 2 = {°} ; U 2 [is] n = rn 2tgm -gm -2£n *2 m 0 2 c ft 0 0 -3 2 U •*r { / } 2 0 0 «*A c (3 U 0 0 u 2 [is]ss = 1"! —2, (2 + f2 ) J7 aft cft/A 0 = 4£m { * 3 3 0 2 C *2 : = ((l + n )s ft + f t - 2 ) U 2 6 = [is] 22 { b m -fi aft cft/A 0 2 2 From the expression of deflection, the following can be concluded : • The deflection in the x-direction, associated with the rc-th mode of the beam, is given by -2 £ ft* sft eft -2 £ n* 2 A2(fi*2 X n • At ft = 0, 90°, 180°&; 270°, • When f2* << c2/?6 + 6 Xn A 4) A n 2 sft, cfo 2(f7*2 2^ c + 1 ) • = 0. 1, the maximum deflection occurs approximately at ft = ± 4 5 ° , ±135°, and is given by (<5 )max = Xn m n /X 2 6 n G b. • The stiffness in the y-direction, associated with the n-th beam mode, is equal to —3m f2 + Xn G / £ m . The geometric stiffness element is constant and negative for 2 6 3 all modes in the y-direction. • The stiffness in the x-direction, associated with the n-th beam mode, is equal to m c /3b ft + 2 G hjH?m. 2 The geometric stiffness element in the x-direction varies with the beam inclination (ft) in the orbital plane. It has no contribution to the beam stiffness at ft = ± 9 0 ° , and a maximum positive contribution at ft = 0, 180° of m H + A G / £ m . 2 4 6 3 126 From the expressions of deflection, the following can be concluded for beams in various orientations : • Deflections in the x and y directions are statically decoupled. • Deflections associated with various beam modes are statically decoupled. • The static deflection in the z-direction is zero for all beam positions in and perpendicular to the orbital plane. The static deflection in the j/-direction is zero for all beam positions in the orbit tangent plane. • The ratio of the contribution of various beam modes to its static deflection under the effect of orbital forces (centrifugal and gravity) when fi* << 1 is given by 2 1 (^1)8 : : ... (-^) : 6 « 1 : 4.07 x 1(T 3 : 1.85 x 10~ The limiting stiffness for a beam in the orbital plane is G B 4 : 2.46 x 10" : ••• 5 > 3 £ m fi /A = G , 3 2 4 0 where GQ is a minimum base beam stiffness. Theoretically, with a stiffness slightly higher than G the beam is statically stable but may suffer a large undesirable static deflection, 0 defining its equilibria, not to mention the effect of the dynamic forces. Therefore, in practice the beam stiffness is much higher than such a calculated value. Let the beam stiffness G used in design be a multiple k\ of the base stiffness GQ, where ky > 1. Thus, B let G = 3 ky £ m fi /A . The maximum values for the beam deflection in equilibrium B 3 2 4 configurations are given by : • for inplane beam : 6 = l/(ki A ) ; 2 Vn • for out-of-plane beam : 6 = 4 £/(3 ky A ) ; 2 yn • for beam in the horizontal plane : 6 The maximum tip deflection is, — £./(3 fci A ) . 2 Xn 127 Let the allowable deflection be a fraction k2 of the final length, where k2 < 1, i.e., ^ - 8 V" 1 From Table D - l : -2 « 0.2844 ; A A x » 0.0454 ; 2 - y « 0.0162 ; A A 3 « 0.0083 4 Therefore, k2 8 x 0.3543 « — — 3 Ki « 0.95 1 -j— « — , «i i.e., fcx , ki « 1 — ^2 This approximate rule can provide a rational start to the design of beam-type appendages. For example, one may require the maximum static deflection to be 1% of the beam length. In that case, the beam stiffness should be approximately 100-times the base stiffness GQ. This can be used as an approximate starting value to be followed by the numerical verification through a dynamic analysis. Note that by considering the displacement field to be linear and keeping only the quadratic terms in the strain energy, as is the case in this application, the calculated geometic stiffness is different from the actual one. The extension of the above analysis to include a nonlinear displacement field is straightforward. Similar approach can be used for plate and other types of appendages. Cases when the Liapunov function is time-varying and bounded may also be studied with the orbit eccentricity and specified deployment history included. These extensions represent useful topics for future work. 6. R E S U L T S A N D D I S C U S S I O N 6.1 Remarks on Computational Approach 6.1.1 Computer program The governing equations are coded using Fortran-IV programming language in a modular fashion, which helps isolate effects of nonlinearity, elasticity, deployment, etc. The computer code simulates the dynamics of a cluster of objects in an elliptic orbit. The cluster consists of an arbitrary number of flexible deployable appendages attached to a central rigid body. The appendages can be of beam or plate type, and are located in any desired orientation. The deployment time history can be sinusoidal or constant with a specified duration. The cluster orientation relative to the orbital frame is described by a set of coordinates which can be treated as generalized, specified, or a combination of both. A specified orientation may arise through a variety of controlled maneuvers including spin, slew, etc. Figures 6-1 and 6-2 show schematically various computational paths in the program, its inputs and outputs. The program input data may be classified as: (a) geometrical and dynamical parameters of the system such as moments of inertia, appendage type, mass, location, orientation, dimension, rigidity, number of modes, etc. (b) time histories of the specified coordinates such as deployment, slew, orbit, etc. (c) initial conditions. The input control parameters serve several objectives which may be divided as: (a) control parameters defining choice of various alternatives such as the method of integration, matrix inversion, optional outputs from the program, etc. 128 129 SPACECRAFT T orbit "I" I circular elliptic components and flexibility — \ rigid I — single b o d y flexible T cluster — plate beam 1 cluster _zr T deployment velocity 1 ] -f-ve 0 1 -ve ; deployment history r constant sinusoidal X J librational coordinates J specified shuttle control ) partially specified r— spinning Figure 6-1 T_ sinusoidal generalized r—^ nonlinear A flow chart showing an approach to system simulation. linear generalized coordinates specified coordinates graphs of stress diagrams system parameters initial computer code for shape conditions system dynamics functions control parameters brief calculations detailed case documentation Figure 6-2 A schematic showing input to and output of the computer program. 131 (b) control parameters defining the system in terms of its configuration, flexibility, deployment, etc. 6.1.2 Computational considerations The equations of motion were analyzed using the main frame computer facilities of the university. The I.M.S.L. subroutine DGEAR was used in integrating the governing ordinary equations. The subroutine finds approximations to the solution of a system of first order ordinary differential equations of the form y' = f(x,y) with initial conditions. There are two built-in approaches available to the user: the implicit Adams method and the backward differentiation procedure, also referred to as Gear's stiff method. For each of the procedures six different iteration schemes are available. The choice is based on the nature of the problem, storage needed, etc., and therefore requires some experimentation. The prime consideration is the stiffness of the problem. For reasons of efficiency and speed, the Adams method is used for nonstiff problems while with deployment and/or flexibility the stiff routine is used. An analytic Jacobian is supplied or calculated internally by the finite difference approach in case of linear and nonlinear problems, respectively. For every step, the DGEAR tests for the possibility that the step-size is too large to pass the error test (based on the specified tolerance), and if so adjusts it automatically. The step-size of integration specified by the user is employed only as a starting value and adjusted automatically by the subroutine. The user is advised to choose a starting step size which is considerably smaller than the avarage value expected for the problem. Values in the range of 1 0 -6 — 10 -8 were used here. The error which is controlled by way of the specified tolerance is an estimate of the local truncation error, that is, the error in a single step with the starting data regarded as exact. This should be distinguished from the global truncation error due to all steps taken to obtain y(x). The later is neither estimated nor controlled by the routine. The user manual for the subroutine advices that, if the problem is mathematically stable, 132 the global error at a given x should vary smoothly with the choosen tolerance limit in a monotonically increasing manner. Different tolerance levels are used here ranging over 1 0 - 6 — 1 0 - 9 depending on the problem. The numerical approach is now applied to study four different problems of contemporary interest in spacecraft dynamics, representing different levels of complexity, thus attesting to the versatility of the formulation and the computer code. The configurations considered are as follows: (a) The Space Shuttle represented as a triaxial rigid body. (b) The Space Shuttle based flexible cantilever beam including deployment. This has some similarity with the experiment, once proposed by Grumman Aerospace Corporation, aimed at construction of structural components in space. As pointed out in the introduction, two astronauts did attempt to integrate a beam in space during their extravehicular activity in 1985. (c) The Space Shuttle based deployment of two long beam type antennas creating a configuration similar to that of Wave In Space Plasma (WISP) study, a scientific project jointly under study by the U.S.A. and Canada. (d) The Space Shuttle based plate type flexible appendage which has some similarity with the SAFE (Solar Array Flight Experiment) which was jointly conducted by NASA and the Lockheed Missile and Space Company in 1984. The configurations are schematically shown in Figure 6-3. The following sections present response results for each of the configurations and associated discussion. The focus throughout is on the librational and vibrational dynamics. Their effect on the orbital motion being negligible for the size of the systems considered, the trajectory is defined by the classical Keplerian relations, i.e., the orbital and librational/vibrational motions are decoupled. Figure 6-3 Schematic diagrams of the configurations showing similarity with the cases studied: (a) Space Shuttle; (b) Orbiter based construction of beam as proposed by Grumman. Figure 6-3 Schematic diagrams of the configurations showing similarity with the cases studied: (c) WISP configuration jointly under study by the U.S.A. and Canada; (d) Solar array experiment conducted in 1984. 135 6.2 The Rigid Space Shuttle Dynamics For the Space Shuttle treated as a rigid body the generalized and specified velocities are f = (a p 7 ); s- T = ( Rc n); where a, /? and 7 are the shuttle libration angles, R defines the position of the Shuttle c center of mass, and 0 is the true anomally (0 = U). The inertia matrix for the Orbiter is taken as follows (Lagrange configuration; x and y along the maximum and the minimum inertias, respectively). 8646050 1 = -8135 1091430 Sym. 328108 27116 8286760 kg m . 2 The Orbiter is assumed to negotiate a circular orbit of 90.3 minutes period, which corresponds to a radius of about 6670 km from the center of the earth ( « 330 km altitude). Figure 6-4 shows librational response of the Orbiter in a circular orbit when subjected to a relatively small disturbance of 0.05° in roll, yaw and pitch simultaneously. The response is evaluated using nonlinear as well as linearized approaches for the Orbiter in three different configurations. Note, except for local details, particularly at large angles of attack, the linearized approach seems to predict the trend towards instability accurately (Figures 6-4a, b, c, d). During the small amplitude bounded motion, Figure 6-4e, linear and nonlinear analyses yield virtually the same response as expected. The Lagrange configuration representing the minimum moment of inertia axis along the local vertical and the maximum moment of inertia axis aligned with the orbit normal is stable as discussed in Section 5.6. Hence, from control consideration, the Lagrange configuration will be less demanding in terms of fuel expenditure. However, in the presence of a relatively large disturbance, the linear analysis would lead to misleading conclusions. This is clearly demonstrated through Figure 6-5. The Lagrange configuration found to be stable under small disturbances is now subjected to a 136 e =0 , 1 0 a M = £ [ o ] = y[0] = 0.05 d e g . 1 1 1 1 2 3 orbit 8 1 4 1 5 Figure 6-4 A comparison between linear and nonlinear responses to a small disturbance with the Orbiter in three different flight configurations. 137 Figure 6-5 Plots showing inadequacy of the linear analysis to accurately predict instability of the Orbiter. 138 roll, yaw and pitch disturbance of 4°. Note, the linear analysis continues to predict bounded motion (Figure 6-5a), while actually the system is unstable (Figure 6-5c). Note that the relatively large deviations from equilibrium in the yaw degree of freedom (Figure 6-5b), become unstable within five orbits with a slightly larger disturbance (Figure 6-5c). To get a better appreciation as to the system dynamics during transition to instability, the Lagrange configuration was subjected to pitch, yaw and roll disturbances separately (Figure 6-6). With a pitch disturbance as large as 30°, Figure 6-6a, the roll and yaw remain unexcited and the system is stable. The same is essentially true with a yaw disturbance as shown in Figure 6-6b. However, even with a reltively small roll disturbance, Figure 6-6c, the diverging yaw oscillations set-in tending towards instability. Thus roll control seems to be a key to ensure stability of the Orbiter in the Lagrange configuration. Figure 6-7, attempts to study the effect of roll control on the librational stability of the Orbiter in the Lagrange configuration. The spacecraft is in a circular orbit and is subjected to an initial disturbance in pitch as well as yaw of 4°, while the roll is assumed controlled, i.e., specified [67]. Here the generalized and specified coordinates are: f Two = (a /? ) ; f different deadband limits are used, i ax m = (i Rc A ). = 1° (Figure 6-7a, b) and i az m = 0.15° (Figure 6-7c, d) to have some appreciation as to the degree of control in roll needed to assure stability. The phase plane response in yaw is included to help judge the velocities involved. It is apparent that the roll control to the extent of 1° is not adequate and the Orbiter becomes unstable in yaw within five orbits. However, with the deadband limit of 0.15°, the system regains stability. It was suspected that such a demanding control for stability may be due to sharp positive peaks in the roll time history caused by firing of thrusters. Considerable extension in the deadband limits can be achieved through smoothing the peaks. This is shown in Figure 6-8, where the Orbiter's control strategy is improved to result in an approximately 139 Figure 6-6 Librational response of the Orbiter to an independent excitation in pitch, yaw and roll. Note the pitch and yaw disturbances lead to essentially uncoupled motions. The system appears to become unstable in yaw through its coupling with roll. 140 e = 0 , celO]-Q[0]- 4 deg. Figure 6-7 Librational response of the Orbiter in the Lagrange configuration with the time history of roll control actually used in practice. 141 sinusoidal roll time history with deadband limits of 6°. It also shows the effect of changing the roll control frequency. In Figures 6-8a,b the roll control frequency is taken to be 2 cycles/orbit which approximately coincides with the natural frequency in roll of the uncontrolled Orbiter in the Lagrange configuration (Figure 6-6). The system is unstable (Figure 6-8a, b), however, with the frequency increased to 6 cycles/orbit the system regains stability (Fig. 6-8c, d). Thus, an increase in the roll control frequency as well as reduction in the sharpness of the peaks in the roll time history appears to promote stability. Therefore at a frequency of around 13 cycles per orbit, which is normally used in the actual practice, the deadband limits can be further relaxed. 6.3 The Orbiter Based Beam Response As pointed out before, the system considered here has some similarity with the experiment carried out in 1985 by NASA, which was aimed at assessing the capability of constructing beam and truss-type structures in orbit. The experiment is idealized here as a uniform deployable beam fixed to the Orbiter. The fully deployed beam has the following parameters: mass length stiffness = 129 kg ; =33 m; = 436 Nm . 2 As the attitude dynamics of the Shuttle has already been studied in the previous section, the focus here is on the beam dynamics. The shuttle attitude now is assumed specified with time histories taken as zero, sinusoidal, or according to the controlled shuttle response. Several different orientations of the beam are considered, and the effect of deployment, orbit eccentriciy and beam stiffness studied. The system generalized and specified coordinates 142 Figure 6-8 Effect of smoothing the peaks and frequency of the roll control on the Orbiter's Librational response. Note, the deadband limits can be relaxed appreciably with the stability assured. 143 in this case are: A: = 1,2,3. Here SXj. and 5 . are the generalized coordinates associated with the j'-th mode of the y beam transverse vibrations in the x and y directions, respectively, and n is the number of modes considered. 6k represents the shuttle librational angles. Note, in the following discussion, we frequently refer to the beam-plane, which is defined by two coordinates: y and z. As pointed out before, the z coordinate is directed along the beam length, while the x and y coordinates define the directions of the beam lateral vibrations. The orbital plane is defined by the coordinates Y0 and Z0 (Figure 6-9). To start with the attention is focused on the vibratory response of the beam with the Orbiter's librations assumed absent, i.e., completely controlled. The beam is considered to be deployed and subjected to an initial disturbance corresponding to tip deflection of 4% of its length in the first mode. Figure 6-10 shows the effect of the beam's inclination, to the local vertical in the orbital plane, on the vibratory response to the above mentioned disturbance. In Figure 6-10a the tip is displaced out of the orbital plane. Note, the inplane motion remains unexcited in the first as well as the second modes leading to no tip motion in the orbital plane. The reverse is true in Figure 6-10b where the tip of the beam is subjected to the inplane disturbance. Thus for any arbitrary orientation of the beam with respect to the orbital plane the transverse oscillations are uncoupled. The orientation of the beam away from the local horizontal or vertical would result in a deflected equilibrium configuration in the orbital plane, as discussed before. The deflected equlibrium configuration acts as an inplane initial excitation. This is clearly shown in Figure 6-11 where the tip is subjected to an external disturbance of 4% out of the orbital plane. Generalized coordinates corresponding to the first and second modes as well as the tip deflection time histories are plotted over an orbit. Note, the small magnitude inplane vibrations as represented by 6Vl , 6y2 and 6y at the beam tip are shifted from zero 144 Z Y 0 X i X 0 Figure 6-9 A schematic diagram showing the orbital coordinates X0,Yo,Z0 coordinates x,y,z and beam with their origins at the instantaneous center of mass. In general the two coordinate systems are not coincident. 145 beam e=0 , L=0 . a=/9 = r=o (a) A j n = 0,90,180,270 Sxl (0)-0.02 I orbit orbit (b) Aj = n normal 30,-30 £ (0)=0.02/ y l \ \ \ \ \ \ \ \ I: *!', *•!; * !', f • ' . I; * I; 8 i; * I ; * " ;'!;;' " I \6yl 8xi i > s ? i'i ii ' i : ? « ? •' > >' S < ' J > >'»i < ' >> -i orbit orbit Figure 6-10 Response of the beam, located in the orbital plane, to a tip excitation in the first-mode of 4% of its length. Note the transverse oscillations S and 6 are x uncoupled. y 146 and occur around the equilibrium position, while the out-of-plane static deflection is zero. As expected the out-of-plane tip disturbance leads to a corresponding out-of-plane response as shown in Figure 6-10a. However, now there is a slight coupling between the inplane and out-of-plane motions, due to the Coriolis force directed along the beam, as indicated by the parametrically defined curve 6 vs. 6 at the beam tip (Lissajous like figure). Thus X y it may be concluded that the lateral vibration of a beam in the orbital plane is statically decoupled, however, slightly coupled dynamically through the Coriolis force along the beam length. Analysing Figures 6-10 and 6-11 in light of the analytical results of Section 5.6 suggests the following: • Neglecting the geometric stiffness due to orbital forces as compared to its material stiffness (as discussed in Section 5.6), the period of transverse vibration of a uniform beam , in the first mode, is 2 mi 3 Tr Thus the beam performs around 29.39 oscillations per orbit, in the first mode. This agrees quite well with the response data in Figures 6-10 and 6-11. • For A j n = 0 , 9 0 ° , 1 8 0 ° , 2 7 0 ° , the deflection in static equilibrium is zero (both in and out of the orbital plane), which is in agreement with the results of Section 5.6. The out-of-plane static deflection for a beam with A j n = 3 0 ° is also zero as expected. The inplane static deflections are around 6 yi = 0.9 x 10~ 2 m, and 6 y2 = —1.5 x 1 0 - 4 m for the first and second modes, respectively. From Section 5.6, for a beam in the orbital plane when fi* << 1 : - tx =0] n At ct = 6y = n 6 £4 m n 2 ^ - ^ 6 sa b ca s b = ~ 2 - 0 0 0 8 3 3 a b A ca 6 —60° (A = 3 0 ° ) , for the first and the second modes: 1 0 - 2 m; 6 y2 b • m 6 Vl = 2.815 x = —1.147 x 1 0 - 4 m . The minor discrepancies may be attributed 147 beam Xi n A e=0 , I --0 <$xl(0)-- 0 . 0 2 / , a = p = Y=o = 30 orbit orbit normal Figure 6-11 Plots showing the deflected equilibrium position of the beam acting as a small initial disturbance, in this case in the orbital plane. Due to an initial tip disturbance in the out-of-plane, the beam response is coupled as shown by the 6x -6y plot. 148 to axial forces caused by the Coriolis effect leading to stiffening and hence a change in equilibrium. • The base stiffness for the case under consideration is given by G 0 = 3 £ 3 m 1.51 N m. 2 fi2/A4 = The beam stiffness used here is about 380 times the base value, i.e., the beam is designed such that its maximum static tip deflection is about 0.263% of the beam length, or approximately 8.68 x 1 0 - 2 m. Figure 6-12 shows tip response of the beam for two different orientations in the plane defined by the local verical and the orbit normal, A o u t = 2 0 ° and 9 0 ° . Note now the two transverse motions, 8 and 8 , are completely coupled with the plane of vibration X y precessing, due to the Coriolis force perpendicular to the plane of oscillation. The uniform precessional rate is governed by the beam inclination angle \ t ou value at \ o u t and reaches a maximum = 9 0 ° . It may be pointed out that the plane of vibration of the beam precesses in one direction only (in this case clockwise) for a given \ 0 U f The effect of beam rigidity is studied in Figure 6-13, where a 9 0 ° out-of-plane beam is considered. The beam stiffness is taken as 436 Nm 2 and ten times this value in Figures 6-13a and b, respectively. From the response plots, the beam frequencies are approximately 30 and 95 cycles/orbit for the two cases. The ratio is around y/l6 as it should be. Figures 6-14 and 6-15 attempt to assess the effect of eccentricity on the tip deflection of the beam and the first two modes contributing to it. Two orientations of the beam are considered, A t n = 3 0 ° and A o u t = 9 0 ° , to facilitate comparison with earlier response data given in Figures 6-11 and 6-12a, respectively. The orbital eccentricity tends to impart a periodic increase in the stiffness, and hence the frequency of the beam reaches a minimum at the perigee and a maximum at the apogee. This is apparent in the tip time history as well as the modes. Such modulations of the frequency, although shown here for two orientations, were found to be present at all A. The frequency condensation effect is related to 1 + ecos# term corresponding to the geometric stiffness associated with orbital rate Omega, appearing in the denominator of the stiffness terms. Furthermore, it is interesting 149 orbit orbit Figure 6-12 Response of the beam to a tip disturbance when located out of the orbital plane. The transverse motions 6 and S are coupled and the vibration plane X Y of the beam tip precesses at a uniform rate. 150 ure 6-13 Response plots showing the effect of a change in the beam flexural rigidity. 151 vertical 16 Figure 6-14 Effect of the orbital eccentricity on vibratory response of the tip for two orientations of the beam. Note the frequency condensation effect is particularly pronounced near the apogee. 152 vertical Figure 6-15 Response plots showing the frequency condensation effect in the first two modes due to an orbital eccentricity of 0.2. 153 to observe, for the planar configuration of the beam (A t n = 30°, Figures 6-14a, 6-15a), modulation of the inplane vibrations (8 , S , 8 ) around the changing equilibrium position y yi y2 of the beam in orbit. The effect of beam deployment on the tip dynamics is studied in Figures 6-16, 617, and 6-18. The initial tip deflection is the same as before. Two cases are considered: sinusoidal time history with different duration of deployment (Figure 6-16), and different time histories with the same duration (Figure 6-17). A comparison between the deployed and deployment cases is considered in Figure 6-18. As can be expected, the frequency of oscillation in the out-of-plane motion gradually decreases with deployment finally attaining a steady state value when fully deployed (Figure 6-16). On the other hand, the beam shows no inplane oscillation during deployment. Instead it drifts in the orbital plane due to the Coriolis effect which, as can be expected, is proportional to the speed of deployment. Note, the out-of-plane oscillations are governed by the initial condition and remain unaffected by the Coriolis effect. Hence they reach the same steady state amplitude, although it is much larger during deployment compared to the deployed case (Figure 6-11). However, this is not the case for inplane oscillations where the higher Coriolis acceleration (case a) leads to a higher steady state amplitude. Figure 6-17 studies response for deployment normal to the orbital plane (X ut = 90°) 0 under two different strategies extending over the same duration of deployment (1/2 orbit). As expected, a reduction in frequency during deployment and coupled charecter of the transverse oscillations associated with an increase in tip amplitude are noticed. It is of interest to recognize that the steady state amplitude here remains unaffected by the strategy for a given time of deployment. As in the previous case, Figure 6-16, the steady state amplitude attained subsequent to the deployment is substantially larger than that in the deployed condition (Figure 6-12a). A direct comparison of the response for a deployed and deploying beam to the same initial conditions is shown in Figure 6-18. The beam is 154 vertical Figure 6-16 Effect of the duration of deployment on the tip response of the beam located in the orbital plane. Note, the beam-tip drifts in the orbital plane due to the Coriolis force. In general, the steady state amplitude of oscillations are higher for the deployment case. 155 Figure 6-17 Effect of deployment strategies on tip response of a beam located normal to the orbital plane. Note a reduction in beam frequency during deployment. The steady state amplitude is essentially independent of the strategy for a given time of deployment. 156 Figure 6-18 A comparison between the response of a deploying and deployed beam, initially aligned with the orbit normal, to a tip disturbance. 157 oriented at 90° to the orbital plane and duration of the deployment is taken to be 0.1 of the orbital period. A s before, the deployment appears to magnify the transient and steady state response. In practice the Orbiter's libration will be controlled to a specified tolerance limit. A typical time history of the controlled Space Shuttle librations, in the Lagrangian configuration, during an orbit is shown in Figure 6-19 [67]. Figure 6-20 shows resonance of the beam due to forced vibration. Figures 6-21 to 6-23 focus attention on response of the deployed beam during such forced excitation of the Orbiter. Figure 6-24 is a comparison between the response of two beams, with different orientations, in an eccentric orbit. Finally, Figures 6-25 and 6-26 compare deployment effects with response of the deployed beam in presence of the excitations caused by thruster firings to control the Orbiter's librations. In Figure 6-20 the frequency of the Orbiter's controlled roll is modified to 30 instead of the normal 12 cycles/orbit (Figures 6-19). The new forcing frequncy of the roll motion is close to the beam's transverse vibration frequency, which is 29.4 cycles/orbit. The beam is oriented along the orbit normal ( A o u t = 90°), hence the inplane and out-of-plane vibrations are dynamically coupled. The inplane beam vibration is controlled by the Orbiter's high frequency roll motion. Even in absence of an initial tip deflection, the beam vibration amplitude continues to grow even after one complete orbit. This clearly demonstrates a close to resonance phenomena. Figure 6-21 shows the tip response as well as the first two modes contributing to it for a beam deployed along the local horizontal in the orbital plane ( A t n = 90°). A t the outset it should be recognized that, for this inplane configuration of the beam, the out-of-plane motion 8 and inplane response 8 are decoupled as seen before (Figure 6-10). X y Hence one would expect the Orbiter's yaw motion (0) to be reflected in the out of plane 8 X motion, while the Orbiter's pitch motion (a) would govern the inplane response (8 ). The y beam response precisely reveals these characteristic features. The 8 vs. 8 plot further X y emphasizes the tip vibration and drift caused by the response to pitch-yaw excitations. 158 local vertical & (yaw) local horizontal /(roll) orbit a (pitch) orbit CO normal X(roll) a (pitch) /3 (yaw) oo 0 Figure 6-19 orbit Representative controlled motion of the Orbiter during a typical orbit. 159 orbit normal o 00 Figure 6-20 Forced response of a beam to the Shuttle librational control. Note the beam frequency is close to that of the roll disturbance. 160 a, /?, 7 specified (Fig. 6-19) e.O,£--o,Xjn= 9 0 ° 0 I -\ ~ Z O x 0 t> i t orbit normal r 1 1 4 r A tip.m x1(H Figure 6-21 Forced response of a beam extended along the local vertical from the Orbiter in the Lagrange configuration. Note, the beam response closely follows the yaw and pitch excitations with modulations at the natural frequencies in different modes superimposed. 161 i H 1 -4 S x 1 0 tip.m x 1 0 1 1 4 3 Figure 6-22 Effect of an orbital eccentricity on the forced beam response of the case considered in Figure 6-21 (beam along the local horizontal). The inplane tip response 6 ^ y increases by an order of magintude even in the presence of relatively small smooth disturbances. still present . The frequency condensation effect is 162 ;ure 6-23 Forced response of a beam, aligned with the orbit normal. Note the effect of the dominant roll excitation at 13 cycles per orbit which is particularly clear in the second mode. 163 vertical Figure 6-24 A comparison between the forced response of a beam, located either along the local horizontal or the orbit normal, with the Shuttle in an orbit of e = 0.2. 164 Figure 6-25 A comparison between the forced response of a beam, aligned with the local horizontal, showing the effect of deployment. 165 £ tip.m x 1 0 3 <$ tip,m x10 ; x Figure 6-26 Effect of deployment on the forced response of a beam located along the orbit normal. 166 The beam being aligned with the roll axis the corresponding disturbance has relatively little effect. The presence of eccentricity (Figure 6-22) does not seem to affect the dominant character of the response as seen in Figure 6-21. However, the inplane tip response 6 is y now increased by an order of magnitude due to the eccentricity induced excitation. The frequency condensation behaviour discussed earlier, Figures 6-14, 6-15, is still retained. Forced response of the beam when oriented normal to the orbital plane (X t ou = 90°), with the Orbiter in the Lagrange configuration, is shown in Figure 6-23. As seen before, Figure 6-12, the beam response in S and 6 for the out-of-plane configuration is coupled. y x Hence the roll and yaw excite both 8 and 8 as evident from the S vs. 6 plot. However, X y x y the roll disturbance, being at a higher frequency and hence with a higher acceleration, appears to be dominant as apparent from the amplitude modulation of the response at the roll frequency (around 13 cycles per orbit). Note that all the modes vibrate at their own characteristic frequencies superimposed on the frequency of the Orbiter's attitude motion. Figure 6-24 compares response of two beams, located in and out-of-the orbital plane, at 90° to the local vertical. The Orbiter is in an eccentric orbit and the excitation is due to the Orbiter's librational control (Figure 6-19). As in the unforced cases (Figures 6-14 and 6-15), the orbit eccentricity tends to impart a periodic increase in the stiffness. The lightly coupled behaviour of the inplane beam transverse vibration as compared to the coupled motions in and out of the orbital plane is still clearly evident for this forced motion. The effect of combined deployment and base forcing of a beam as compared to the forced case alone is shown in Figures 6-25 and 6-26. The beams are oriented at 90° in and out of the orbital plane. Even in the presence of forcing, the deployment is able to suppress the beam lateral vibrations. When the beam is completely deployed, the response tends to be similar to that for the deployed case. It was thought useful to study the effect of the number of modes considered to represent the beam dynamics. To this end the Orbiter was taken in the Lagrange configu- 167 ration with the beam normal to the orbital plane (A t = 90°). The 4% tip deflection 8 in X ou the first mode was applied with the Orbiter free of librational motion. In absence of any external forces one would expect the beam to vibrate in the first mode alone. However, due to orbital motion giving rise to centrifugal, Coriolis and other forces, the beam deflection deviates from the first mode requiring higher modes to evaluate its dynamics. To assess contribution of the higher modes, the beam response was analysed using the first, the first two and the first four modes. The results are presented in Figure 6-27. Two important conclusions can be drawn: (i) time history of the beam deflection remains essentially unaffected by the number of modes used; (ii) amplitudes of the generalized coordinates associated with the first mode {8 ,6 ) Xl 6 ,8y X2 2 yi are around three orders of magnitude higher than of the second mode and, of course, still larger compared to the higher modes. In the case of the beam excited by the controlled librations of the Orbiter, Figure 6-19, the dominant mode will be governed by the resonance condition. It is apparant that over the range of frequency of the forcing functions the first mode will have the major contribution (Figures 6-20,6-26). A word concerning the accuracy of the results presented here would be appropriate. To this end, permissible error during numerical integration of the governing nondimensional differential equations was varied systematically to assure reliability of the results. The case presented in Figure 6-23 was studied using error tolerance levels of 1 0 , 1 0 - 4 - 6 and 10~ . 8 The results are shown in Figure 6-28. Both components of tip deflections and generalized coordinates are compared. It is apparent that the tolerence level of 1 0 -4 is inadequate and gives misleading response. However, the results obtained using the tolerance levels of 10 -6 and 1 0 -8 tolerance of 1 0 6.4 are essentially the same. Hence during the numerical integration an error -8 or lower was used. Dynamics of the W I S P Type Configuration Here a configuration similar to the WISP system is considered. The Shuttle is in 168 (a) first mode orbit orbit (b) first two modes Figure 6-27 A typical case of beam response showing the effect of mode truncation. It appears that the higher modes have very little effect on the time history of the beam deflection: (a) first mode ; (b) first two modes ; 169 orbit orbit (c) first four modes Figure 6-27 A typical case of beam response showing the effect of mode truncation. It appears that the higher modes have very little effect on the time history of the beam deflection: (c) first four modes . Figure 6-28 Plots showing the influence of a permissible error during numerical integration. The allowable error of 1 0 - 8 or less resulted in reliable data. 171 the Lagrange configuration with two beams extending in the orbital plane, in and opposite to the direction of the orbit tangent. The system parameters are as follows : • beam mass = 40.25 kg • beam length = 150 m • beam rigidity = 1676 Nm 2 • the Shuttle inertia matrix = 9608110 0 0 0 0 1227612 0 0 9204755 kg m? In the following, the system dynamics is studied by introducingflexibility,deployment, orbit eccentricity, etc., progressively to help idenitfy their distinctive contributions. In particular, the system is first considered as rigid and the effects of orbit eccentricity as well as deployment on the attitude dynamics examined. This is followed by the study where booms are treated as flexible beams, and the libration-vibration interaction dynamics analysed for both circular and elliptic orbits. 6.4.1 System treated as rigid To start with, the effect of orbit eccentricity on the WISP configuration, treated as rigid, is studied. For the purpose of comparison two cases are considered: (a) the Shuttle without beams (Figure 6-29); and (b) the Shuttle with two rigid beams fully deployed (WISP configuration, Figure 6-30). In Figures 6-29a and 6-29b, the Shuttle orbit has an eccentricity of e=0 and 0.2, respectively. The WISP configuration is considered in Figure 6-30 with three different values of eccentricity: e=0, 0.1 and 0.2. In all the cases the system is given an initial disturbance of 12° in pitch. The plots show that, in general, an eccentric orbit imposes a periodic force to which the system responds in pitch depending on its inertia parameters. Figures 6-29b and 6-30b display an amplitude modulation in the pitch response, while in Figure 6-30c the pitch motion is driven to a magnitude higher than 90°, causing the shuttle to end up rotating about the orbit normal. Note, for all the cases considered in Figures 6-29 and 6-30, the roll and yaw remain unexcited. Thus during 172 a(0) = 12°, /J(O) = 7(0) = 0 / i i ' I i » \ t i \ I / > t i i I I \ \ / \ / (a) e=0 o O. ~ \ hi: <=>' o ( b ) e=0.2 o.o 0.6 n 1.6 i 2.4 3.2 4.0 orbits Figure 6-29 A comparison showing the effect of orbital eccentricity on the libration response of the Orbiter in the Lagrange configuration with an initial pitch disturbance. 173 1 a(0) = 12°, m = 7(0) = 0, 1 = 0 orbit o. CM ,' > , o o" ' 1 ' cf / « ' 1 ' \ ' , v ; , / 1 \ x , > ' / X X > > 1 V (a) e=0 o I i i i i /' * \\ o o" ' / . / ' s X \ / w \ ' // < 1 1 * \ /• > \ / x X \ 1 X » / / ' ' ' (b) e=0.1 0.0 1 ] 1 1 o.e 1.6 2.4 3.2 4.0 orbit Figure 6-30 A comparison showing the effect of orbital eccentricity on the libration response of the WISP system, treated as rigid and in the Lagrangian orientation, with an initial pitch disturbance. 174 pure pitch excitation, the pitch response is uncoupled from roll and yaw. To help understand the results, the forces associated with deployment maneuvers are discussed first. It should be noted that reaction forces due to deployment cancel in the cases under consideration because of the system symmetry. In absence of any initial disturbances in yaw and roll, the system responds in pitch to the exciting torques, which in this case, lie in the orbital plane. During the deployment of the two beams, three types of forces are involved: • Conventional orbital forces (centrifugal and gravity), which in the case of the Lagrange configuration act to stabilize a as long as the pitch angle is less than 90°. For 180° > a > 90°, the torque on the system switches its direction towards another equilibrium position in which the system is upside down. • Coriolis forces due to the system rotation and deployment of the beams lead to a torque along the orbit normal with direction depending on the the involved motions. Note, the system rotation is the resultant of the orbital and librational motions (the pitch rate in this case). The magnitude of the torque is a function of the angular and deployment rates as well as the beam mass, i.e, the length of the beam at a given instant. • Forces due to rate of change of inertia during deployment. Figure 6-31 studies the effect of deployment velocity on the stability of the pitch motion. Here the two beams are deployed, each extending to a final length of 140 m, using a sinusoidal velocity history. Note, as the beam length is fixed (140 m) and the deployment time as well as the time history specified (sinusoidal), the peak deployment velocity in each case is different; however, the areas under the velocity curves are the same (Figures 6-31a, b and c). The system is subjected to an initial pitch disturbance of 12° as before. The deployment periods are taken to be 0.4, 0.8 and 1.5 orbit in Figure 6-31. In all the cases, the pitch response indicates instability. In the first case (Figure 6-31a) the system ends up rotating about the orbit normal while in the other two cases (Figures 175 6-31b,c) it undergoes a steady state oscillation about a new equilibrium position (where the system is upside down) with an amplitude governed by the conditions at termination of the deployment. Numerical experiments indicated that the permissible computational error is a crucial factor in obtaining reliable data during deployment. Here, the specified numerical error was taken to be 1 0 . A higher value of allowable error often predicted, -8 incorrectly, stability. This is attributed to the delicate balance between extremely small forces involved, particularly during the transient state as the pitch angle approaches zero. Next, to study the effect of the extended mass, the beams are deployed for 25 m using the same deployment history and duration but with initial lengths of zero (Figure 6-32a) and 100 m (Figures 6-32b). As expected, the impulse associated with deployment being proportional to the deployed mass, amplitude of the steady pitch oscillations for the larger initial beam length is higher (60° against 40°). Note, the similar trends are observed in Figures 6-33a and 6-33b, where 50 m of beam is deployed from initial lengths of zero and 100 m, respectively. In fact, in the latter case the system becomes unstable due to a larger initial deployed length resulting in a higher impulse and a corresponding higher torque in pitch. A comparison between Figures 6-32b and 6-33a where 25 and 50 m of beam lengths are deployed from a starting length of 100 meters and zero, respectively, indicates a larger impulse associated with the first case. Inspite of the higher deployment velocities associated with the second case, the resulting impulse in the first case is larger because of the larger deployed mass. The effect of deployment time history is studied in Figure 6-34a, where the beam is deployed to 22 m using a sinusoidal displacement function. Similarity with the results in Figure 6-32a suggests that the character of the deployment time history has virtually no effect on the response. Thus the parameter of importance is the impulse associated with the deployment which, in turn, is governed by the initial deployed length, deployment velocity and the duration of deployment. 176 Figure 6-31 Effect of deployment velocity on the librational response of the WISP configuration, treated as rigid and in the Lagrangian orientation, with an initial pitch disturbance. Duration of deployment: (a) 0.4 orbit . 177 Figure 6-31 Effect of deployment velocity on the librational response of the WISP configuration, treated as rigid and in the Lagrangian orientation, with an initial pitch disturbance. Duration of deployment: (b) 0.8 orbit . 178 Figure 6-31 Effect of deployment velocity on the librational response of the WISP configuration, treated as rigid and in the Lagrangian orientation, with an initial pitch disturbance. Duration of deployment: (c) 1.5 orbit . a(0) = 12°, 0(0) = 7(0) = e = 0 o a: faO »N / O 0.0 0.4 n 0.8 orbits • r— 1.2 l .6 2.0 Figure 6-32 Effect of the initial deployed mass on librational response of the WISP configuration with an identical deployment time history and duration: (a) £(0) = 0. 2.0 orbits a(0) = 12°, 0(0) = 7(0) = e = 0 CD . I 80.0 0.0 80.0 deg. -a- o Figure 6-32 Effect of the initial deployed mass on librational response of the WISP configuration with an identical deployment time history and duration: (b) £(0) = 100 m. O O 0.0 0.4 2.0 orbits 0.0 0.4 2.0 00 o c*(0) = 12°, 0(0) = ( 0 ) 7 = e = 0 00 Figure 6-33 Influence of the initial deployed length on librational response of the WISP configuration. Note, the deployment time history and duration are identical in the two cases: (a) deployment from zero to 50 meters. -SJ0.0 0.4 2.0 orbits 0.0 2.0. a(0) = 12°, (3{0) = ( 0 ) = e = 0 7 o bb Figure 6-33 Influence of the initial deployed length on librational response of the WISP configuration. Note, the deployment time history and duration are identical in the two cases: (b) deployment from 100 to 150 meters. 0.0 0.2 1.0 orbits 0.0 0.2 1.0 a(0) = 12°, = (fj) 7 = e = 0 CD O O TO o o" 03 O 0.0 -40.0 a 40.0 deg. Figure 6-34 Effect of sinusoidal deployment history on the librational response of the WISP configuration: (a) c*(0) = 12°. Note, a comparison with Figure 6-32a shows that the time history of deployment has virtually no effect on pitch response. a{0) = -12°, 0.0 0.4 0(0) = 7(0) = e = 0 1.2 0.8 1 .6 ORBIT 50.0 2.0 orbits Figure 6-34 Effect of sinusoidal deployment history on the librational response of the WISP configuration: (b) a(0) = ^0 2.0 orbits -12°. 0.0 0.4 2.0 185 A similar case is considered in Figure 6-34b but with a negative initial pitch disturbance. Here, the initial torque tends to move the system towards equilibrium briefly, however, the steady pitch response is slightly larger than that for the case of positive pitch disturbance (Figure 6-34a). Thus the initial system configuration for the same magnitude disturbance affects the steady state response. 6.4.2 System treated as flexible As a second step, towards a more realistic model, the beams are considered to be flexible. The objective here is to examine the libration-vibration interaction dynamics for the system in circular and eccentric orbits. Figure 6-35a shows the system response to a disturbance in the form of inplane symmetric deflection of the two beams, while keeping zero librational initial condition. With only the first mode considered, the two beams continued to vibrate symmetrically at their own natural frequency. As expected, the libration degrees are not excited by such a symmetric mode of motion. Actually a pitch response of an order of 10~ 15 was recorded, which was contributed by numerical errors, as confirmed through additional computational runs (not shown). This also indicates the high level of accuracy demanded during the numerical computations to assure precise information. Figure 6-35b shows the response for the same situation, however, with e = 0.2. As before, with the eccentricity induced excitation, the system responds in pitch only with the vibrations confined to the orbital plane. Thus the coupling between inplane and out-of-plane degrees of freedom continue to be absent as long as the disturbance is in the orbital plane. Of particular interest is the frequency modulation of the response due to variation in the geometric stiffness of the beam along the orbital path as explained earlier. In Figure 6-36 the two beams are initially deflected asymmetrically in the orbital plane in the first mode with the librational disturbance taken to be zero. As can be expected, even for e = 0 now the pitch motion is excited and modulated by the beam 186 beam 1 : ^(0) = 0, S {0) = 0.0146£ , Vl beam 2 : 6 (0) = 0, 6 (0) X = 0.0146£ , yi orbit a(0) = 0(0) = 7(0) = 1 = e = 0 . CO o / A - / X \ \ ^ / to \ / N / o o" * * / / "V * \ / \ \ / \ v _ . V / \ • \ \ \ \ -T \ INJ beam 1 o IT) ' o CL, O" •I -H beam 2 \ AAA AAA o in. 00 t- J V \! 0.2 r MM] j V M l 1 'i w i V 0.6 0.4 0.8 "v" 1.0 orbit Figure 6-35 Libration-vibration response of the WISP system, in the Lagrangian configuration, to a disturbance in the form of inplane symmetric deflection in the first beam mode: (a) e=0. 60.0 ^ A ^ f\ i 'i ii 'ii A ^ A A /A A Mi; A ii i A A A ' 1 o o' f A 1f A ( 1 A orbit beam 1 : 6 {0) = 0, <5 (0) = 0.0657£ , 1 X o beam 2 beam 2 : 6 {0) = 0, 6 {0) = 0.0657£ , X a I yi (0) = 0(0) = 7(0) = t = 0, e = 0.2 . A ,1 /i Ou <=>' o 0.0 yi 0.4 —j !— 0.8 orbits 1.2 1 .6 2.0 Figure 6-35 Libration-vibration response of the WISP system, in the Lagrangian configuration, to a disturbance in the form of inplane symmetric deflection in the first beam mode: (b) e=0.2 . 188 frequency while the two beams continue to oscillate according to the initial asymmetric condition (Figure 6-36a). Response of the beams in the first two modes is presented in Figure 6-36b. As is usually the case, the dominant contribution to the beam motion is from the first mode where the generalized coordinate exhibits its characteristic frequency. The second mode responds at its frequency which is modulated by the first mode frequency, and its contribution is several orders of magnitude smaller. Figure 6-37 focuses on the effect of orbital eccentricity with the WISP configuration subjected to the initial conditions same as above. Since the pitch response to the eccentricity induced excitation is dominant, the modulations due to the asymmetric motions of the beams are not visible in Figure 6-37a compared to those in Figure 6-36a. On the contrary, the beam response is clearly modulated by the Shuttle's pitch motion. Note, the modulations are reflected in both the beam modes as shown in Figure 6-37b. The frequency condensation effect due to orbital eccentricity can also be discerned. So far the attention was directed towards disturbances in the plane of the orbit. The next step would be to consider out-of-plane excitations. Figure 6-38 shows response of the system to such a disturbance with beams deflected symmetrically in the first mode out of the orbital plane. Obviously, to this initial condition, the system libration and the beam inplane vibration should not respond. Hence, the two beams continue to oscillate symmetrically, while the Shuttle remains unexcited (not shown). Next, the out-of-plane disturbance is introduced (symmetrically) once again but in the second mode (Figure 639). The WISP antennas (beams) predominantly respond in the second mode with an insignificant contribution from the first mode, which is now modulated at the second harmonic. Note also that the system response is dominated by the elastic force as against the inertia or gravity effects. The effect of orbital eccentricity is studied in Figure 6-40. The system is subjected to an out-of-plane symmetric disturbance in the first mode with the tip deflections as specified on the diagrams. Two values of eccentricity are considered: e = 0.1; 0.2. As can 189 beam 1 : ^(0) = 0, 6 (0) = - 0 . 0 1 4 6 £ , yi beam 2 : 6 (0) = 0, 6 (0) = + 0 . 0 1 4 6 £ , X yi a(0) = /?(0) = 7(0) = l = e = 0 . orbit Figure 6-36 Libration-vibration response of the WISP system to an inplane asymmetric deflection in the first beam mode: (a) libration and beam tip response . beam 1 r i 0.0 1 0.2 1 0.4 1 0.6 orbit 1 0.8 1 1.0 o i beam 2 i 0.0 —i 0.2 i i 0.4 0.6 1 0.8 1 1.0 orbit Figure 6-36 Libration-vibration response of the WISP system to an inplane asymmetric deflection in the first beam mode: (b) beam-mode response. to o / \ \ _» / / / / —. s \ \ \ \ / / beam 1 -100.0 60.0 beam 1 : 6X(0) = 0, 6yi{0) = -0.0146£ , beam 2 0.8 orbits 1.2 1 .6 2.0 beam 2 : ^(0) = 0, 6yi(0) = +0.0146£ , a{0) = /?(0) = 7(0) = t = 0, e = 0.2 . Figure 6-37 Effect of an orbit eccentricity on the libration-vibration response of the WISP system in the Lagrangian configuration. The system is subjected to an initial inplane asymmetric deflection in the beam's first mode : (a) libration and beam tip response. o beam 1 orbits beam 2 0 orbits Figure 6-37 Effect of an orbit eccentricity on the libration-vibration response of the WISP system in the Lagrangian configuration. The system is subjected to an initial inplane asymmetric deflection in the beam's first mode : (b) beam modal response. 193 Figure 6-38 Vibratory response of the WISP system in the Lagrangian configuration to an out-of-plane symmetric excitation in the beam's first mode. The Shuttle remained unexcited in the librational mode (not shown) . beam 2 if) I 0 Figure 6-39 orbit Vibration response of the WISP system in the Lagrangian configuration to an out-of-plane symmetric deflection in the beam's second mode. The beam 1 : 6 (0) = - 0 . 0 0 7 3 £ , 6 (0) = 0 , librational response was found to be zero beam 2 : 6 (0) = + 0 . 0 0 7 3 £ , 6 {0) = 0 , (not shown). a(0) = 0(0) = -y(O) = 1 = e = 0 . X2 X2 y y 195 a(0) = 0(0) = 7(0) = I = 0 . beam 1 and 2 : symmetric initial deflection in Sx , 8y = 0 deg. LO \ *- " O cf / \ O CM _ 1 \ \ \ ^ 100.0 1 ?~ / \ /' e=0.2 0.0 0.4 0.8 1.2 orbits Figure 6-40 A comparison showing the effect of orbital eccentricity on the librationvibration response of the WISP system, in the Lagrange configuration, when subjected to an initial symmetric deflection (first mode) in the out-of-plane direction : (a) librational response for the orbital eccentricities of 0.1 and 0.2. 196 beam 1 : S (0) = +0.0073£, S (0) = 0 , Xl y beam 2 : 6 (0) = -0.0073£, 6 (0) = 0 , Xl y a(0) = /3(0) =' -/(O) = I = 0, e = 0.2 . beam 1 6 tip X m beam 2 £ tip x m Figure 6-40 A comparison showing the effect of orbital eccentricity on the librationvibration response of the WISP system, in the Lagrange configuration, when subjected to an initial symmetric deflection (first mode) in the out-of-plane direction : (b) vibrational response for the orbital eccentricity of 0.2 . 197 beam 1 : SXl (0) = +0.0657£, 6y (0) = 0 , beam 2 : 6Xl {0) = -0.0657£, 6y (0) = 0 , a(0) = P{0) = i(0) = £ = 0, e = 0.1 . beam 1 <5 tip x m beam 2 6 tip X m Figure 6-40 A comparison showing the effect of orbital eccentricity on the librationvibration response of the WISP system, in the Lagrange configuration, when subjected to an initial symmetric deflection (first mode) in the out-of-plane direction : (c) vibrational response for the orbital eccentricity of 0.1 . 198 be expected, the attitude response is essentially of the same character as observed earlier for corresponding rigid andflexiblecases (Figure 6-40a). Note, the orbital eccentricity induced forcing excites the pitch motion, which in turn causes the inplane beam vibration (6 ). y Therefore, the latter follows the character of the pitch response quite closely, while the outof-plane motion (6 ) proceeds with the frequency condensation effect due to eccentricity. X A comparison between the results in Figures 6-40b and 6-40c indicates that the beam vibrations in and out of the orbital plane are decoupled, as seen earlier in Figures 6-10 and 6-11. It may be pointed out that a large tip deflection of 6.6% is taken in Figure 6-40c to emphasize this behaviour. The small inplane motion in Figure 6-40c as compared to that in Figure 6-40b indicates that the inplane vibrations are governed by the orbit eccentricity. Response to an impulsive pitch disturbance, corresponding to a torque along the orbit normal, is presented in Figure 6-41. Such impulsive disturbances can arise during asymmetric firing of the on board thrusters, release or retrieval of payloads, slewing maneuvers of the CANADA A R M , and a variety of other sources. The Shuttle responds at its own pitch librational frequency while the beams vibrate asymmetrically at their characteristic frequencies modulated by the pitch response. Response of the system to an initial pitch disturbance of 12°, while negotiating a circular or an elliptic trajectory, is studied in Figure 6-42. Here, as in the previous case (Figure 6-41), the Shuttle responds at its own frequency, while the beams vibrate (the first two modes) in an asymmetric fashion at frequencies modulated by the pitch motion. An initial yaw disturbance of 12° is considered in Figure 6-43. The libration response given in Figure 6-43a is essentially the same as that for the rigid case studied earliear (Figure 6-6b). Note, through coupling, the small amount of pitch and roll motions are also excited which, in turn, lead to both the in and out-of-plane vibrations of the beams. The beam response, as before, is at its characteristic frequency modulated by the Shuttle's librations in pitch and yaw (Figure 6-43b). The roll axis being initially aligned with the beam suggests a negligible effect of the roll on all the beam vibrations. The Lissajous- d(0) = 1.5 x 1 0 - 2 deg/s , A y beam 2 : 6 (0) = 0, 6 {0) = 0 , X 1 " ^--__L~--^ 1 beam 1 : ^(O) = 0, 6 {0) = 0 , O * y a(0) = /?(0) = -y(0) = t = e = 0 . CO \ bi a/ \ • \ t \ \ / \ / / \ r- ° \ / \ o \ fin / \ \ / V / \ / \ \ o 00 i / s i / i i beam 1 H. / 1 1 1 8 X 1 beam 2 / \ 0.0 0.2 0.4 0.6 S x 0.8 1 .0 orbit Figure 6-41 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an impulse along the orbit normal : (a) libration and beam tip response. Figure 6-41 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an impulse along the orbit normal : (b) beam modal response. K5 o o 201 a(0) = 12°, 0(0) = 7(0) = I = t = 0 beam 1 and 2 : 6 {0) = «„(0) = 0 . X o beam 1 A A a ° \ /*\ A / / \ / \ 1 1 v I \/ \ I \ J \J . r- . A / \/ V \' V v V i i i i beam 2 o A / \u/ \ / \ A A" \ / v " / \ 1 \ / \ o 0.0 i 0.2 V i 0.4 0.6 • V\/\ \/ i 0.8 1.0 orbit Figure 6-42 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an intial pitch disturbance: (a) libration and beam tip response, e=0. orbit .0 o r b i t Figure 6-42 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an intial pitch disturbance: (b) beam modal response, e=0. bib <u ^ ° \ r—S \ \ „' a t f 1 1 1 J. \ \ \ / / / V / 1 v ' 1 i 1 beam 1 o 1 o in i i i 40.0 a(0) = 1 2 ° , e = 0.2 , 0(0) = 7 (0) = t = 0 , beam 1 and 2 : S {0) = 6 {0) = 0 t y Figure 6-42 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an intial pitch distur0.0 0.8 1.6 orbits 2.4 3.2 4.0 bance: (c) libration and beam tip re- sponse, e=0.2 . P(0) = 12°, a{0) = 7(0) = I = e = a deg 0 deg 7 deg. Figure 6-43 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an initial yaw disturbance: (a) libration time history . 205 Figure 6-43 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an initial yaw disturbance: (b) beam tip response . 206 like representation shows that the tips of the beams move in opposite directions at all instants. This is expected, since the inertial force on the beams, due to system librations, is dominant in this case. Guided by the results of Figure 6-6, two initial roll disturbances (7 = 3° and 9°) are considered in Figures 6-44 and 6-45, respectively. With 7(0) = 3° (Figure 6-44), the system is stable both in librational and vibrational degrees of freedom. However, with an increase in roll disturbance to 9°, the Shuttle becomes unstable in yaw through coupling (Figure 6-45). Note, the growth of yaw is relatively slow here as compared to the case of the Orbiter considered in Figure 6-6. This is because of the change in the ratio of roll-yaw moments of inertia brought about by the presence of the two beams. The beams respond in asymmetric fashion as expected with the amplitude of motion higher by an order of magnitude due to librational instability. 6.5 The Orbiter Based Plate Dynamics The NASA sponsored Solar Array Flight Experiment (SAFE, Figure 6-3d) was con- ducted in 1984 to help develop the capability of on-orbit extension/retraction of the light weight, deployable solar panel. It was also aimed at measurement of vibration modes and structural damping of the solar panel using video cameras and laser ranging instrumentation. A solar array, around 31 m in length and capable of generating 13.7 kW of power, was used in the experiment. The array can serve to augment the Orbiter's power supply. An idealized representation of the experiment may be a flexible plate-type appendage assumed fixed to the Orbiter. In this preliminary study the attitude motion is linearized about the nominal equilibrium position, i.e., the system librations are considered to be small. The plate mass is taken to be 100 kg, and its dimensions are the same as in the SAFE expriment: 31 m length and 4 m width. It was located at the Orbiter cm. The Orbiter's moments of inertia are identical to those used in Section 6.2. In Figures 6-46 and 6-47 the flexible plate is deployed along the maximum moment o a deg. 0 deg. ^ deg. Figure 6-44 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an intial roll disturbance of 3 ° : (a) libration time history. to o ~4 208 6 , tip mxlO £*tip mxiu Figure 6-44 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an intial roll disturbance of 3° : (b) beam tip response. o to bb (0) = 9 ° , a(0) = 0(0) = £ = beam 1 and 2 : Sx{0) = 6y{0) 7 CD O o i -20.0 20.0 0.0 « d e S- ' -200 200 p deg. 600 I -20.0 0.0 7 deg. Figure 6-45 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an initial roll disturbance of 9° : (a) libration time history. Note the instability in yaw through a coupling with roll. 20 210 beam 1 6 X t i p m beam 2 $ x t i p m Figure 6-45 Libration-vibration response of the WISP system, in the Lagrangian configuration, to an initial roll disturbance of 9 ° : (b) beam tip response. 211 Figure 6-46 Effect of the Orbiter's flight configuration on the transient librational response during deployment of a flexible plate-type appendage along the maximum moment of inertia axis. 212 of inertia axis during three different orientations of the Orbiter. In Figure 6-46 the attitude response of the Space Shuttle to an initial small disturbance i n the librational angles is presented. Figure 6-47 compares time histories of the generalized coordinates associated with the first four modes for the three different orientations of the Orbiter considered. The plate vibration was described by superposition of the first two modes of a cantilever and a free-free beam. The mode combinations are represented as Hij where i corresponds to the cantilever mode and j to the free-free mode. Thus H2 \ designates the generalized coordinate associated with the second cantilever and the first free-free mode term. Note, the plate is located in the plane of symmetry of the Orbiter. A s expected the system is unstable both in librational and vibrational modes i n the flight configurations (a) and (b), while the Lagrangian configuration and the associated panel vibrations continue to be stable (case c). Of course, the large amplitude response i n this case should be treated as qualitatively suggesting trends. Note that the vibration generalized coordinates Hij in Figure 6-47 are nondimensionalized with respect to a characteristic length (the system radius of gyration). Similarly, the deployment velocity is nondimensionalized using the characteristic length and time, which is taken to be the orbital period. Figures 6-48 and 6-49 attempt to study interactions between the librational dynamics, flexibility of the appendage and its deployment. The Orbiter is in the Lagrange configuration with a plate-type flexible appendage deployed along the orbit normal. The plate is oriented in two mutually perpendicular directions, out of and i n the plane of symmetry of the Orbiter. The results correspond to four different conditions of the plate in orbit. Both transient and post-deployment phases are considered. Initially the Orbiter is in the nominal equilibrium configuration corresponding to the undeployed appendage. Thus the applied disturbance corresponds to a change i n the initial equilibrium brought about by the deployment of the appendage. As before, the appendage vibration is represented by superposition of four modes corresponding to cantilever and free-free beams. A uniform deployment rate is used for a duration slightly less than 0.1 of an orbit. In Figure 213 Figure 6-47 Typical vibrational response of the flexible solar panel during and after the deployment for three different flight configurations of the Orbiter. 214 e=0 . E-0.1 , Plate Out of the Plane of Symmetry Initial Condition : Undeployed Equilibrium Configuration Figure 6-48 The Orbiter-appendage coupled dynamics during deployed and transient deployment conditions of the plate out of the plane of symmetry . 215 e=0 . E^0.1 , Plate in the Plane of Symmetry Initial Condition : Undeployed Equilibrium Configuration Deployed M-l '0 1 1 1 2 3 orbits Deploying , L = 0.075 m/s 1 1 4 1 5 1 '0 1 1 1 2 3 orbits i 4 5 Figure 6-49 The Orbiter-appendage coupled dynamics during deployed and transient deployment conditions of a plate in the plane of symmetry. 216 6-48 the plots show that the effect of deployment is to increase both the librational (yaw in this case) and vibrational motions. Note, the plane of the flexible plate is coincident with the yaw plane. The plate stiffness in this mode being large, it behaves essentially as a rigid body during the deployed configuration with very small values of generalized coordinates Hu, H12, H21, M22 associated with the plate vibrations. As before, during deployment the plate drifts because of the Coriolis effect. However, after deployment the plate is excited in its characteristic transverse mode, which is superimposed on the yaw motion. Note, for the solar panel oriented in the plane of symmetry of the Orbiter (Figure 6-49) all the librational degrees of freedom are excited upon termination of the deployment. The vibratory motion in this case is, relatively, larger. In the actual experiment involving the solar panel deployment, the Orbiter's librations will be controlled to a specified tolerance limit as indicated before (Figure 6-19). In the following, attention is focused on the plate vibration during a deployment maneuver in the presence of such controlled libration of the Orbiter. To start with the plate response to deployment or the controlled shuttle is considered separately as shown in Figure 6-50. As in the beam case, the vibrations are supressed during deployment with a drift in deflection, due to the Coriolis force. The plate vibration starts at termination of the deployment. In general, the deployment leads to the vibratory response with relatively large amplitudes. High frequency oscillations for the relatively rigid plate (nondimensionalized rigidity E = 400; with respect to plate mass, length and P orbital period) are as expected. Next, the plate response due to the controlled librations of the space shuttle is considered. All four modes of the plate are excited corresponding to the Orbiter's roll motion. A study of the combined effects of deployment and forces due to controlled libration of the Orbiter follows. Figures 6-51 to 6-53 represent vibratory response of the plate, deploying in or out of the plane of symmetry of the Orbiter in three flight configurations. Together they represent five different orientations of the deploying solar panel in space. At 217 e=0 . Out-of-Plane L= 0.075 rn/s , E =400 p '0.0 Deployed , E =26 p 0.2 0.4 0.6 orbits 0.8 orbit-*— p 1.0 '0.0 0.2 0.4 0.6 orbits 0.8 1.0 Figure 6-50 Vibratory response of a plate to deployment or specified motion during the Orbiter's flight in the nominal configuration. 218 Figure 6-51 Influence of the deployment rates on the vibratory response of a plate oriented tangent to the orbit. 219 Figure 6-52 Response of the flexible plate-type appendage as affected by its inplane and out of the plane of symmetry orientations during the Orbiter's flight in the nominal configuration. 220 Figure 6-53 Vibratory response of the deploying plate as affected by its orientation during the Orbiter's flight in the Lagrange configuration. 221 the outset it is easy to recognize that all the modes are stable. It is interesting to see how the Orbiter's controlled roll motion at around 13 cycles per orbit dominates the character of the appendage behaviour, which depicts the forced response occasionally modulated at its characteristic frequency. In particular, Figure 6-51 compares the plate response for two different durations of deploymnet (0.175 and 0.75 of an orbit). In both the cases, the response exhibits essentially the same characteristics. During deployment the vibrations are supressed but there is a drift in the deflection. Upon termination of the deployment all the four modes of the plate respond with a frequency corresponding to the Orbiter's roll motion. The similar trends are discernible in the other cases (Figures 6-52, 6-53). 7. C O N C L U D I N G 7.1 R E M A R K S Multibody Mathematical Model The introduction of computers in the mid-fifties renewed interest in classical dy- namics with the goal of efficient and effective development of mathematical models for complex dynamical systems. One of the outcome of these efforts is a specialized discipline of multibody dynamics and control, primarily motivated by developments in robotics and proposed large flexible space stuctures. In spite of the ongoing efforts at multibody formulation procedures during the past three decades, there is a considerable scope for lasting contributions to improve upon the present state of the art. The present contribution establishes a rather comprehensive approach to mathematical modelling of multibody dynamics applicable to a large class of earth as well as space based systems. Unlike the existing methodologies, the approach developed here leads to governing equations in a standard recognized format, which is well suited to stability analysis and control design. Importantly, it provides better physical appreciation of contributing forces in the equations of motion and their role in governing the dynamics. The generality of the approach in terms of system geometry, inertia distribution, choice of coordinates, discretization scheme, etc., is the central feature of this development. Characteristics of the versatile formulation have been listed earlier at length (p. 19-20). It is important to emphasize here that the procedure is applicable to a large class of diverse problems, of contemporary interest, likely to be encountered in the design of the next generation of highly flexible communications satellites, the Space Shuttle based experiments (Control Of Flight Structure-COFS, tethered satellite systems, etc.), the proposed Space Station, the Mobile Remote Manipulator System (MRMS) on board the Station, and many others. Indeed the thesis has attempted to tackle a difficult problem using an 222 223 elegant approach which leads to a versatile mathematical model. The following few salient features of the procedure help emphasize this point: • The model accounts for an arbitrary number of interconnected bodies forming branched and closed loop topology. In general, members of the system are flexible and deployable with interconnecting joints allowing six degrees of freedom. • This is the first time the classical Lagrangian procedure is applied to this class of problems. The key feature has been a development of the mass matrix associated with the quadratic form of the kinetic energy for such a multibody system. It leads to: (i) governing equations, in familiar and favourable form free of any holonomic constraints, are reducible to their minimum set by embedding the nonholonomic constraints if present; (ii) important system functions such as the Lagrangian, Hamiltonian, etc., which aid in the study of system stability and control. • The use of relative motion between members in the development of the kinematics permits separation of coordinates as specified and generalized during the dynamical formulation. This innovative approach has both a mathematical and physical basis, and vastly increases applicability of the model. • Using the concise matrix format, the formulation procedure approaches the problem making the model progressively more complex: starting with aflexiblebody followed by a system of rigid bodies with flexible terminal members and finally a cluster of flexible bodies. The procedure clearly identifies the system character during the stages of development thus imparting elegance to the mathematical model. • Since the multibody system under investigation is kept arbitrary, i.e., not identified in advance, the associated kinematics and dynamics are developed in terms of a matrix of geometry which, together with the constraint relations, completely define the system configuration. The kinematics can be in terms of coordinates which are relative, inertial or a mixture of the two, leading to the Lagrangian or the Newton-Lagrangian formulation. 224 • The formulation allows for alternatives in the choice of system coordinates thus permitting a varying degree of sophistication in modelling. For example : (i) flexibility is described using a general displacement field, which allows linear or nonlinear representation; (ii) flexible character can be modeled through discretization using lumped parameters, admissible functions, finite elements, or modal synthesis; (iii) rotational motion is described in general terms, i.e., a variety of representation such as Euler angles, Euler parameters, etc. can be used; (iv) the system reference point can be its center of mass or any specified point. Thus the mathematical model represents innovation at virtually every step. 7.2 Examples Illustrating Application of the Formulation Application of this versatile mathematical model is illustrated through a set of ex- amples of considerable current interest. They are based on a system consisting of a central rigid body with an arbitrary number of beam and plate type deployabale appendages. Equilibrium configurations of the system are established and their stability studied . It is found that for attitude stability the dynamic potential has a local minimum only when the flexible system is in the Lagrange configuration. Relations governing the equilibria and stability of the flexible appendages are also obtained in terms of the appendage attitude and flexibility parameters as well as the geometric stiffness due to orbital forces. The case of a cantilever beam in a circular orbit is studied in detail for three different orientations: (a) beam in the orbital plane; (b) beam in the plane defined by the orbit normal and local vertical; (c) beam in a plane defined by the orbit normal and the local horizontal. The results show that: • For all the three cases the transverse beam deflections are statically uncoupled. • A beam aligned along the local vertical, the local horizontal, or the orbit normal undergoes no static deflection. A beam oriented, approximately, at ±45°, or ±135°, (in the above mentioned planes) has the maximum deflection. 225 • By considering the maximum permissible deflection of the beam, an approximate first estimate of the beam flexural rigidity can be obtained as G = 3£ mf2 /6A , b 3 2 4 where 6 is the ratio of the allowed maximum beam tip deflection to the beam length. The governing equations of motion for the orbiting cluster were coded for numerical integration. The fortran program is applicable to a relatively wide class of gravity gradient spacecraft configurations. Its versatility is demonstrated through application to both rigid and flexible deploying bodies. The program, developed in a modular fashion, isolates the effect of flexibility, deployment, nonlinearity, etc., and was used to analyse spacecraft dynamics involving different levels of complexity. Four configurations involving the Space Shuttle and the Space Shuttle based experiments were considered. The system configurations analyzed are as follows: (i) The Space Shuttle idealized as a rigid body. (ii) The Space Shuttle based deployment of a beam type member. This attempts to simulate manufacture of structural members on board the Orbiter for construction of the proposed Space Station. (iii) The Space Shuttle based deployment of long (150 m each) twin beam type antennas for the proposed joint experiment by Canada and the U.S. called the Wave In Space Plasma (WISP). (iv) A 33 m long plate-type solar panel deployed (or deploying) from the Space Shuttle thus approximately simulating the Solar Array Flight Experiment (SAFE). Based on the results obtained, the following general observations can be made: (i) Linearized analysis when applied to this class of problems can lead to grossly misleading conclusions and hence should be used with caution. (ii) Stable equilibrium orientation of a gravity gradient flexible spacecraft in a circular orbit is defined by the Lagrange configuration with the minimum moment of inertia axis aligned with the local vertical and the maximum moment of inertia 226 axis along the orbit normal. The configuration, being relatively stable, is better suited to missions involving deployment of flexible members. This may be relevant to the Orbiter based manufacture of structural components for construction of the proposed Space Station. For the Lagrangian configuration, the Shuttle roll motion must be controlled within ±0.5°, to prevent instability. This is attributed to the coupled yaw-roll motions which are dependent on the associated moments of inertia ratio. (iii) For a beam type appendage located in the orbital plane, inplane disturbances lead to inplane response only. Similarly an out-of-plane excitation result in the outof-plane vibrations only. Thus, for inplane orientations of the beam, transverse vibrations are essentially uncoupled. On the other hand, for a beam located out of the orbital plane the transverse oscillations are coupled. The plane of oscillation precesses due to the Coriolis effect at a rate governed by the beam inclination and reaches a maximum when the beam is oriented along the orbit normal. (iv) The eccentricity of the orbit affects the system response in several ways. First, it causes a periodic variation in stiffness of the flexible appendage. This leads to the modulation in the beam frequancy reaching a maximum at apogee and a minimum at perigee. It also imposes a force on the flexible element (beam) which causes a cyclic variation in the inplane equilibrium position of the beam. The orbit eccentricity also induces a cyclic torque about the orbit normal which affects the pitch attitude motion of the whole system. This could lead to librational instability depending on the inertia parameters of the system and the magnitude of the orbital eccentricity. (v) The frequencies of beam and plate oscillations diminish during deployment attaining stationary values upon their full extension. In general, the deployment leads to a larger steady state amplitude of tip oscillations. The impulse due to deployment is an important parameter causing an inplane drift in the deflection of an elas- 227 tic element. Furthermore, it can induce a torque destabilizing the the librational motion. (vi) Response of the beam and plate-type members to forced excitations, caused by librational control of the Orbiter, follows the roll, yaw and pitch time histories closely with the inherent high frequency modulations superposed on them. Clearly, all natural frequencies of the structural members in the system should be far from the control frequency. The mode truncation study showed that the beam response in the first mode was always dominant for all the cases reported in the thesis. (vii) Some numerical experimentation was required to select an acceptable level of the computational error. This was dictated by the choice of the step-size. It was observed that inadequate specification of the permissible error can lead to misleading results. At times it failed to give a true picture of the system dynamics due to supression of the higher harmonics. 7.3 Recommendations for Future Study (a) The ultimate goal of this work is to develop a computer program for a general multibody system according to the mathematical model presented here. As pointed out before, such a program would prove applicable to a wide range of problems. Development of this general code should, therefore, form the highest priority. (b) The procedure used to analyse the static equilibrium and stability of a beam in orbit can also be applied to plate analysis. Cases where the Liapunov function is time-varying but bounded (eccentric orbit, deployment, etc.) are also logical areas of future extension. (c) More detailed analysis of the plate dynamics is recommended considering the important role of large solar panels on the proposed Space Station. (d) As the matrix of geometry can be used for physically connected as well as separated 228 systems, it could find application in problems such as the simulation of docking dynamics. (e) Although the equations of motion developed here account for geometric nonlinearities, the numerical computations were performed considering linear stiffness and displacement fields. The inclusion of nonlinearities in the computational procedure may provide additional insight into the system dynamics. (f) Similar to the deployment response, an analysis of the retrieval dynamics should be of considerable interest. (g) With a firm foundation for studing dynamics of a large class system established, the next step would be to focus attention on their control. (h) Dynamics and control of flexible structures in the presence of enviromental forces has received little attention and represents a major challenge for space engineers. BIBLIOGRAPHY [1] Williams, C.J.H., "Dynamic Modelling and Formulation Techniques For Non-Rigid Spacecraft," ESA Proceedings of a Symposium on Dynamics and Control of Non-Rigid Spacecraft, Frascati, Italy, May 1976, pp. 53-70. [2] Kane, T.R., and Levinson, D.A., " Formulation of Equations of Motion for Complex Spacecraft," J. of Guidance and Control, Vol.3, No. 2, March-April 1980, pp. 99-112. [3] Hollerbach, J.M., "A Recursive Lagrangian Formulation of Manipulator Dynamics and a Comparative Study of Dynamics Formulation Complexity," IEEE Transactions on Systems, Man, and Cybernetics, Vol. SMC-10, No. 11, Nov. 1980, pp. 730-736. [4] Ness, D.J., and Farrenkopf, R.L., "Inductive Methods for generating the Dynamical Equations of Motion for Multibodies Flexible System, Part I, Unified Approach," Synthesis of Vibrating System, ASME, New York, 1971, pp. 78-91. [5] Ho, J.Y.L., and Gluck, R., "Inductive Methods for Generating the Dynamic Equations of Motion for Multibodies Flexible system, Part 2, Perturbation Approach," Synthesis of Vibrating System, ASME, New York, 1971. [6] Levinson, D.A., "The Derivations of Equations of Motions of Multiple-Rigid-Body Systems Using Symbolic Manipulation," AIAA/AAS Astrodynamics Conference, August, 1976, Paper No. 76-816. [7] Broucke, R., and Garthwaite, K., "A Programming System for Analytical Series Expansions on a Computer," Celestial Mechanics, Vol. 1, 1969, pp. 271-284. [8] Rom, A., "Mechanized Algebraic Operations (MAO)," Celestial Mechanics, Vol. 1, 1970, pp. 301-319. [9] Jeffreys, W.H., "A Fortran-Based List Processor for Poisson Series," Celestial Mechanics, Vol. 2, 1970, pp. 474-480. [10] Bayazitoglu, Y.O., "Methods of Automated Analysis of Three Dimensional Mechanical Dynamical Systems with Application to Nonlinear Vehicle Dynamics," Ph.D. Dissertation, University of Michigan, 1972. [ll] Bayazitoglu, Y.O., and Chace M.A., "Dynamic Analysis of a Three Dimensional Vehicle Model Undergoing Large Deflections," Informational Automotive Engineering Congress and Exposition, Detroit, 1977, SAE Paper No. 770051. [12] Langrana, N.A., and Bartel D.L., "An Automated Method for Dynamic Analysis of Special Linkages for Biochemical Applications," Journal of Engineering for Industry, 229 230 ASME, May 1975, pp. 566-574. [13] Huston, R.L., and Passerello C.E., "On the Dynamics of a Human Body Model," Journal of Biomechanics, No. 4, 1971, pp. 95-102 [14] Sheth, P.N., and Uicker Jr. J.J., "IMP (Integrated Mechanical Programs), A Computer Aided Design Analysis System for Mechanisms and Linkages," ASME Paper No. 71VIBR-80, 1971. [15] Chadwich, C.H., "On the Dynamics of Systems of Rigid Bodies Connected in Open Chains by Spherical and (or) Revolute Joints," Ph.D. Dissertation, Stanford University, 1972. [16] Bodley, C , Devers, A., Park, A., and Frisch, H., "A Digital Computer Program for Dynamic Interaction Simulation of Controls and Structures (DISCOS)," NASA Technical Paper 1219, May 1978. [17] Frisch, H.P., "A Vector-Dyadic Development of the Equations of Motion for N-Coupled Flexible Bodies and Point Masses," NASA TN D-8047, August 1975. [18] Ho, J., and Herber, D.R., "Developments of Dynamics and Control Simulation of Large Flexible Systems," Journal of Guidance, Control, and Dynamics, Vol. 8, May-June 1985, pp. 374-383. [19] Singh, R.P., VanderVoot, R.J., and Likins, P.W., "Dynamics of Flexible Bodies in Tree Topology - A Computer Oriented Approach," AIAA paper 84-1024, 1984. [20] Wittenburg, J . , Dynamics of Systems of Rigid Bodies, B.G. Teubner Publisher, Stuttgart, 1977. [21] Modi V.J., "Attitude Dynamics of Satellites with Flexible Appendages - a Brief Review," Journal of Spacecraft and Rockets, Vol. 11, November 1974, pp. 743-751. [22] Bainum P.M., "A Review of Modelling Techniques for the Open and Closed Loop Dynamics of Large Space Systems," Presented at the 15th International Symposium on Space Technology and Science, Tokyo, Japan, May 1986, Paper No. c-6-1. [23] Kane, T.R., and Levinson, D.A., Dynamics: Theory and Applications, McGraw-Hill Inc., New York, 1985. [24] Hooker, W., and Margulies, "The Dynamical Attitude Equations for an n-body Satellite," Journal of Astronautical Sciences, Vol. 12, 1965, pp. 123-128. [25] Roberson, R., and Wittenburg, J., "A Dynamical Formalism for an Arbitrary Number of Interconnected Rigid Bodies with Reference to the Problem of Satellite Attitude Control," Proceedings of the Third International Congress of Automatic Control, Lon- don, 1966, Butterworth and Co., Ltd., London, 1967. 231 [26] Roberson, R., "A Form of the Translational Dynamical Equations for Relative Motion in Systems of Many Non-Rigid Bodies," Acta Mechanica, Vol. 14, 1972, pp. 297-308. [27] Hooker, W., "A Set of r Dynamical Attitude Equations for an Arbitrary n-Body Satellite Having r Rotational Degrees of Freedom ," AIAA Journal, Vol. 8, No. 7, July 1970, pp. 1205- 1207. [28] Ho, J., "Direct Path Method for Flexible Multibody Spacecraft Dynamics," Journal of Spacecraft and Rockets, Vol. 14, February 1977, pp. 102-110. [29] Hooker, W., "Equations of Motion for Interconnected Rigid and Elastic Bodies: A Derivation Independent of Angular Momentum," Celesstial Mechanics, Vol. 11, 1975, pp. 337-359. [30] Likins, P., "Dynamic Analysis of a System of Hinge-Connected Rigid Bodies with Nonrigid Appendages," International Journal of Solids and Structures, Vol. 9, 1973, pp. 1473-1487. [31] Hughes, P.C., Spacecraft Attitude Dynamics, John Wiley and Sons, New York, 1985, pp. 76-83. [32] Ho, J . , "The Direct Path Method for Deriving the Dynamic Equations of Motion of a Multibody Flexible Spacecraft with Topological Tree Configuration," AIAA Mechanics and Control of Flight Conference, Anaheim, California, 1974, Paper No. 74-786. [33] Hooker, W., "Equations of Attitude Motion of a Topological Tree of Bodies, the Terminal Members of Which May Be Flexible," Lockheed Missiles and Space Company, Palo Alto, California, Technical Rept. LMSC-D354938, Nov. 1973. [34] Russel, W.J., "On the Formulation of Equations of Rotational Motion for an N-Body Spacecraft," Airospace Corp., El Segundo, California, U.S.A., TR-0200(4133)-2, 1969. [35] Hughes, P.C., "Dynamics of a Chain of Flexible Bodies," Journal of the Astronautical Sciences, Vol. 27, 1979, pp. 359-380. [36] Jerkovsky, W., "The Structure of Multibody Dynamics Equations," Journal of Guidance and Control, Vol. 1, No. 3, May-June 1978, pp.173-182. [37] Keat, J . , "Dynamical Equations of Multibody Systems with Application to Space Structure Deployment,"Ph.D. Dissertation, MIT, 1983. [38] Boland, P., Samin, J., and Willems, P., "Stability Analysis of Interconnected Deformable Bodies in a Topological Tree," AIAA Journal, Vol. 12, 1974, pp. 1025-1030. [39] Boland, P., Samin, J., and Willems, P., "Stability Analysis of Interconnected Deformable Bodies in a Closed Loop Configuration," AIAA Journal, Vol. 13, 1975, pp. 864-867. 232 [40] Wittenburg, J., "Nonlinear Equations of Motion for Arbitrary Systems of Interconnected Rigid Bdies," Proceedings of the Symposium on Dynamics of Multibody Systems, IUTAM, Munich, Germany, August-Sept., 1977, Editor: K. Magnus, Springer-Verlag, 1978. [41] Wittenburg, J., "Dynamics of Multibody Systems," IUTAM (International Union of Theoretical and Applied Mechanics), 1980, pp. 199-207. [42] Ohkami, Y., Okamoto, O., Kida, T., Yamaguchi, I., and Matsumoto, K., "Dynamic Formulation and Simulation of Multi-Body Space Structures," 87th Congress of the International Astronautical Federation, Innsbruck, Austria, October, 1986, Paper No. IAF-86-238. [43] Ohkami, Y., Okamoto, O., kida, T., Yamaguchi, I., and Matsumoto, K., "A Unified Matrix Approach Applied to Dynamic Formulation of Complex Space Structures with Nonlinear Hinge Forces and Torques," 88-th Congress of the International Astronautical Federation, Brighton, U.K., October, 1987, Paper No. IAF-87-348. [44] Huston, R.L., and Passerello C.E., "Multibody Structural Dynamics Including Translation Between the Bodies," Computers and Structures, Vol. 12, 1980, pp. 713-720. [45] Vigneron, F.R., "Stability of Freely Spinning Satellite of Crossed-Dipole Configuration," CASI Transactions, Vol. 3, No. 1, March 1970, pp. 8-19. [46] Vigneron, F.R., "Comment on Mathematical Modeling of Spinning Elastic Bodies for Modal Analysis," AIAA Journal, Vol. 13, January 1975, pp. 126-127. [47] Cherchas, D.B., and Gossain, D.M.,"Dynamics of Flexible Solar Array During Deployment From a Spinning Spacecraft," Trans, of the Canadian Aeronautics and Space Institute, Vol. 7, No.l, March 1974, pp. 10-18. [48] Lips, K.W., and Modi, V.J., "Transient Attitude Dynamics of Satellite with Deploying Flexible Appendages," Acta Astronautica, Vol. 5, 1978, pp. 797-815. [49] Hughes, P.C., "Deployment Dynamics of the Communication Technology Satellite - A Progress Report," ESRO Proceedings of a Symposium on Non-rigid Vehicle Dynamics and Control, Frascati, Italy, May 1976. [50] Ibrahim, A . M . , and Misra, A.K., "Attitude Dynamics of a Satellite During Deployment of Large Plate-Type Structures," Journal of Guidance, Control and Dynamics, Vol. 5, September-October, 1982, pp. 442-447. [51] Novozhilov, V . V . , Foundations of the Nonlinear Theory of Elasticity, Graylock Press, Rochester, N.Y., 1953, pp. 1-47, 177-216. [52] Hurty, W.C., "Dynamic Analysis of Structural Systems Using Component Modes," 233 AIAA ./., Vol. 3, No. 4, April 1965, pp. 678-685. [53] Tabarrok, B., "Complementary Variational Principles in Elastodynamics," Computers and Structures, Vol. 19, No. 1-2, 1984, pp. 239-246. [54] William, C.W. Jr., and Earl, C.S., "A New Matrix Theorem and its Application for Establishing Independent Coordinates for Complex Dynamical Systems with Constraints,'' NASA TR R-S26, October 1969, pp. 1-13. [55] Meirovitch, L., "On the Effects of Higher-Order Inertia Integrals on the Attitude Stability of Earth Pointing Satellites," Journal of Astronautical Sciences, Vol. 15, No.l, 1968, pp. 14-18. [56] Kaza, K.R., and Kvaternik, R.G., "Nonlinear Flap-Lag-Axial Equations of a Rotating Beam," AIAA Journal, Vol. 15, No. 6, June 1977, pp.871-874. [57] Likins, P.W., Barbera, F.J., and Baddeley, V., "Mathematical Modeling of Spinning Elastic Bodies for Modal Analysis," AIAA Journal, Vol. 11, No. 9, Sept. 1973, pp. 1251-1258. [58] Rosenberg, R.M., Analytical Dynamics of Discrete Systems, Plenum Press, New York, 1977, pp. 211-217. [59] Hemami, H., and Weimer, F.C., "Modeling of Nonholonomic Dynamic Systems With Applications," Journal of Applied Mechanics, Vol. 48, March 1981, pp. 177-182. [60] Meirovitch, L., and Oz, H., "Modal-Space Control of Distributed Gyroscopic Systems," Journal of Guidance and Control, Vol.3, No.2, March-April 1980, pp. 140-150. [61] Meirovitch, L., and Baruh, H., "Optimal Control of Damped Flexible Gyroscopic Systems," Journal of Guidance and Control, Vol.4, No.2, March-April 1981, pp. 157-163. [62] Likins, P.W., and Roberson, R.E., "Uniqueness of Equilibrium Attitudes for Earth Pointing Satellites," Journal of Astronautical Sciences, Vol. 13, 1966, pp. 87-88. [63] Auelmann, R.R., "Regions of Libration for a Symmetrical Satellite," AIAA Journal, Vol. 1, No. 6, 1963, pp. 1445-1446. [64] Pringle, R. Jr., "Bounds on the Librations of a symmetrical Satellite," AIAA Journal, Vol. 2, No. 5, 1964, pp. 808-912. [65] Pringle, R. Jr., "Tumbling Motions of an Artificial Satellite," AIAA Journal, Vol. 3, 1965, pp. 1087-1093. [66] Pringle, R. Jr., "On the Stability of a body with connected Moving Parts," AIAA Journal, Vol. 4, 1966, 1395-1404. [67] Budica, R.J., and Tong, K . L . , "Shuttle On-Orbit Attitude Dynamics Simulation," AIAA/A AS Astrodynamics Conference, San Diego, California, August 1982, Paper 234 AIAA-82-1452. [68] Young, D., "Vibration of Rectangular Plates by the Ritz Method," Journal of Applied Mechanics, Vol. 17, No. 4, Dec. 1950, pp. 448-453. APPENDIX A : DERIVATION OFKINETIC A.l PARAMETERS Displacement Field The position and velocity vectors of a generic point in Bi relative to the local frame Li are given by Eqs.(2-1) and (2-2) as i f = p + i w, f = hrt + F ii. iT 1 Let the deformation w be expressed by the linear and quadratic terms in the deflection generalized coordinate as hence, from Eq.(2-3), piT = $ t T + 2[ ] . Admissible functions used in the discretization of flexibility are the mode shapes. Such mode shapes are, generally, functions of the form <& = ${pj() = $(/?), and hence \I> = \P(/3) , where I and p are, respectively, the appendage length and coordinate in the £ deployment direction. Therefore, from Eqs.(2.4) and (2.5), rf = di + (i - ^ ) { i i and P I' 235 + yt-}, 236 where: p = position vector to a generic point on Bi in nominal undeflected position ; 1 {p\} = projections of p* on Xi,yi,Zi , = 1,2,3 ; $* = matrix of admissible functions associated with the linear displacement term ; ty) = matrix of admissible functions associated with the second-degree nonlinear displacement field , / = 1,2,3 ; 8i = generalized coordinates ; ti = time rate of deployment ; dj = unit vector in the direction of deployment ; £t- dt- represents the deployment velocity vector . The procedure implies that the time-varying displacement can be suitably approximated by a finite expansion. Each term in the series consists of a product of a time-varying coefficient and a spatial function. The spatial functions are assumed given and describe the shape of the deformed body, while the time-varying coefficients are the amplitudes of the corresponding shape functions, and serve as generalized coordinates for the flexibility. While grouping the generalized coordinates associated with Bi in a vector Si of length n , the shape functions can be arranged in a matrix form $ l T s dimensions 3 x n l s and n l s and ty) (/ = 1,2,3), of x n\, respectively. Thus, the finite expansion associated with the elastic deflection is conveniently expressed in a matrix form rather than using the summation convention. To facilitate the subsequent development, it is desirable to develope alternative forms for $ l and ty). The rig x 3 matrix of shape functions $ l has the general form ^ = [{^'}y*]. j,fc = 1 , 2 , 3 , 3 = [{J2^M 3=1 = [K\ , (A.la) 237 where {«*}!* = { Whk) {0} {0} {/c'} * = { {4?} }; 2 J {0} { r ak \ ; {K*} 3 * = J {0} I {0} l{*''}8* and 4 = E4t> * = 1,2,3. Similarly, *!' = [A]' ], Aj = £ / c j fc fc f c y , / , * = 1,2,3.- (A.lb) 3=1 Using the above equations, expressions for several quantities of interest associated with a flexible body undergoing axial deployment are developed in the following sections. A.2 First Moment and the Center of Mass By definition the first moment about the body coordinates is / —i —i j £ drrii — i b , J n%i m = j {? + * = rrii{ c + H l i i T T 6i + 6j St} drrn , ' -T St + { 6< H} 6 } } { where c* = position of the center of mass of body Bi treated as rigid , = — I H i T t m I — [ = — drrii 5 $«' rfro ; r i i Jmi H)= i Jm, m ¥J drrii • 238 A.3 Linear Momentum By definition the linear momentum relative to the body coordinates is given by f fdmi=f {£\r' + [$*r = mi{lite + H i tAdrrn, + 2[£f i T i } . The terms on the r.h.s. represent contributions from the deployment and elastic vibration of the member, respectively. Therefore, m,- h and m - Hg define mass vector and mass e t matrix associated with the deployment and vibration, respectively. Here where : Note that H l c and H l e are due to convection and time variation of the shape function associated with the linear term, respectively. Also, H_ has similar contributions associated l el with the second-degree nonlinear term. The matrix associated with the vibration contribution in the linear momentum expression is H l T f = — i ir + 2[*f *j]|dm,-, Jmi l m = H [* i T + 2 [ 6f Hj } . The two matrices on the r.h.s. were defined before. 1 239 A.4 Moment of Momentum The moment of momentum in the body coordinates is given by f drm = j f x f + [ $l'T + 2 [ 6 i f x { i i ¥ = m {li pi + Pl ty) T dmi , ti} . T { Once again, the two terms represent contributions from deployment and vibration, respectively. Here, m,- p^ and m - Pg define the vector and matrix of mass moments associated t with deployment andflexibility,respectively, where P = — / m- J . e t T X ffi dm , { m = —[ { P + ^ i i T Si + Sj ty)-Ji } x { di+ ( 1 - } ) { $' P 6i + 6j ty'/ 6i } } dm , lT { i Jmi m = ~d d i + [^[ Pi - P l f - l H ] 6i + { 6i i T [ I [ P^ - T ] + P., ] }. The first term represents contribution due to deployment of a member treated as rigid while the other two terms, which are linear and quadratic in deflection, account forflexibilityof the member. Note : p*T = A . and Pl T = — mi f ^ Q'iT d fJ . A ~P * ' { i T m . ) dmi , m = [ { *jki mi ±rf J. A P) A"* }f 1 • m The nonlinear contribution in p^ involves terms like P^, P*j and P_lj. Their expressions can be developed as follows : fi PiJi = — m i f {ty t 6i} x { & i T 6 } dm , { Jm { = —ejki I m iT Jm, ^ At- ~£ k % T dmi , { 240 hence, U f mi Jm = — mi A* ~£ T k drrii ; * / {$ J. 6i} x {p\ & iT Si) drrn , i T m = —tjkl / Pi S { A•A \ T 6 dmi , { therefore , Pli = — eyjfei / i Pi A* A' dm T k { J m .• m Similarly, {S i Pii6i} T = — f {p x(l-f){6j Wit Jm,- <• = {^ ¥$» 6i} + { 6 i i } x 4,} T *t / (p) dm , { J (i - f ) ( y f * f y<) + ( y f * j * ) * ( * ) ) ^ } , hence " * t 7 m.- *-t L The moment of momentum term associated with the flexibility contains the Pg matrix, which can be expressed as ps l T = — f mi V m . = — / e [*<T + 2 [ « ? *j] dmi , L { />' + * ; i T Si + fi 9) St } X f* < T + 2[?f *}] dmi 5 The first term on the r.h.s. of the above equation is independent of deflection, while the other two terms are linear and quadratic in deflection. The contributing matrices PgS and 241 Pli can be expressed as u ~™~iLi~ $ pi , J drrii , p) A* f c = [ { ejki ^~ I dm )f 1 S i and —T Si Ph = — f { * mi Jmi = — / i T 6i}** drrn , i T eyfcz S [ A • A { T k ]j dm,- , i.e., Similarly : —T Si P}i = J- / ? [ ?f — 6* £ i ™ * ] / ey/c/ p) V drrii I k = — I \ 2 {<&tT y,- }x [ yf * j ] + {yf *• y,. }><drrii $ iT, m i Jmi L fl = 2( -J- / T )( )+( ( ^ * ^ * i ^ *5 * ^ ) } ejki mi Jm; A.5 ] dm,-, [ 2A l A * + V) A* (m) T m fc T drrii , dm,- Kinetic Energy /" £ • Z drm = [ { If rf + [ $ i T +2 [ tfj ] ] 5,- V { * if + [ $* + 2 [ sftfj] R = m< (£? < + 2 ^ « 4 j j f wU i+ ti ) , T ] ^ } dm,- , 242 mi WgS, and m, w\ where mi w , ee represent the mass, mass matrix, and mass vector s associated with deployment, flexibility, and their coupling, respectively. Now v\e = — [ Vi Vi dmi , iT = ±fm { * + ( i - f ) { *' si+ Jj *;»• i i } } r {li + i i - 1 ) dmi, {<&'iT ^ + * } } i + | d f [ffj-#:r The first term represents contribution of the deploying member treated as rigid, while the other two are linear and quadratic in deflection. Here : = £- f w i " cc = k. { " ce rrii 1 w { " ee w { rrii 1 drm ; j p\ <&'*' $ / pf 2 rrii lT Jmi rrii w $' $ ' Y, I d i(l) / , T dm, ; $ ' t T dm, ; [ il-Pt) J m *V dmi i The coupling term, wJ^5, has the form <6 T = — ™i [ Jmi rTi F dmi, T iT dm,dj WS T + Sj +{s; m r. Si The first term is independent of deflection, while the other two represent its linear and 243 quadratic contributions. The matrices W* , W\ , W\i and W} can be expressed as : s i W! cS WU mi — = i m 9 J m.- l IJmi T { . m { J m.- i m 2^ m Si y -[ mi ( l - ^ ) { ^ J m± | - / » c i S i } T T [ 6 j ^ ] d m i + *t ± - f mi Jm. "H Jmi dmi ; m mi j {s; dmi ; d i [ ? f *j ] dm , J. mi i T Pi *" f S*Wi = — s { l - P li i ) { S (i-^) i * T S i } T E ( s T '"•i Jm y y $ i T d m i , . - ) ( ^ * y ) ^ - + c i T t y Similarly, the term associated with the mass matrix can be expressed as ts = — mi = f Jmi ^-[ F F { [ $ t T 1 1 dm { , + [yf *y]] 2 T [* iT Hence : Jm, dmi ; + 2[S~J ¥}] dm,- , 244 Si EG, = — f m, y m i f [ L 6? »}.] + [ * j r Si } * j i T 1 dm,- , i.e., dm,- wi IL-Sl m,- 7 ^ / m m i £ ( i .• j Jm 4}i * T ) ( * ? m t ; 4m ) dm,- , i.e. ^ m = ^- f A.6 D 4 4 i ^ T M a s s M o m e n t of Inertia The mass moment of inertia is defined as in= [ = [f fu-ff }d , T T mi f Ci dmi. iT J m^ Since the pseudo-inertia is defined as [J ] = J % m . f* dm,-, it is apparent that the mo- ment of inertia can be obtained once an expression for the pseudo-inertia is known. Clearly, for a flexible member, the moment of inertia matrix should have contributions corespond- ing to its rigidity and flexibility. The interest here is to find the expression due to deflection only, since the matrix of inertia for a rigid member is obtained for a specific geometry either 245 analytically or experimentally. Now the pseudo-inertia can be written as = / { ? + * i T Si + Sj ty) St } { = rm [ [JJ] + [ { A j k +A. } k T (? + * i Si + Sj ty) Si } T drm , T + Jijk Si] + [ ST { r 3jk }6i}}. The first term on the r.h.s. is due to the member treated as rigid while the other two terms, which are linear and quadratic in deflection, can be further developed as follows : [ i 2}jk J r ^r\f 6i\ = = i-[/ Jmi i Jmi Si} T dmi] , iT {p)^, i m i* Si dmi], i.e., {A}jk = m Pj A** dmi } • Here, the Kronecker delta, Sjg, is equal to one if j = t and zero otherwise. On the other hand, Sjg_ is equal to zero if j = I and unity otherwise. Similarly, [ *? [ 4 + Ji jk jk i * i=^ / j { $ i r h } { $ + { If , r I i } T + f *!' St } p { iT * ' * } T ] dmi , = - £ - / [ {*f Sj-} {Al ?,-} + ? f [ pj- * l + pl ty) ) r mi jmi i.e., and Si ) dmi , 246 The expression for mass moment of inrtia can now be developed as [/'] = [ / f fu-ff }d , T T mi L As before, the rigid and flexible contributions can be readily recognized. Here : =/ [in r dm{ ? = moment of inertia of body, treated as rigid, about local coordinates ; 3 {P } 2 = jk m i { 6 3 £$>,,„{ jk J2 +4 l lm }}~mi{ m l Ji. k + J* . } ; k 1=1 m = l 3 [lihk = rm [ 6 The jk 3 Y £ 1=1 m=l *im 4 m ] - mi [4 } . k kinetic parameters such as momentum, moment of momentum, etc., for a flexible deploying body have been obtained before (Eq.3.4) as indicated below : t = c* + H i T 6i + { lj #!Si } ; K = *i + j;[ HI -Hi\ Si T HI = H + 2 [ H} l p* = Ci di + r TjiT 8 Si } ; T jAH-Pi _ p iT , — SS + sUP Sl t r + { Sj Hi, Si } ; } ~di H T iT ]s { Si [ + Pil}]+\Sj P} m P i - P., } + P }Si}; T i+ el Si}', [p[wi -2wi +wi ]}+w Si ; i c i weS T d; Hg + T si j-{wi -wi ]+m s s l e + {sfwj Si } ; T 247 Wig =Wl + [Wi 6 ] + [ sf Wj T l i St } ; m I = Z + [ {J5>L *••] + [ * f [h)im ^ ] • 1 They are further developed in terms offlexibilityintegrals here. Terms with underbars in the above equations are those arising from the quadratic term in displacement. The terms on the r.h.s. of the equations are defined by Eqs.(A.2) and (A.3) below : c* = position of mass center of body Bi treated as rigid ; H= f — { m y I * m #" drrn = [ { E Jmi ' = ^r H {A-2b) H ]; L i Hl = ^m $ drrn = [ { E (A.2a) [ i * p J m i $ , d t i m = [{ e <Cjk }k ] ; (A2c) y m = I{ £ >* 1 ; ( A 2 d ) Jm{ {tS iU i £ + liS ) it jt }, ] ; (A.2e) m ^5 = [ { *y*i — m / i P i Py A*** rfm - } f ] , t Jmi = [ { ey« (tS je + tiSit) E (A.2f) }F ] ; m p l * = £ e * " / ^ m P T £ 66 H = [( = [ {ey H i = I ' £ W E i E T = [ EJKL ^m* (A.2g) n T m Pj K 1; E E .^m* 1; i I m m m «' ^ > / ^ ^ ^ = d t m P ^ n drrn }f } , Jm{ (rs + ttSjt) E jt >F 1; (^20 i Jm.- = & f m w << m,' ~ w J = £ m „, $ J * " ' t = ^r[ < = nT wi m t ^ p i m i t f m "' $ » • = [ £ j iT / / /4 t m i £ i; ^ S * k ( 1; i = EEE [ d m i y i; *' ** 2 j = [£ r Jm.- £ „ jk ) T fc rm { 6 A M (^2») *" *' ^ = [ E E E n 2 ) (^M ( A 2 2 1: E • 15 i (^- P) i = moment of inertia of body i as rigid about local coordinate ; {P hk = A 2 f c ( 1; EE E fc - £ * » * = [ £ £ £ T $ , i$ i T ' = ^r I w » <*'* * " ' & = ~ m * i $ « w T - m / { , {A.2q) n J2 £ Sim { 4 l m + 4 }}-rm{ m l 4 j k + 4 k j }; {A.2r) i=0 m=0 n n £ £ *Zm [ /=0 m=0 [/ ]jfc = 3 {Ja}y* = — [ m t / » dlilm ] ~ i[ m 1 = ( t h + Pi 4 Jm- [ sj. s i 3jk J £ Jmi I4hk = ±-f m + + ±%jk }; 1 ! (^- ) 2b . {AM) i T ] dm,- = [ £ £ i Ft m 1; (^.2«) °) 249 and /£)'=—/ drrn = }Ai; (A.3a) i Jm m { Sii = ^-f r (i - i Jmi m = — eyw i = e ^ / ( A 3 6 ) [ * ( * ) * y + P) (1 - T ) 1 d m * ' U,- + B{ - i f B J 1 ] ; J fc ; (A.3d) J mi kl I2 * S J-T + * y Si(m)]dm,- = 2 ? A j + ?Aj-; r a r (A.3e) m; *(0 / E i i m ( "*t Jmi c *J t • = 2 E - 7^ ) c / 2?A? - 1 * i * ( 0 I \Bi - j = 2 E *W i i A « iBf ] ; (A.3/) *t l ! 'U-M ?UJ< ] ; + ^ ' - I (A.3*) * E l S y W * J + Ay. A y T ] drm = 2 E ( I W>i=^ 1 Jmi S= 2 E — = Ef i dmi = 2 ?A p) f jU € = ^ (A.3c) 2 = ™- > f m B f - e ^ " m.-i m Bl *t [ jkl i™ =< ~ T** J mi m P d m i y fA} 1 + ? A J , r ] ; (A.3*) y = E I ^ / £ } i AJ-™t ^ 4 ] = My? ; U-3i) i Jm m y { [ 2 4 ]y* = ^ m t / [ Py *i + Pi * y 1 ^ Jm-i = [ My + Mi 1 • (^-3*0 Eqs. (A.2) and (A.3) are denned in terms of body integrals involving elastic deflec- 250 tions in spatial coordinates. These integrals are to be calculated once only for input to the computer program. The quantities defined below represent all the necessary integrals for calculation of the kinetic energy in case of flexible deployable body using a second-degree nonlinear displacement field. Integrals Associated with a Linear Displacement M o d e l {ifjk} = — i m {iCjk} = {«t}y* dmi ; I i K'i Jmi — m / i }jk dmt = Jmi f ( jt + &3I JT ) — c = — iEjf] / i m {Ki}jk K}mn p) 8 mi iF$?\ {AAb ; 1 {i Vinn} (AAa dmi J j . (AAd m {Ki}mn (AAc dm{ ; (AAe Jm{ = A . f { }i Jmi {^ = — {^mn Ki k K mn d m i ; m iC?k n\ m I {<}ik Jm{ i = ^rf i m iDTk\ {i9mn) £ d m i (AAg ; Jm{ Pi Wi)jk I Pi Jm{ Wi}jk I dmi; (AAh dmi ; (A At m j . t iNTk n\ (AAf dmi ; ^~ i m m 77T- {«(•}£„ TTT / P) i i}mn K dm{ ; {AAj (AAk m; [t*6jt + ZiOjl) Jm{ Note that eight of the eleven integrals defined above would be absent if their is no deploy- 251 ment. Four of the eight intgerals due to deployment originate from the convective term. They are (AAb), (AAe), (A.4/), and (AAt). The other four integrals contributed by deployment are (A.4c), (AAg), (AAh), and (AAk). The latter set appears in the dynamics as the shape functions are time dependent. Of course, in absence of deployment the elastic body may vibrate, translate and rotate. The three integrals (AAa), (AAd), and (AAj) contribue to the dynamics in this case. Integrals Arising from the Quadratic Term (A.5a M y 1 - mi J„ * j dmi ; U) -/ \= n mi J M ? ]= (A.56 p\ 9) dmi 5 i r I A* AJ^ dm ; T { mi J 1 M j r i - mi ./ A}, A y r m ^ It p£ Ay A ^ M ? ]= m{ J T dmi ; (A.5d dm,- ; (A.5e m M f ] == i / *j. Al(/) dm ; mi J f { \Bj ] = (A.5/ (A.5 *y dm, ; l mi J ff n IB* } = i / mi J,"i IB* ] = mi J, JBf (A.5c i *y ' dm,- ; (A.5h pjpjtfjt*dm,- ; {A.5i 1 ] == i f *y 'A!(/) dm -; mi J m • 1 t {A.5j 252 [f?BH ] = — / p\ dm,- . (A.5A) Evaluation of the kinetic energy can now be summarized as follows : • Choose the displacement model, i.e., 7cj.y, A^. and Aj.y are known (Eq.A.l). • Evaluate the integrals due to body flexibility ( Eqs. A.4 and A.5 ) analytically or numerically. • Calculate components of the characteristic kinetic quantities of a single body as defined by Eqs.(A.2) and (A.3). • Use Eqs.(3.4) to calculate the characteristic kinetic quantities. APPENDIX B : RELATIONS INNONLINEAR ELASTICITY As mentioned before, we are considering problems with geometric nonlinearities only. This is the case of a deformed thin beam or a thin shell with large rotations as compared to the strains which are considerted small with respect to unity. The expressions for the nonlinear strain are as follows [51]: en 6 1 1 1 + 2 + (^ 12 + Uez) 2 + {^eis - W ) 2 22 + (^ 12 - We3) 33 + (^ 13 + W ll e e 622 e2 2 + \ ^33 ess + - ^23 e3 + e 2(^e ei3 eis + n(^ i3 + w ) + e (ie e + en(^e e e 2 e 2 + (^e + Wei) 2 23 + (^e - W i) 2 23 1 2 C 1 2 i 2 e 2 2 3 ) e - w i) e e 2 2 e (^e + 3 3 e2 1 2 3 3 1 3 -w ) + e (ie e3 2 2 + w i) + ( ^ e i - w ) ( ^ e i + w ) ; 2 3 e e e e 2 e 3 - w ) + (^e + w ) ( ^ e e2 i2 e 3 + w ) + (^ei - w 1 2 e3 3 e 2 e2 3 -u ) ; 2 3 )(ie2 e l 3 +w ). el Comparing the above equations with Eq.(3.33), {£} = { e} + { e T Ei e } + { e T E{ u e } + { ^ E 3 3 u e 3 = 1,2,...,6 } , gives 0 0 0 0 0 0 0 0 0 \E\\ = \ 0 0 0 0 \ 0 0 0 0 0 0 0 0 1 0 0 0 0 £ 0 0 0 0 0 x O O 4 0 0 0 0 0 0 0 0 0 i I 0 0 0 0 0 0 0 is 0 0 0 0 0 0 0 0 0 0 1 0 0 0 4 0 Ls 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 i 0 is J 4 0 0 [25?] 0 0 0 0 0 0 0 0 0 S Also , "0 0 0 0 0 .0 0 0 0 0 -1 0 0" 0 0 0 0 1. "0 0 0 1 0 .0 ; \El\ = \ 0 0 0 0 0 0 253 0 " 0 0 0 0 -1. ' 0 0 0 -1 0 . 0 ; \El\ = \ 0 0 0 0 1 0 0" 0 0 0 0 0. 0 | i 0 \ 0 0 \ 0 0 0 0J 254 •0 -1 1 0 0 .0 0 0 0 0 0 o • 0 0 0 1 2 0 . 1 0 -1 0 0 0 \Et\ 1 2 •o 0 0 0 0 0 0 1 2 0 0 \El\ -r l 0 0 0 o . 1 2 .0 and 0 \E$\ S 0 0 1 0 1 0 0 0 0 0 0 1 \El] o ° \ s o 0 S I 0 0 \El\ 0 0 1 0 0 [Et] 0 0 The linear strains and the elastic rotations are expressed in the form { we } = L „ u; , { e} = L w , e where w is the displacement vector. L and L e d dai 1 Hx 1 HXH2 2 d<xx 1 HXH3 aH a«! 1 H2 H1H3 da3 d da2 H2H3 9H9. dac 9H3 da2 1 H2H3 3 J__d_\ ~ H^ at I 3 dH da2 H1H2 9H Hi (d(X/Hl) H y < <*« *3 , in their expanded forms, are defined as u H d(l/H ) H V ax 0 3( , 1 cM H a3 ' H d(l/H ) H v . I d ) Hr, i2 aa, x •I 2 3 glf8('/fll) H y 2 « | _1 ~T H a2) 1 2 o 2( 2 2 1 a HH 2 3 3 ~ H x H ,d(l/H ) H \ a3 2 d da3 1 H3 i fdHn V a3 H ,d(l/H ) H \ a2 3 3 J_d_\ I a i i Q\ ^ 3 a1 3 2 n i / 9H3 TT J M H2H3\ a2 - 3 ; 1 l dH3 i II d\ H~TH^\~o7 ~^' 1 1 3 a~[) + n a HXH3\ a3 ,3H 1 l° Hl V~H-[H^\ a2 + Hi-?-) i IT d\ + n l a ^ ) 0 dHr, 1 (°H 2 | HjW^\ai + 1 as 0 where if, "da. 'dctj da.j 3 = 1,2,3. 2 255 Here : Hj = Lame coefficients ; ctj = curvilinear coordinates . For cartesian coordinates, Hi = Hi = Hz = 1 , and the expression for L and L w reduce e to : APPENDIX C :DIFFERENTIATION OFTHE MASS Cl MATRIX Partial Differentiation of the Mass Matrix From Eq.(3.29), the syst em mass matrix has the form M = T j [ Mi + A B A] T , T v where Mi, A, B and T are as defined in Eqs.(3.27) and (3.28). With j and k ranging v over all the system coordinates, the differentiation of the mass matrix with respect to the state parameter q can be written as r M, qr = 2 T„T , qr [ Mi + A T B A] T + i f [ 2 A , T v qr BA+A B, T A]T , qr V where : '0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 W0,gr 0 r>qr E 0 E0>qr 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 .0 0 0 0 0 S]k{Bk,q Gk + BkGk,q J wL + m ]>qr "ei •3 T m,-•3 weS r r T ' \ m* ie>qr W 3 V 0 0 hisi w m3- 9i ,qr Ci+ kT • 3 9i kT C 3k,qr m 0 3'9r G m m °3k 3>qr 3 m W 3SS + WjS,qr hi si m 3 G S S k 256 k + C . >?r k + C my G 3s kT C 3k,qr 8}i$klH 3k,qr .Sym. T 257 A 1 — 0 lk L '<7r The partial differentiation of the B matrix with respect to a system coordinate is n3 -6, kO r c+ 3 k . ™ikC k ,qr °kl + m £*m/>«r k-1 G —m-V'«'- fc-i + ti.j'"J fc-lHr w m i k i q k L rnjkC 3k,qr m G + —SnSki r 3k c-.y-i r- & ciz\, k + m G i k'Vr m G ii ' i k-l'Qr 6*k m ,gr ClZ\+ C i>1r l 6 ik • -I qr m 3 ik,qr C k + m ikCk ,qr m —SnSkoSki i ^ r m G i k>9r C l + • SnSkoSioSjiSkiCl,qr . sym. Here <5ycf = 6j ( l + ^ko )> and <5y£ = 6fco<5yi6fci • The underbarred values are the only terms 0 in the B matrix to be considered when the system reference point is not its mass center. 258 C.2 C o m p o n e n t s o f the M a s s M a t r i x rrij = mass of Bj rrijk = rrij rrik n m? — X] Sj rrii t=i fc fc rrij = rrij m m? = m? m k k n ™ jk = £ ji s s m i i=i L k j k = JZ i m si 3 s i ~ °i k °i bi i=l n t'=l Hf k = SJ Si C{ [ Ii - mi b* b iT + m fb< * 6 « T ] C\ { \u\ = [ii] + [{/.}jfc *,•] + [ yf [m ^" Jk [JJ] = moment of inertia of body B{, treated as rigid, about the local coordinates n {4hk n = mi{ 6 £ £ *, {4 jk ro lm + 4 }} ml 1=0 m=0 n [4\jk = rrn{ Jk £ 6 n £ Si 1=0 m=0 m Jl \ lm -\J 3 ) X jk - {4 jk + 4 .} k 259 Ke = 1 + f < ii K -^'] = dJ W +±-6j T s = Pj + Rj Uk + w S} tf = S}*ai + k + j ^ f T [ - 2W' + Wi ] Si cl [Wi-Wi] T s &k [ ^ k ; Cf?i V = <? + Hg Si T K = di+j.[Hi-Hi] Si T pi = cidi+[-[P -pi] -diHi" i T c s <+ —T = P 68 + &i Ph l T m ^ = [ { e y „ , ( r ^ + tiSjt) E m r - {rtn + Wi) [ {*iH {J&y* = + t a £ m l k >r ] igL }f 1 } £ = j? + uJ*' Si T FT & + 2[ s) ty\ } dp{d(3^ dp 3 >d0> dtjJok j:{ 6t l ii- ei]6i} T p pt 260 \dvJ- " 1 Bk = U + Sjk — r L dp J O 3 f c Note that the following quantities are constants as defined by Eq.(3-4) and (A.2) : Hi, Hi, Hi, p^, pi„ PI, p* , wi, wi, wi, wu, wu, wu, [j*] . f C.3 Partial Differentiation w.r.t. Mass mjk,mr m\ = 8jr rrik + $kr rrij = Sj r mr m k,mr = Sj m + Sk rrij m r k r , m = Sj m + Sk rrij k rn jk,mr = S rj S rk L kj k irn QT o r /-ik kZr /~tr — &j &k W r ° °Jfc r H3 ,mr — Hi k,mr=s rj Sj Sk H} ,m r s rk ct [ i , r Titty -i r i rT+fb r mr {{I^jk * - ] ' + [ t [IrUr = [IlUr +^~[ C.4 jk ii rT] ci [P3}jk 6 ] r L Partial Differentiation w.r.t. the Deployment Coordinate n £fc '*r = _C i t=l J m S} S lk C k b l,ir k C lk 261 H jk,tr=irs)siHi\ ir i=i Hi , = S} Si C\ [ 6 [ I , k lr ir +m fib** i fb r + iT £r -m r fb , iT r [IrUr^[I riUr+[{I 2}Jk,tr ir [b , b r ] er rT +b r b , rT &r ci **] n n 1=0 m=0 2 Sj a,i,L — Sj E l l=i k s r ak,i = Sjr n , r S j s X Cf_ l ai, ir j = k— I k j k '*r = j j Utr +Sri C? b' ,£ b S b ,e = d r k a r K,,=—-\H:-H:\ t* s t 1 Pi * = r c 1d r e - ^ [P 'c ] 6 + ± { P ' - P- } T r c T >€f . 6 r 262 prT pr Ur = S Pc l T it T [ { e jkl £ iv j mk }l } m = iL [ { Jkl £ l \tr P 7 6 e r d 3 mk }f m Pw.«r = j t I { ifcJ S e _ ! rglnk m {J^jkUr kUr F = Sji { £ - r gjk } tjr d ^ 2 "kUr C.5 )T ] d$k dw 3 Partial Differentiation w.r.t. the Translational Coordinate t'=i t=i Sj- a,i k , s u — SJ 5* 1 ar 263 r - 1 Sj i C.6 ci j b k— j i S a Partial Differentiation w.r.t. the Deflection Coordinate t=i ,8? = )<k >6? K + 6>6° = farf* H\ Gt ri S +6 T y«r { K [ fa Hg , s r ri »«• +Pe +P , s T r } ] rT 6 6 6 n H j ,8? = T , J k l > v i=l k s Hi ,6? k H s k = SJ Si C\ [ 6 ri [I -m r V b rT ] r t S . +m [ \V fb< T ] , . ] C\ + i Sj Si C{ , s [ It - rm V b iT + rm *6* a s [Iris; = [ ( J s Ii-miVb^ 3 r 5 ; Jo,-,*. = SJ' £ k,J r = 6 Clss k r n rlT [l£\ + [J£] jk 1 5/ 1=1 a J * element in {I^jk = f[ ,6 } k7i 6] r { i s ( « ) }jk = s row in - { <e + ) j k ] + 2[{I (s)}j )jk= ( 2i ) mipW jli sjsicl s w ,j k r r+ p d C?.!a [ ^crc - 2 ^ + ^ . +S 1 M 1 S Ct^s r k e ] *r a } t 264 Sj k ".s? = Sj )ai s +6 b iS }= —rT Pe >$r s b ; s Y \H -H:} r r [j-rl ]= C * b' , s +C*, s ri c c- P p e] + H d ]+l { r 6 r i [ P*i + Pcl c T - Pll - PllT Ps , s=[{P S^)}f} T l S {Ph(s)}7= -rowofP^ j jfi* — W - l ° r dw 2 -Bfc,gs — Sjk C.7 .8* <^j j d8?dp 3 Partial Differentiation w.r.t. the Rotation Parameter j ,6r — W - l ° r ><?r re rrt c * rik \ -"fc ^8 — 2-^/ t'=l r { 3 0?>*r = { S'^t'^r ^ Gfi ,6r 3 , 6 r } e [ fartr H\ = i k T ] re H ,0 3k r = Y !i i i >°r S i=l S H k kli * , ^fc °fc + ° t i " * felt rii , /'-(fc fcli s~ti °fc + ° t t °fc><?r 0 } 6r SJ Si \c{,e [ii-rm r + C\ [ ^ - rm [V b iT +rni Cj [ n j^2 si S j s _1 c fb iT t-i*er 1=1 Sj jGiidr [6* b iT + Gi,gr b . *l + [V *?T]]CJ + ># fb iT] } Cle + ># fr*' ,gr ] &k APPENDIX D : CHARACTERISTICS D.l OFB E A M MODES Fixed-Free Beam The characteristic functions of a fixed-free beam is given by ip {z) = [cosh(^) - cos(^i)] - o- [sinh(^) - s i n ( ^ ) ] , n n where £ is the length of the beam, X are the roots of the transcendental equations n 1 + cosh A cos A = 0 , and an are constants dependent on A . The values of A and on are shown in Table (D-l) n n below for the first four modes [68]. Table D - l The first four eigenvalues and associated constants for a fixed-free beam n A n 1 0.7340955 1.87510410 2 1.01846644 4.69409110 3 0.99922450 7.85475740 4 1.00003355 10.99554070 The modes, ipn, of a fixed-free beam satisfy the following relations : • boundary conditions V n = ^ • orthogonality conditions J Q = 0 at z = 0; and V'nV'a dz = 6 ns I; • ij)n and their derivatives satisfy the relation dd^4 = 266 = ipn ; = 0 at z = t ; 267 • with primes representing differentiation with respect to z VnW=2(-l)'1-1 ; = (-i) «» Y ; n+1 ^'(0)=2 ( ^ ) <(0)=-2« 2 ; (^) B 3 . Denoting cr /A by /? , using the above relations and integrating by parts gives : n n re j J r/>„ dz = 2 { 0 / ^ 1 D.2 iP'ndz * } ; = 2{ dz = J^Tpn^m ^ n 2 K } ; 8mn ; dz = 2 { ( - I ) - * --!.}; Free-Free B e a m The characteristic functions of a free-free beam are given by : <Pi{x) - 1 ; <p (x) = V 3 ( l - £ E ) ; 2 iPin(i) = [cosh(^p) + cos(^^)] - <5 [sinh(^^) + sin(^^)] . m m> 3. 268 Here d represent the length of the beam and /z m are the roots of the transcendental equation 1 — cosh fj, cos fi — 0 . As before, 8 are constants dependent on /^ m . The values of /z m and 8 are given in Table m m D-2 below for the first four modes [68]. Table D-2 The first four eigenvalues and associated constants for a free-free beam m Mm 1 0 2 0 3 0.98250222 4.7300408 4 1.00077731 7.8532046 The mode shapes, <p , of a free-free beam satisfy the following relations m • boundary conditions J2 " = m rf 2 • orthogonality conditions • ip m J3 ^3 = 0 , at x = 0, d ; <Pm<Ps dx = 8 d ; ma 4 .4 and their derivatives satisfy the relation = ^-<p , m > 3 ; m • with primes representing differentiation with respect to x : / <Pm = d 8ml d2 d Jo Z £>m = •/ 0 <Pm <Pn /' Jo /' Jo (Pm <Pn — ; \/3d2 e 0ml dx = d 8 Tl dx = ll m 8 mn r Om2 Some of the useful integrals, evaluated numerically, are presented in Table D Table D-3 Some numerically evaluated integrals - n: 4 2 3 0.85824 -11.74322 27.45315 -37.39025 1.87385 -13.29425 -9.04222 30.40119 1.56451 3.22933 -45.90423 -8.33537 1.08737 5.54065 4.25360 -98.91821 f 4,'m 4'n dz 4.64778 -7.37974 3.94150 -6.59312 -7.37974 32.41667 -22.35200 13.58119 3.94150 -22.35200 77.29883 -35.64684 -6.59321 13.58119 -35.64684 142.89893 2 f z <p'm Pn dz 0.00000 0.00000 18.58910 0.00000 0.00000 0.00000 0.00000 40.59448 0.00000 0.00000 -12.30262 0.00000 0.00000 0.00000 0.00000 46.05002 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 13.85641 0.00000 0.00000 49.48082 0.00000 0.00000 13.85641 0.00000 8.92453 2.00000 -4.75938 3.78433 -4.11964 0.75948 1.99990 -6.22209 3.38317 0.21566 2.22219 1.99998 -8.16832 0.11983 0.61653 4.16851 1.99974 1.50000 -4.10448 4.01302 -4.00323 0.10453 1.49990 -4.58469 4.13722 -0.01303 0.58479 1.49999 -5.40835 0.00340 -0.13750 1.40854 1.49981 3.07690 -6.95745 5.01359 -5.72001 -6.95745 23.76962 -24.24208 17.22463 5.01359 -24.24208 52.34668 -43.98531 -5.72001 17.22463 -43.98531 91.43974 2.26110 -6.00807 5.52626 -5.50034 -6.00807 19.46038 -22.69086 19.40503 5.52626 -22.69086 41.57977 -42.19564 -5.50034 19.40503 -42.19564 70.78823 x f < I 1; l £ dz e 2! 1.3 i 4 ; t 2 e f Vm K e dz \f z^ i>' dz t m n r APPENDIX E : INTEGRALS ASSOCIATED FLEXIBILITY The WITH integrals required due to flexibility of the appendages, corresponding to Eq.(A.4), are listed here. Note that for a beam j, k = 1,2 and for a plate j, k = 2 : = — mi = —i 1 J m = drrii 1 f Jm. { }yy d m i ; mi ( <W + L« 77 ) mi £ + iiSmi) Jm, / 1 mi (£*6mi J + tiSmi) = — / { t}yy t Jmre Pm m Wi\j3 Pm { Ki}jj Pz drrii ; t m 1 i — — f — m i = — m i Wi}jj i Ki}kk drrii ; Jmi / Pz x i Ki}j3 { K'i\kk d m i ; Jmi = ~T / mt i Ki}kk dm, ; Jmi mi = i *}33 K jm. I Jm{ P 3 {«<}yy i Ki^33 270 i ; {«5}yy Pm m — ~ d m t Jm t {«t}fcfc i; dm t / 1 mi (£*6mt / rfm »- i K'i}kk ; i dm ; rfm * 271 O n using the data in Appendix D, the values of the above integrals are as follows: BEAM /fx = 2 ( /?„ 0 ) ; cf x = 2 ( 3? T 1 / M= 2(0 0 ) ; / 2 2 = 2 { « 5 m l f3n } ; n 0); = 2(1/A» PLATE 0), 1 c? 2 = 2 ( 0 (-1)- ^ 1/An); 2 = 2(0 *yy ?y = • ^yy = «yy 7yy ; = F 0 _ ^ 1 c 2 2 = 2 {6ml (-l)"" } ; ); 2 dyy = dyy = 0 , dyy = cyy - ffyy . = 0,f v ™=^ ^22 = ^ 22 d = 0 22 = 2 {5 m l /A n } ; ~ Tr^K-i)"- }, ^ 2 2 = ^ 2 2 - 722 ; 1 ; ^22 _ 022 > 2 ^22 = 0 > d 2 2 = c 2 2 - fff2 . Here the subscript n is associated with the n-th mode of a fixed-free beam. The subscript m is associated with the m-th mode of a free-free beam. 272 BEAM 77>22 ^11 U 0 0 0 0 0 u 7?22 *22 *22 0 T 0 5 0 1^22 _ ^11 - r ll °11 u — u vA 0 75 0 0 0 0 0 «^mn 0 0 22 ^11 D u 0 U 0 0 0 0 0 J /o22 °22 0 0 [U] — l^sr J \ mn 0 Jl. C|| — [<5sr J ] mn r ll 0 V 22 o22 11 VV22 0 D 22 6 1 r>22 _ mn] ^ 2 2 — [<*sr -7, 6 D 22 11 0 7 0 0 0 0 0 p22 -^22 0 JT-22 -^22 22 0 J Nil 0 °22 0 T 0 nil -^ll 0 -&22 7r22 0 m 0 0 72 8 1 0 ^22 0 0 PLATE J mn 0 AT 2 2 0 0 0 Ji. mi Where the subscripts m and n are associated with the ra-th and n.-th modes of a fixed-free beam. The subscripts r and 5 are associated with the r-th and s-th modes of a free-free beam. APPENDIX F : BASIC KINETIC QUANTITIES The basic Kinetic quantities arising from the flexibility of the body, as defined by Eq.(A.2), are : PLATE BEAM [7u 722 o] m=[ tC22 0 ] H* = [0 Hi = [ivtl i^22 0] Ht = [0 Pc = l pi- r r pi ^i t ^ l l o] PI = [ -U if pi = 2 K = [ i ll G wi =-" - ^l2!2 1 iD\l ] o] S3 Pt = [ iFg - iF 22 ] 2 + = P e2 = Pis = [0] [o] l 0 Pk = P S2 = PL wL = [ l iC\l} wL = [ ^ ] wi t e 2 2 2 2 ] = [ iV% ] wis = [ iAl ] wis = [ ^11 ] wis == [ iPR + iFg ] wi = [ 6 273 d Pil + iDll ] e5 0 I w == [ iD\\ i d = P C2 = pis = II ] T jrll _i_ J?22 1 i + t^22 J 0 pii Ph = [-tii922 + tC 2 ] i = [ iVii 22 i^22 ~ii 22 pis = \ »'-^22 " = w == [ ii + " ce s Pis = [ t ^ [ ~l% 1022 ^» : " cc 77-11 ^ a = [0], _ i 0] [ — A ' t'^22 pi — W iUf 2 Pi = [ — pi 0] o] _ el c22 ^i — pi t 0] [ ~~^i t^22 _ cl 722 Hs = [0 ] = <- [o] [J*] = 0 0 0 0 Mi tgi '[FX] [Fin . [o] [0] [Js} = 2 [Fin "o o" 0 0 [0]" [0] w = [J»] = 0 0 d 922 d g\i '[0] [0] [Fin .[o] [0] [0] o" 0 0 [0] [0] [0] APPENDIX G : MODIFIED EULERIAN ROTATIONS The transformation matrix corresponding to rotation from frame LA to frame LB in terms of the modified Euler angels with the 2-3-1 rotation sequence (i.e., (3,7, a) has the form c/?c7 57 —cac/3s^ + sas/3 cacy sac(3s^ + cas/3 —C75/? casf3s^j + sac/3 —sac^f —sas/3si + cac/3 The angular velocity can be written as u = W 0, where 9 is a vector of the rotation angles and W 1 57 0 0 cac7 sa 0 —sacy ca The partial differentials of W with respect to the rotation angles are : 1^1 oa 0 0 0 0 —sac7 ca cacy —sa .0 [—] 0 1$ 0 0 cy 0 —cas^i 0 0- 50:57 The direction cosines of the x,y and z coordinates of a frame LA in frame LB is given by the columns of the c/?c7 n = { sas/3 — cac/3s^/ V ; cas/3 + sac/3si J denoted by n, 2 and m, respectively : ( 57 I = < 00:07 ( > ; I —50:07 > 275 —c^s/3 m = < cas/3s^j + sac/3 \ —sas(3si + cac(3 Partial differentials of n, I, m with respect to a, /3 and 7 are as follows 0 ' dn .dn .dn x x —s/3c) cas/3 + sac/3s^/ sac/3 + —sas/3 + cac0s"f 'dm dm dm Xf Xf x —cf3s^i cases') cac/3 — sasj3s^/ 0 0 cy -sac) 0 —casi -cacq 0 sas^ . —cac/Sc) sac/3cq 0 — eye/? s^s/? -sas/3s^f + cacjS cac/3s~f — sets/? cas/3c^ .—cases') — sac(3 —sac/3s*j — cas/3 {|^} = \ -sas/3 + cac/3sq > , dff} -cas/3 — sac/3s^j ) -sas/3 + coccus') i c a c ^ _ sasf3si K —sac/3 — cas(3s~i 0 { — } = —sas/3c). —c/?C7 n. -sas/3 + cac/3s^i . —cas/3 — sac(3s~j , -cas/3 — sac/3s^ ^ 7 ^ dl dt dl dp d/3dy 2 dad/3- 2 0 dm . dm 0 2 2 f 2 —cases') — sac/3 x sas/3si — cac(3 dm _ ^dadj ~ -sasfis^ + cac/3 —sac/3s7 — cas/3 —cacfis^ + sas/3 cys/? 0 2 x ^da~d/3 {—} —cas/3s^j — sac/3 sas(3si — cac/3 -cas(3si — sac/3 , {—}
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Mathematical modelling of flexible multibody dynamics...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Mathematical modelling of flexible multibody dynamics with application to orbiting systems Ibrahim, Ahmed El-Hady M. 1988
pdf
Page Metadata
Item Metadata
Title | Mathematical modelling of flexible multibody dynamics with application to orbiting systems |
Creator |
Ibrahim, Ahmed El-Hady M. |
Publisher | University of British Columbia |
Date Issued | 1988 |
Description | A relatively general formulation for the governing equations of motion, applicable to a large class of flexible multibody systems, is developed using a concise matrix format. The model considered consists of a number of arbitrarily connected flexible deployable members forming branched and closed loop configurations. Joints between bodies are permitted up to six degrees of freedom in translation and rotation. To be effective, the matrix-Lagrangian formulation necessitates development of the kinetic energy expression in a quadratic form in terms of the system velocities. The mass matrix associated with such a quadratic form is known for simple systems such as a collection of point masses, a group of connected rigid bodies, and a discretized flexible structure. However, for a multibody system, where the contributing forces arise from system's translation, rotation, elasticity, deployment, and their interactions, such an expression is not available. To fill this gap, multibody kinematics is developed in terms of the elements of the geometry matrix, which uniquely describes the configuration of branched systems. The characteristic dynamical quantities, i.e., elements of the mass matrix, are identified and the formulation is approached in an increasing order of complexity. The concept of specified and generalized coordinates together with established procedures of analytical dynamics lead to characteristic quantities ( Lagrangian, Hamiltonian, etc. ) and finally result in governing equations of motion which are new to the multibody dynamics. To account for flexibility in a consistent manner, a second-degree nonlinear displacement field is permitted. Alternatively, a linear displacement field can be used if the nonlinear terms up to the fourth-degree are preserved in the strain energy. An algorithm for calculating the stiffness matrix of a flexible element is developed, where terms up to the third-degree of nonlinearity in displacement are retained. Application of this versatile formulation is illustrated through a set of examples of contemporary interest. They pertain to a spacecraft comprising of a central rigid body with attached flexible appendages. The configuration corresponds to a large class of present and planned communication satellites. It can also represent the Space Shuttle based deployment of beam and plate type appendages aimed at scientific experiments or construction of the proposed Space Station. The system static equilibrium and stability are discussed. A computer code is developed and specialized to the specific cases in hand. Typical results of an extensive parametric study are presented for two particular situations : (i) the Space Shuttle based deployment of a beam or a plate type structural member; (ii) the configuration similar to the Waves In Space Plasma (WISP) experiment jointly proposed by Canada and the U.S.A. The problems are analyzed systematically, through progressive introduction of complexity, to help appreciate interactions between librational dynamics, flexibility, deployment, inertia parameters, orbit eccentricity, initial conditions, appendage orientation, etc. The information is fundamental to the missions concerned and essential to help develop appropriate control strategies. |
Subject |
Equations of motion Artificial satellites -- Orbits |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-09-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0098059 |
URI | http://hdl.handle.net/2429/28840 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1988_A1 I27.pdf [ 16.46MB ]
- Metadata
- JSON: 831-1.0098059.json
- JSON-LD: 831-1.0098059-ld.json
- RDF/XML (Pretty): 831-1.0098059-rdf.xml
- RDF/JSON: 831-1.0098059-rdf.json
- Turtle: 831-1.0098059-turtle.txt
- N-Triples: 831-1.0098059-rdf-ntriples.txt
- Original Record: 831-1.0098059-source.json
- Full Text
- 831-1.0098059-fulltext.txt
- Citation
- 831-1.0098059.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0098059/manifest