Stress Stiffening and Dynamic Stress Computation in Flexible Multibody Dynamics: Formulation Ragnar Ledesma Center for Automotive Stirtural Durability Simulation Technical Report UM-MEAM-95-06

Stress Stiffening and Dynamic Stress Computation in Flexible Multibody Dynamics: Formulation Ragnar Ledesma Center for Automotive Structural Durability Simulation The University of Michigan 2216 G.G. Brown Ann Arbor, MI 48109-2125 Abstract A formulation for the dynamics of flexible multibody systems is presented in this paper. This formulation relies on the use of floating reference frames to describe the configuration of the multibody system. Component flexibility is described in terms of the fiite element deformation coordinates. The equations of motion are derived through the generalized d'Alembert principle, and the resulting set of differential-algebraic equations are reduced to ordinary differential equations through the use of the augmented Lagrangian penalty method. Modal reduction is utilized to reduce the dimension of the deformation coordinates. Stiffening effect is included through the use of a stress stiffening matrix, which is computed efficiently as a linear combination of constant stress stiffness matrices with time-dependent scalar coefficients. The same formulation is used to devise an efficient method for computing the dynamic stresses which include the effects of inertia due to gross motion and elastic deformation.

1. Introduction The field of dynamics of flexible multibody systems has been gaining importance in the CAE community, especially in the aerospace and automotive sectors. This discipline has a broad base of applications, such as precision mechanisms, large space structures, automotive suspension systems, high-speed robotics, and aircraft structures. The increasing popularity of this field is found evident by the enormous amount of papers dealing with various aspects of multibody dynamics during the last decade alone. Recent advances in this field include developments in the solution of multibody dynamics differential-algebraic equations of motion1,2, inverse dynamics of flexible multibody systems3,4, nonlinear material behavior characterizations,6, modeling of geometric stiffening effects in flexible multibody systems7-10, modal reduction techniques in multibody dynamics11'12, modeling of contact and impact conditions in flexible multibody dynamics13-5, dynamic simulation using parallel computer architecturesl6-18, modeling of composite materials in multibody systems19'20, assignment of modal damping in flexible multibody dynamics21, among others. This paper presents a formulation of the dynamics of flexible multibody systems wherein stress stiffening effects are included in the formulation. The same formulation is used to devise a method for efficiently computing the dynamic stresses in each component of the flexible multibody system. The kinematic quantities are derived through the use of floating reference frames associated with each component of the multibody system. The deformation field in each body is described in terms of finite element nodal deformations, and the dimension of the deformation coordinates is reduced through standard modal reduction schemes. The generalized d'Alembert principle is employed to formulate the equations of motion, and the augmented Lagrangian penalty method is used to reduce the system of differential-algebraic equations to a set of ordinary differential equations. 2. Finite Element Formulation In this section, the equations of motion for a flexible multibody component are derived through the generalized d'Alembert principle. A floating reference frame is employed to describe the deformed state of the multibody component. The finite element method is employed for the spatial discretization of the flexible multibody component and for describing the component flexibility through the nodal elastic degrees of freedom. A practical consequence of this approach is twofold: multibody components of any geometry can be easily modeled; and existing commercial finite element codes for structural analysis

can be easily modified in order to make them applicable for the simulation of flexible multibody systems. Furthermore, the elastic forces are expressed as linear functions of the nodal deformations in order to make the method applicable to structural models that have already been linearized in a separate finite element code prior to input to the flexible multibody dynamics code. This premature linearization leads to a spurious softening of the structure during the dynamic simulation. Corrections to the prematurely linearized equations of motion are presented in Section 3. 2.1 Kinematics Consider an n -body flexible multibody system such as that shown in Fig. 1. A typical multibody component, say body i, is shown in Fig. 1 along with the floating reference frame associated with that body. The generalized coordinates consist of rigid body coordinates qi which describe the position and orientation of the floating reference frame associated with each multibody component, and deformation coordinates q1 which describe the deformation of the flexible body with respect to its floating reference frame. The rigid body coordinates qi consist of Cartesian coordinates Ri which describe the position of the origin of the floating reference frame associated with body i, and a set of Euler parameters &i which describe the orientation of the floating frame. The decision to use Euler parameters among several choices of orientation coordinates is based on the fact that Euler parameters form a minimal set of singularity-free orientation coordinates. The deformation from the nominal (rigid body) configuration is assumed to be small, so that the principle of local superposition22 is applicable. With the above choice of coordinates, the position of an arbitrary point P in body i is given by ri = Ri + Ai ui (1) where Ri is the location of the origin of the body axes with respect to the inertial frame, ui is the location of point P with respect to the body axes, and Ai is the rotation transformation matrix from the body axes to the inertial frame. In the three-dimensional case, the rotation transformation is given by [ 2(00+01)- 1 2(012 - 003) 2(13 + 0002) A= 2(0102 +0003) 2(0+02)- 1 2(20 -000) I (2) 2(0103- 0002) 2(0203 + 0001) 2(02+02)- 1 1

where the orientation coordinates are represented by four Euler parameters O', O 1, O2, and 0( which satisfy the following identity:23 j(ik)2=1. (3) k=O The position vector with respect to the body axes, ui, can be decomposed into Ui = ui + ui (4) where ui is the position vector of point P in the undeformed state with respect to the body axes, and uf is the deformation vector of point P with respect to the body axes. Differentiating Eq. (1) with respect to time yields the velocity vector of point P i'= Ri +i' ui +Ai tif (5) where the superposed dot represents differentiation with respect to time. In deriving Eq. (5), we have used the fact that vector ui is constant in time. The time derivative of the rotation transformation matrix can be expressed as the matrix product Ai = Ai (6) where the tilde symbol refers to the skew-symmetric matrix representation of the vector product operator, and xi is the angular velocity of the body axes whose components are measured with respect to the body axes. The angular velocity vector oi can be determined from the time derivatives of the Euler parameters through the following relation: oi= Gi 0 (7) where Gi is a matrix that depends linearly on the Euler parameters and is given by -01 00 03 -02 2 -02 -03 00 01 (8) -03 02 -1 00o The deformation vector u} can be expressed in terms of the nodal deformations by using a finite element discretization scheme, i.e., u= Ni qi (9)

where Ni is the shape function matrix and qi is the nodal deformation vector. Since the shape function matrix is time-invariant, the time derivative of the deformation vector becomes Uif= Ni qf (10) where qf is the nodal elastic velocity vector. Substituting Eqs. (6) through (10) into Eq. (5), we obtain the following expression for the velocity vector of point P: i' = Li Qi (11) where the velocity coefficient matrix Li is given by Li= [ 13 - A u G A N ]i (12) and qi is the vector of generalized coordinates for body i, i.e., qi= 0 (13) Differentiating Eq. (11) with respect to time yields the acceleration vector of point P namely, i'=L _ + A[ci]2ui+2A coui (14) The kinematical quantities derived above will be used in formulating the dynamics equations of motion in the next sub-section. 2.2 Dynamics In this section, we use the generalized d'Alembert principle to derive the equations of motion for a flexible multibody component. The generalized d'Alembert principle for an unconstrained flexible multibody component can be written as fi.Sri dV- I ii Sri pi dV=O (15) where f' is the vector of active forces (including elastic forces) per unit volume acting on the flexible multibody component, pi is the mass density, and Vi is the volume of the

component, and 8ri is the virtual displacement of the generic point P. The virtual displacement of point P is the first variation of the position vector which can be induced from the expression for the velocity vector in Eq. (11). Hence, the virtual displacement vector can be expressed as -ri = Li 8qi (16) The first integral in Eq. (15) represents the virtual work done by the active forces. This term can be decomposed into the sum of the virtual work done by the applied external forces and the virtual work done by the internal elastic forces, i.e., fi.Sri dV=f fiext6ri dV+I fint ri dV (17) If the applied external force vector is expressed in terms of components along the body axes, the virtual work done by the applied external forces can be written as J fiext ri dV = {(q iT [Li]T Ai fxt dV (18) The virtual work done by the internal elastic forces can be further decomposed into the sum of the virtual work done by the elastic forces associated with the structural stiffness and the virtual work done by the elastic forces associated with the geometric stiffness, i.e., I fint'6ri dV = -({8qT (kif + k q (19) where kff is the linear stiffness matrix of the flexible multibody component corresponding to the nodal defonrmation vector q and kg is the geometric stiffness matrix which depends on the geometry, displacement field, and state of stress of the elements representing the deformable multibody component. A discussion on the efficient computation of the motioninduced geometric stiffness is deferred until the section 3. We note that the elastic forces are due only to the elastic deformation and are therefore independent of the configuration of the floating reference frame. Combining Eqs. (12) through (19), and taking into account the fact that the generalized coordinates qi are independent for an unconstrained multibody component, the equations of motion for an unconstrained flexible multibody component are thus obtained:

MRR mR mRf ii FRi 03i QR i moR min mof = Fo I - 041 +1 QI (20) mfR mfo mff J Ff Pf QfJ or, alternatively, the equations can be written in the more compact form Mi qi = Fi- pi + Qi (21) where M1 is the component mass matrix, Fi is the external force vector, Pi is the elastic force vector, and Qi is the quadratic velocity force vector which absorbs the Coriolis and centrifugal force components coming from the second and third terms on the right hand side of Eq. (14). For convenience, the superscript i will be dropped in the remainder of the formulation. The elements of the mass matrix can be expressed in terms of the so-called finite element invariants which depend only on the undeformed geometry of the body and on the assumed displacement field characterized by the shape function matrix Ni. The finite element invariants, which are computed only once in a pre-processor module, are expressed by the following integrals:24 l = xk p dV;k=1,2,3 (22) Z = xj xk p dV; j,k=1,2,3 (23) Zk= NkpdV; k=1,2,3 (24) Zk =xj N k dV; j,k=1,2,3 (25) zJIk = NT Nk p dV;j,k=1,2,3 (26) 5 v where xk is the component of the undeformed position vector ur along the kth body axis, and Nk is the kth row of the shape function matrix. With the finite element invariants shown above, the components of the mass matrix are given by the following expressions: mRR = m I3 (27)

mRo =-AS1 G (28) mRf= A Z3 (29) moo = GT [Jo + J + J2] G (30) mof = GT [Jfc + Jft] (31) mff= Z1 + z522 + Z533 (32) where m in Eq. (27) is the total mass of the component, and the tilde symbol above the vector S1 in Eq. (28) refers to the skew-symmetric operator. The vector SI and the matrices Jfc, Jft, Jo, J1, and J2 are given by S1 = Zi + Z3 qf (33) z23 - z32 Jfc = Z31 - Z3 (34) z2- z21 qT? (z23 - z2) Jft = qfT (z31 z3) (35) LqT (z12 z21) (z222 + z23) _44 _2 z3 Jo = Z22 (Z23 + Z41) -Z23 (36)'31 -z4 (Z74 +z~~ (p22 + p33) p12 2p13 J1 = p21 (p33 + pll) p23 (37) p31 _p32 (p11 + p22) (q22 + q33) _q12 _q13 J2 = q21 (q33 + qll) q23 (38). q31 q32 (qll + q22) in which

pij = (Z + Z) qf;i,j = 1,2,3 (39) and qij = qT Zi qf;ij = 1,2,3 (40) The quadratic velocity force vectors are given by QR =- A [c S + 2Z3 qf] (41) Q = - GTS2 2GTJr (42) 3 3 Qf= E aij [Z i + Zij 4f I i=1 j=1 (43) - 2 o1(Z32 - Z23) + o,(z3 - Z 1) + _o3(Z21 - z2)] qf where aij are components of the (3 x 3) matrix Co, and the matrix Jr is given by (r22+r33) -r12 -r13 Jr =. r21 (r11+r33) r23 (44) r31 - r32 (r11+r22) in which rij q zji + q f Zij qf; i,j=1,2,3 (45) and the components of the vector S 2 are given by 3 (S2)i = ( Z + pij + qij) oj; i=1,2,3 (46) j=1 The generalized external force vector will depend on the type of load that is applied to the flexible multibody component. For example, if a distributed force vector fb, whose components are measured with respect to the body reference frame, is applied to the flexible multibody component, the generalized force vector due to this load can be determined directly from Eq. (18), which yields FR = A IfbdV (47) FO=GTfI u fb dV (48)

Ff= NTfbdV (49) Similarly, if a concentrated load vector fp, whose components are measured with respect to the body reference frame, is applied to a particular point P in the flexible multibody component, then one can obtain the generalized force vector due to the concentrated load by direct application of the principle of virtual work to yield FR = A fp (50) F = GT upfp (51) Ff = NTfp (52) where up is the position vector of point P measured relative to the body reference frame, and Np is the shape function matrix evaluated at point P. The generalized force vectors due to other types of loads such as distributed moments or concentrated moments can also be derived from the principle of virtual work. 3. Stress Stiffening Effects The premature linearization of the elastic forces in the equations of motion leads to a spurious softening of the structure during the simulation. The spurious softening is present for all angular velocities and causes numerical instability in the simulation when the angular velocity reaches a critical value. This critical value, of course, does not exist in reality but is a mere consequence of the premature linearization of the elastic forces. A simple and inexpensive procedure for avoiding spurious softening is to introduce stiffening effects that are due to an initial state of stress in the structure. The stiffening effects on the structure can be expressed as elastic forces that are obtained as the product of the initial stress stiffness matrix and the elastic deformation. Hence, the total internal elastic forces can be expressed as the sum of elastic forces associated with the linear structural stiffness matrix and elastic forces associated with the stress stiffness matrix, Pf = kff qf + kg qf (53) where kff is the linear structural stiffness matrix and kg is the initial stress stiffness matrix. The initial stress stiffness matrix is a function of the reference state of stress and the assumed displacement field for the elastic deformation, and is defined by the following integral:25

kg aNT a* aN dV (54) where the space-dependent functions aN arise from the nonlinear terms in the straindisplacement relation and the (9 x 9) matrix a* is a function of the reference state of stress. The reference state of stress is chosen as the state of stress in the structure due to the distributed inertial forces, constraint reaction forces and external forces that are present when the structure is subject to the nominal rigid body motion. This nominal motion is obtained by neglecting terms involving elastic displacements in the expression for the acceleration vector in Eq. (14). Since the reference state of stress due to the nominal motion varies with time, determining the stress stiffness directly from Eq. (54) will be computationally prohibitive. However, if we assume that the total stresses in the structure remain in the linear elastic regime, we can use the principle of superposition to compute for the time-varying stress stiffness matrix in a very efficient manner. This procedure is presented in the remainder of this section. To start the discussion on the efficient computation of the stress stiffness matrix, we note that the stress stiffness matrix depends only on the assumed displacement field and on the state of stress in the structure. The stresses in turn depend on the elastic deformation present in the structure. If we choose the stresses in the structure as those due to the nominal motion of the structure itself, then the forces that act on the structure are forces that arise from inertia, forces coming from the constraints, and externally applied forces. These forces form an equilibrated set of forces with the internal stresses in a quasi-static equilibrium. We can therefore determine the nominal deformation and nominal stresses in the structure through a sequence of linear quasi-static analyses. Furthermore, if we consider only one flexible multibody component at a time, the principle of virtual work applied to the flexible structure can be expressed as ffint.6r dV +f fext.r dV- anom. r p dV = (55) where fxt refers to the external forces applied to the structure and anom is the nominal acceleration vector obtained by dropping the terms due to the elastic deformation in the expression for the acceleration given in Eq. (14), i.e., anom = iro + a ur + c o( ur (56)

where r0 is the acceleration of the origin of the body reference frame, co is the angular velocity of the body reference frame, a is the angular acceleration of the body reference frame, and Ur is the position vector of a generic point in the structure. We note that the external forces include those coming from the reaction forces at the joints, so that the second term in Eq. (55) represents the virtual work done by the reaction forces and applied forces. By using a finite element discretization scheme, and assuming linear elastic behavior, the principle of virtual work expressed by Eq. (55) can be transformed into the following set of linear quasi-static equilibrium equations: kfflunom= P + R (57) where kff is the conventional linear stiffness matrix, Unom is the vector of nodal displacements due to the nominal inertia forces, applied forces and reaction forces, R is the time-varying nodal force vector due to the applied forces and reaction forces at the joints, and P is the time-varying nodal force vector due to the distributed nominal inertia forces and is given by P = - NT (i + a Ur + x Ur) p dV (58) The time-varying nodal force vector P can be expressed as a product of a constant matrix and a time-varying force vector, i.e., P=[- NTBpdV]{c} (59) where the space-dependent matrix B is given by B =[13 B1B2B3] (60) in which 13 is the (3 x 3) identity matrix and the sub matrices B1, B2, and B3 are given by 0 X3 -X2 B1 = -X3 0 xl (61) x2 -xl 0 O -X 1 -XI1-I B2 = -x2 -x2 (62) -x3-x3 0 [X2 0 x3-1 B3= xx3 0 (63) _Ox2 x1

where xl, x2, and X3 are the coordinates of a generic point in the structure, measured with respect to the body reference frame. The time-dependent vector { c I is defined as {C}T= [rl r2 r3 ( 2 (X3 (O2 c02 0o0CO2 (02 03 (03C01l (64) where the components of iP, a and wco are measured with respect to the body reference frame. Considering the linear quasi-static equilibrium equations written in Eq. (57), the nodal force vector P can be expressed as a linear combination of constant nodal force vectors 12 P = X ci Pi (65) i=1 where each of the constant nodal force vectors Pi (i = 1,... 12) are obtained from Eq. (59) by setting the element in the ith row of { c } equal to unity and all other elements to zero. Likewise, the nodal force vector R can be expressed as a linear combination of constant nodal force vectors m R = X i Ri (66) i=1 where Ri is the constant nodal force vector corresponding to a unit value of the ith concentrated force component, Xi is the value of the ith concentrated force component obtained during the dynamic analysis, and m is the number of concentrated force components associated with the mechanical constraints and external forces applied on the flexible multibody component. Herein, the concentrated force components also include concentrated moments associated with applied torques or torques that are generated by mechanical constraints. By assuming linear elastic behavior in the entire structure and by using the superposition principle, the deformation vector and the stress state vector due to the inertia forces and reaction forces can also be expressed by linear combinations similar to that of Eqs. (65) and (66), 12+m Unom = Ci (67) i=1 12+m anom = X Cib i (68) i=l where c12+i = Xi (i = 1,... m ), Unom is the time-varying nodal deformation vector due to the nominal motion, Gnom is the time-varying vector of stress states due to the nominal motion, and ui and ai (i = l,... 12+m) are constant nodal deformation vectors and

constant stress state vectors corresponding to each of the constant nodal force vectors Pi and Ri. By virtue of the fact that the stress stiffness matrix as defined in Eq. (54) depends only on the state of stress and on the assumed displacement field, the stress stiffness matrix associated with the nominal inertia forces, reaction forces and external forces can be computed as the linear combination of constant stress stiffness matrices, i.e., 12+m kg= ci (kg)i (69) i=l 1 where each of the constant stress stiffness matrices (kg) (i = 1,... 12+m) is evaluated from Eq. (54) with the (9 x 9) stress matrix ai, which is associated with the stress state ai, replacing the (9 x 9) stress matrix ac* which is associated with the reference stress state ca. These constant stress stiffness matrices are computed only once in a preprocessor module before the solution of the dynamic equations is started. Eq. (69) is based on the assumption that the deformation remains small. It is important to recognize the tremendous savings in computation time that Eq. (69) provides when compared to the direct computation of the stress stiffness matrix from Eq. (54). Eq. (69) allows us to determine the time-varying stress stiffness matrix without having to evaluate the integral in Eq. (54) at every time step. Instead, all that is needed to compute for the time-varying stress stiffness matrix is the summation of constant matrices, each of which is multiplied by a scalar quantity which varies at every time step. 4. Augmented Lagrangian Penalty Method of Solution Consider now the entire multibody system as an assembly of the individual multibody system components. The configuration of the multibody system can still be described by the generalized coordinates (R, 0, q f) associated with each multibody component, and the generalized d'Alembert principle expressed in Eq. (15) is still valid. However, the generalized coordinates (R, 0, qf) are no longer independent because the motion of specific points in different bodies are related according to the type of mechanical joints that interconnect them. This interdependence of the generalized coordinates can be expressed in the form of a vector of kinematic constraint equations, such as <b(R, 0, qf, t)=O0 (70) Because of the fact that the generalized coordinates are not independent, the generalized d'Alembert principle of Eq. (15) can not be used directly to generate the equations of

motion as in the case of a single multibody system component. To generate the equations of motion for the entire flexible multibody system, we use the augmented Lagrangian penalty method26. The augmented Lagrangian penalty method of formulating the equations of motion of constrained multibody systems is essentially a synthesis of the regular penalty method and the classical Lagrange multiplier method. In this method, a "dynamic penalty system" associated with each constraint equation is appended to the equations of motion corresponding to the unconstrained system, and Lagrange multipliers are also added as correction terms to the penalty system. The resulting equations of motion for the constrained flexible multibody system will thus take the form similar to that of Eq. (21) but will also include penalty terms and Lagrange multipliers associated with the constraint equations: Mrq c +T a[+ 2[.Qt + Q]+ ~q X = F- P+ Q (71) where q is the transpose of the constraint Jacobian matrix, X is the vector of Lagrange multipliers, ac is a diagonal matrix of penalty numbers, and the quantities inside the brackets characterize the dynamic penalty system associated with each of the constraint equations. The first and second time derivatives of the holonomic constraints are given by ~ -= Dq + Dt (72) ID = (D~q q + (Dq l + (1t (73) where the subscript denotes partial differentiation and the superposed dot denotes the total time derivative. The augmented Lagrangian penalty method requires an iteration procedure to compute for the correct value of the Lagrange multipliers. This iteration can be performed simultaneously with the Newton iteration or modified Newton iteration during the solution of the nonlinear equations at each time step in a numerical time integration scheme. The iterative equation for the Lagrange multipliers is given by A(i+l) =(i)++ 2 Q4 2+ ] (74) The iterative process described by Eq. (74) involves only a few additional operations during each iteration but it significantly improves the convergence of the solution of the dynamics equations of motion compared to that of the regular penalty method. 5. Modal Reduction and Component Mode Synthesis

The finite element formulation of the dynamics of flexible multibody systems as presented in the previous sections is computationally expensive because of the large number of elastic degrees of freedom in the finite element model. Except for very trivial cases, finite element models of multibody components normally contain in the order of thousands to tens of thousands degrees of freedom. It is therefore necessary to formulate a scheme for reducing the number of degrees of freedom in order to implement a flexible multibody dynamics code in a manner that is economically viable. Towards this end, a modal reduction scheme is utilized to minimize the number of elastic degrees of freedom in the flexible structure. Although at first glance, the superposition of different modes of vibration in a nonlinear problem appears to violate the principle that mode superposition is valid only for linear problems, the validity of mode superposition is in fact not violated in the present formulation because the elastic deformations are measured from a floating reference frame, and therefore the deformations remain small throughout the motion. The mode superposition principle can therefore be considered as a similarity transformation that implements a change of basis from the finite element deformation coordinates to the modal deformation coordinates. It should be noted that for highly nonlinear problems, the instantaneous frequencies of the system exhibit significant changes continuously, and the use of dynamic mode superposition would not be as efficient as the direct time integration methods. However, for dynamic systems with local nonlinearities, the use of mode superposition in conjunction with the pseudo-force method offers an attractive alternative to direct time integration methods27. This is particularly true in the dynamics of flexible multibody systems where most of the nonlinearities are confined to the rigid body coordinates and at the finite element nodes located at the attachment points. Each multibody component is naturally treated as a substructure with interfaces at the mechanical joints, and the behavior of the dynamic system can be accurately predicted through component mode synthesis. Furthermore, since the elastic deformation of each component is known to remain small throughout the motion, the mode shapes based on the initial tangent stiffness can be used as the basis for the elastic displacements throughout the motion. In cases where geometric stiffening effects are significant, the mode superposition procedure can be used with the pseudo-force method where only the residual forces due to geometric nonlinearities need to be transformed at each time step. There are two methods that can be used in implementing modal reduction in flexible multibody dynamics, depending at which stage of the formulation that modal reduction is introduced. In the first method, modal reduction is introduced after the finite element

invariants have been computed. In the second method, modal reduction is introduced at the outset, hence, modal invariants are computed in lieu of finite element invariants. The first method is preferred since the finite element invariants are based on a formulation that is consistent with the displacement interpolation used in generating the structural stiffness matrix. The result is that the consistent mass matrix and consistent quadratic velocity force vectors are preserved in the modal reduction. The advantage of this method is that the rotary inertia properties and the exact mass distribution are captured automatically. On the other hand, the disadvantage of this method is that it can be used only if the analyst has access to the coding of the finite element analysis software. More common in practice, the analyst does not have access to the coding of a commercial finite element software. In this case, the first method is not an option, and the analyst has to use the mode shape vectors obtained from the finite element software to approximate the modal invariants. In addition, the mass distribution in the flexible multibody component can be characterized only through the lumped masses at the nodes. This approximation will yield valid results as long as the highest frequency that is excited during the motion has a wavelength that spans across several nodes. Otherwise, this limitation can be mitigated by making a reasonable estimate of the rotary inertia lumped at each node. Details of the two methods of implementing modal reduction are discussed in the remainder of this section. 5.1 Modal Reduction: Consistent Mass Formulation Consider again the equations of motion for the flexible multibody system as expressed in Eq. (71). The equations of motion can be written in the following partitioned form: [fr m ][ q~~f ]+ ( D X ( + a [ 2 + 2 CL ( D + n2+]= ] [ ]+[ (75) mfr mff rq T Ff Pf where the generalized coordinates (qr, qf) denote the rigid body coordinates and finite element deformation coordinates, respectively. We now introduce a change of basis from the n finite element displacement coordinates top modal coordinates, with p <<n, through the following linear transformation: qf= [-d s] [ ] (76) where the columns of the sub-matrix id consists of the retained dynamic modes and the columns of the sub matrix yrs are the appropriate static correction modes, and rj is the

vector of modal deformation coordinates. The above combination of dynamic and static modes define the new basis for the elastic deformation. These modes are computed only once at the start of the simulation and are based on the undefonnrmed configuration of each of the flexible multibody component. Applying the transformation given by Eq. (76) to Eq. (75) and pre multiplying the second set of equations in Eq. (75) by the transpose of the modal transformation matrix yields the reduced set of equations of motion in terms of the rigid body coordinates and the modal deformation coordinates, [ WTm, WT ] I[.. ] ~+ [ X + a D[ + 2 4 Q (D + Q = v Tnfr vt Tmfff J LTT JL T(Ff-Pf+Qf) (77) It is important to note that modal mass matrix is a full matrix when static correction modes are used. Furthermore, the different modes are coupled through the constraint reaction forces and through the force vectors on the R.H.S. of Eq. (77). Therefore, direct integration methods are still needed to integrate Eq. (77), but the reduced system is much smaller than the original system. Finally, we observe that the consistent mass formulation in the mass matrix and in the quadratic velocity force vector is preserved because the modal transformation is applied after the finite element invariants have been computed. 5.2 Modal Reduction: Lumped Mass Formulation A common situation that a multibody dynamics analyst finds in practice is that the structural data for a flexible multibody component is available only in the form of finite element output from a previous linear dynamic analysis. The structural dynamics analysis will normally provide the natural frequencies and corresponding eigenmodes, as well as some selected static correction modes. At best, the dynamics analyst can obtain the mass distribution in the form of lumped masses at each of the nodes in the finite element model. As a consequence, the finite element invariants given by Eqs. (22) - (26) can not be evaluated exactly. One alternative is to approximate the finite element invariants by the socalled modal invariants. The modal invariants are computed by a summation of discrete variables instead of the integration of continuous variables as shown in Eqs. (22) - (26), with the lumped mass at each node replacing the continuous mass distribution and with the values of the discrete mode shape matrix evaluated at each of the nodes replacing the continuous shape function matrix. The modal invariants are given by the following summations:

N ZI = a xk mi; k=1,2,3 (78) i=l N Zjk = - xj xi mi; j,k=1,2,3 (79) i=l N Z= ikmi; k=1,2,3 (80) i=l N Z x ik mi; j,k=1,2,3 (81) i=1 N 7Z3k =, ~Tik mi; j,k=1,2,3 (82) i=l where N is the number of nodes in the finite element model, mi is the lumped mass at node i, xk is the coordinate of node i along the Wh axis of the body reference frame, and Vk is the (1 x p) row vector extracted from the mode shape matrix and corresponding to the deformation of node i along the kh axis of the body reference frame. The elements of the mass matrix and the elements of the quadratic velocity force vector are computed in a similar manner as presented in Section 2, with the modal invariants replacing the finite element invariants and with the modal deformation coordinates replacing the finite element deformation coordinates. The equations of motion are formulated by using the augmented Lagrangian penalty method and the system of equations are reduced by pre multiplying the equations associated with the elastic motion by the transpose of the mode shape matrix. The resulting reduced set of equations can be written as [rr m+[ q r + 2 (Fr+Qr) m[r mrr rl 1 IJr+r] +a + 2 LJ + nJ]=[F) () where the holonomic constraint equations are now expressed in terms of the modal deformation coordinates, i.e., O(R, 0, rl, t) = 0 (84) The elements of the generalized force vector F due to external forces can be derived from the application of the principle of virtual work, which results in the following:

FR = A j fi (85) i=l Fo = GTE ui fi (86) i=l N Fl=X hiTfi (87) i= where fi is the force applied at node i and whose components are measured with respect to the body reference frame, v' is the (3 x p) matrix extracted from the mode shape matrix and associated with the translation deformation of node i. In Eq. (86), ui is the instantaneous position vector of node i whose components are measured with respect to the body reference frame. The coordinates of this vector are determined from the equation Ui = X i + Nfi r (88) where xi contains the coordinates of node i in the undeformed configuration. Expressions similar to Eqs. (85) - (87) can be formulated for applied moments at the nodes. Returning to Eq. (83), we also note that the modal mass matrix man is a full matrix because the mode shape matrix v is a combination of dynamic modes and static correction modes. The generalized force vector due to internal elastic forces is given by P = [As + Ag] 11 (89) in which the generalized structural stiffness A, is a constant, block-diagonal matrix given by As = VTkff N' (90) and the generalized stress stiffness matrix Ag is a time-varying, full matrix given by A -= Tkg t (91) The modal reduction of the elastic forces associated with the stress stiffening matrix is valid as long as the elastic deformation remains small. By using the procedure described in Section 3, we can efficiently compute for the time-varying generalized stress stiffness matrix by expressing it as a linear combination of constant generalized stress stiffness matrices with time-varying scalar coefficients, i.e.,

12+m Ag = Ci fT (kg)i l (92) i=l Finally, we observe that the reduced set of equations shown in Eq. (83) are fully coupled through the system mass matrix, the constraint force vectors, and through the generalized force vectors. Therefore, direct time integration methods are needed to integrate the reduced system of equations. 6. Dynamic Stress Computation Dynamic stress computation has often been treated as a procedure that is separate from multibody dynamics. The current CAE practice in industry is to post-process the results of multibody dynamics by taking the time-varying reaction forces between rigid or flexible multibody components and using them as inputs to the structural dynamics analysis of each multibody component via the finite element method. The structural dynamics analysis of each component may be geometrically linear or nonlinear, depending on the magnitude of the motion that is imposed on some attachment points in the component. Structural dynamics analysis using large finite element models is computationally prohibitive except for very short duration events. This is especially true for complicated components which require a very fine mesh. As a consequence, the analyst usually resorts to performing a quasi-static analysis of a few selected load cases taken from the reaction force history provided by multibody dynamics. This procedure for stress computation is adequate for cases where inertial forces due to component flexibility are not significant. However, for cases where inertial forces are important, the stresses resulting from quasistatic analysis are not accurate when compared to a full dynamic stress analysis. An alternative to the aforementioned procedures for computing dynamic stresses is to make the dynamic stress computation an integral part of flexible multibody dynamics where the dynamic stresses are computed along with the flexible body motion at every time step during the simulation. Again, there are two methods of implementing the stress computation within the flexible multibody dynamics analysis, depending on whether the analyst has access to the coding of the finite element analysis software. If the analyst has access to the finite element source code, then the dynamic stresses can be determined directly from the stress-displacement relation a=CB rj (93)

where a is the stress vector, C is the constitutive stress-strain relation matrix, and B is the strain-displacement relation matrix for a particular element. The computation of the stresses at selected points of interest may be performed as a post-processing procedure at every time step once the elastic deformation has been determined. The operation described by Eq. (93) may be implemented as function calls to some modules in the finite element code. If the analyst does not have access to the finite element source code, which is the more common case, then the operation described by Eq. (93) can not be implemented. In this case, an effective solution for the computation of dynamic stresses is to employ a quasi-static procedure similar to the one described in the section pertaining to the stress stiffening effects. This procedure is based on the assumption that, in a quasi-static sense, the stresses in the flexible multibody component form an equilibrated set of forces together with the inertia forces, reaction forces, and applied forces. Hence, the quasi-static equilibrium equation for the flexible multibody component can be expressed in the form kffUf= fI + fR + fA (94) where Uf is the time-varying nodal deformation vector, and f I, fR, and fA are time-varying force vectors due to inertia, constraint reaction, and applied loads, respectively. The force vector due to inertia may be decomposed an inertia force due to the nominal motion of the body which is obtained by neglecting terms associated with the elastic deformation in the expression for the acceleration, and an inertia force due to the elastic deformation, i.e., f[ = fnom + fdef (95) where the inertia force due to the nominal motion is given by N fnom = - 2 mi (r'O + a ui + c C uO (96) i=l and the inertia force due to elastic deformation is given by fef = _ m i {(a + co o) Ski Ti + 2 c Wli 4 + /i (97)} i=l The lumped mass formulation is used in Eqs. (96) and (97) because the mode shape matrix provides the deformation only at the nodes. It can be shown that the inertia force due to the elastic deformation may be expressed as a product of a constant matrix and a time-varying vector, so that Eq. (97) may be written as

N fdef X = _ m,i i4 (b ) (98) i=l where Nig is an augmented (3 x 9p) matrix of mode shapes at a particular node, p is the number of modes, and the elements of the (9p x 1) time-dependent vector are given by bjkl = [a + C 2]jk rli + 2 []jk il + Sjk ql; j,k=1,2,3; l=l,...,p (99) where the oprator [. k denotes the element of the matrix in the jth row and kth column, and 6jk is the Kronecker delta operator. By assuming linear elastic stress-strain behavior and small elastic deformations, the time-varying force vectors in Eq. (94) can be expressed as a linear combination of constant force vectors with time-varying scalar coefficients. Furthermore, the constant force vectors due to nominal inertia forces, constraint reaction, and applied forces are the same as the constant force vectors previously derived in Section 3, so that the time-varying force vector in the quasi-static equilibrium equations of Eq. (94) may be written as 12+m 3 3 P fIom + fR + fA+ fef = i Pi + bjkl Pjkl (100) i=l j=1 k=l 1=1 where m is the number of concentrated force components due to applied forces and reaction forces, and p is the number of modes used to describe the elastic deformation. The time-varying scalar coefficients ci and bjkl are determined from the results of the flexible multibody dynamics analysis, and the constant force vectors Pi and pdkl are computed in a pre-processor module before the dynamic simulation is started. By using the principle of superposition, the dynamic stresses can be obtained in a similar manner, i.e., 12+m 3 3 P ((t)-= Ci =i + c Z E bjkl:k (101) i=l j=1 k=l 1=1 where the constant stress influence matrices ai and ej fare constant stress influence matrices which are computed only once in a pre-processor module. The operation described in Eq. (101) can be performed as a post-processing procedure at every time step during the dynamic simulation after the multibody dynamics equations of motion have been solved and reaction forces have been determined. Hence, the dynamic stress computation simply requires a superposition of (12+m+9p) stress influence matrices at each time step during the dynamics simulation. In contrast with stresses obtained from a conventional quasi-static analysis, the stresses given by Eq. (101) include inertia terms that arise from the gross

nominal motion and inertia terms that arise from the elastic deformation. Finally, to save on memory space and computation time, only the stresses at a few selected points of interest need to be computed. These selected points of interest may be determined from a previous quasi-static analysis wherein the applied forces are those coming from rigid multibody dynamics. 7. Summary A formulation for the dynamics of flexible multibody systems has been presented. The formulation is based on the use of floating reference frames to describe the configuration of the multibody system. The deformation field in each multibody component is discretized by using the finite element method. The equations of motion are formulated through the generalized d'Alembert principle, and the differential-algebraic system of equations describing the entire multibody system is reduced to a set of ordinary differential equations by using the augmented Lagrangian penalty method. The dimension of the resulting equations are reduced through a modal reduction of the deformation coordinates. The formulation also includes an efficient computational method for the stress stiffening effect which is needed to correct the spurious softening effect that arises from the premature linearization of the elastic forces. The stress stiffness matrix is computed as a linear combination of constant stress stiffness matrices which are computed only once in a preprocessor module. The same formulation is employed to efficiently compute the dynamic stresses in each component of the flexible multibody system. The dynamic stress are computed from a simple superposition of constant stress influence matrices which have been previously computed from a series of linear static analyses. The dynamic stresses thus obtained are comparable to those obtained from a full dynamic stress analysis but provides a substantial savings in computation time. References: 1. E. Bayo and R. Ledesma,'Augmented Lagrangian and projection methods for constrained multibody dynamics', Proceedings of the 23rd ASME Biennial Mechanisms Conference, Design Engineering Division, vol. DE-72, 237-247, 1994. 2. D. Bustamante, A. Kurdila and R. Menon,'Convergence and attractors of augmented Lagrangian formulations in the dynamics of multibody systems', Proceedings of the 34th Structures, Structural Dynamics and Materials Conference, 2076-2086, 1993.

3. R. Ledesma and E. Bayo,'A non-recursive Lagrangian solution of the non-causal inverse dynamics of flexible multibody systems: the planar case', International Journal for Numerical Methods in Engineering, vol. 36, 2725-2741, 1993. 4. R. Ledesma and E. Bayo,'A Lagrangian approach to the non-causal inverse dynamics of flexible multibody systems: the three-dimensional case', International Journal for Numerical Methods in Engineering, vol. 37, 3343-3361, 1994. 5. M. Xie and F.M.L. Amirouche,'Treatment of creep and nonlinearities in flexible multibody dynamics', AIAA Journal, vol. 32, 190-197, 1994. 6. J. Ambrosio and P.E. Nikravesh,'Elasto-plastic deformations in multibody dynamics', Nonlinear Dynamics, vol. 3, 85-104, 1992. 7. 0. Wallrapp and R. Schwertassek,'Representation of geometric stiffening in multibody system simulation', International Journal for Numerical Methods in Engineering, vol. 32, 1833-1850, 1991. 8. A. K. Banerjee,'Block-diagonal equations for multibody elastodynamics with geometric stiffness and constraints', Journal of Guidance, Control, and Dynamics, vol. 16, 10921100, 1993. 9. J. Ryu, S. Kim and S. Kim,'General approach to stress stiffening effects on flexible multibody dynamic systems', Mechanics of Structures and Machines, vol. 22, 157-180, 1994. 10. Z. Boutaghou, A. Erdman and H. Stolarski,'Dynamics of flexible beams and plates in large overall motions', ASME Journal of Applied Mechanics, vol. 59, 991-999, 1992. 11. H. Wu and N. Mani,'Modeling of flexible bodies for multibody dynamic systems using Ritz vectors', ASMEjournal of Mechanical Design, vol. 116, 437-444, 1994. 12. A. Liu and K. Liew,'Non-linear substructure approach fro dynamic analysis of rigidflexible multibody systems', Computer Methods in Applied Mechanics and Engineering, vol. 114, 379-396, 1994. 13. F. Amirouche, M. Xie, N. Shareef and M. Valco,'Finite element modeling of contact conditions in multibody system dynamics', Nonlinear Dynamics, vol. 4, 83-102, 1993.

14. H. Lankarani and P. Nikravesh,'Continuous contact force models for impact analysis in multibody systems', Nonlinear Dynamics, vol. 5, 193-207, 1994. 15. S. Ko and B. Kwak,'Frictional dynamic contact analysis in deformable multibody systems', Finite Elements in Analysis and Design, vol. 12, 27-40, 1992. 16. A. Eichberger, C. Fuehrer and R. Schwertassek,'Benefits of parallel multibody simulation', International Journal for Numerical Methods in Engineering, vol. 37, 15571572, 1994. 17. J. Chiou, K. Park and C. Farhat,'Natural partitioning scheme for parallel simulation of multibody systems', International Journal for Numerical Methods in Engineering, vol. 36, 945-967, 1993. 18. S. Ider,'Finite-element based recursive formulation for real time dynamic simulation of flexible multibody systems', Computers & Structures, vol. 40, 939-945, 1991. 19. J. Kremer, A. Shabana and G. Widera,'Application of composite plate theory and the finite element method to the dynamics and stress analysis of spatial flexible mechanical systems', ASME Journal of Mechanical Design, vol. 116, 952-960, 1994. 20. A. Ghazavi, F. Gordaninejad and N. Chalhoub,'Dynamic analysis of a composite material flexible robot arm', Computers & Structures, vol. 49, 315-327, 1993. 21. A. Lee,'Component modes damping assigment methodology for articulated, multiflexible body structures', Journal of Guidance, Control and Dynamics, vol. 16, 11011108, 1993. 22. R. Nickell,'Nonlinear dynamics by mode superposition', Computer Methods in Applied Mechanics and Engineering, vol. 7, 107-129, 1976. 23. P. Nikravesh and I. Chung,'Application of Euler parameters to the dynamic analysis of three dimensional mechanical systems', ASME Journal of Mechanical Design, vol. 104, 785-791, 1982. 24. A. Shabana, Dynamics of Multibody Systems, Wiley, 1989. 25. R. Cook, D. Malkus and M. Plesha, Concepts and Applications of Finite Element Analysis, 3rd edition, Wiley, 1989.

26. E. Bayo, J. Garcia de Jalon and M. Serna,'A modified Lagrangian formulation for the dynamic analysis of constrained mechanical systems', Computer Methods in Applied Mechanics and Engineering, vol. 71, 183-195, 1988. 27. K. J. Bathe and S. Gracewski,'On nonlinear dynamic analysis using substructuring and mode superposition', Computers & Structures, vol. 13, 699-707, 1981.

ri k~~~~~~~~~~~~~~~~q xi X3 Fig. 1 Multibody system and reference frames r i 1 X2 k I X2 3~~X 1%,~X x r 33 X3 Fig. 1: Multibody system and reference frames