A Nonlinear Viscoelastic Bushing Element in Multibody Dynamics Ragnar Ledesma Zheng Dong Ma Gregory Hulbert Alan Wineman Center for Automotive Structural Durability Simulation Technical Report UM-MEAM-95-05

A Nonlinear Viscoelastic Bushing Element in Multibody Dynamics Ragnar Ledesma Zheng-Dong Ma Gregory Hulbert Alan Wineman Center for Automotive Structural Durability Simulation The University of Michigan 2216 G.G. Brown Ann Arbor, MI 48109-2125 Abstract This report presents a formulation for nonlinear viscoelastic force elements in multibody systems. The formulation of nonlinear viscoelastic force elements is based on the assumption that the relaxation function can be expressed as a sum of functions which are nonlinear in the deformation and exponentially decreasing in time. These forces can represent elastomeric mounts or bushings in automotive suspension systems. The numerical implementation of the nonlinear viscoelastic bushing model into the a general purpose rigid multibody dynamics code is described, and an extension of the formulation wherein component flexibility is included is also presented. The validation of the formulation is carried out by developing two models for bushings, one using a nonlinear elastic model, and the other using a nonlinear viscoelastic model, and comparing the simulation results of the two models with experimental data. For this purpose, a simple system consisting of a lower control arm supported by bushings is built at the Center's test facility and subjected to a simulated road load event. Comparison with test results shows that the viscoelastic bushing model gives a better load prediction than the conventional nonlinear elastic bushing model.

1. Introduction The use of multibody system dynamics as a viable tool in computer-aided engineering has recently received wide acceptance in many engineering and manufacturing industries, especially in the automotive industry. The recognition of the potential of multibody dynamics as a CAE tool is due to tremendous strides in the development of this field of study as evidenced by the enormous amount of publications on the subject during the last decade. Recent developments in this field cover a wide range of issues including the choice of coordinates which describe the dynamic system1, the choice of methods for the formulation of the equations of motion2, the development of algorithms for the efficient numerical solution of the equations of motion3 4, the modeling of material properties and mechanical constraints5,6, the proper modeling of component flexibility7-l0, among others. The unifying objective of these developments is to enhance the ability of multibody dynamics simulation to predict the actual behavior of the dynamic system that is being modeled The ultimate goal is to have a general-purpose multibody dynamics code that can predict the behavior of a dynamic system during the design stage so that the perfonrmance of the system can be evaluated before expensive prototypes are built. In the ideal situation, the multibody dynamics code can help engineers determine the optimal design for the dynamic system and only one prototype of the final design is needed to verify the performance measures that have been predicted by the computer simulation. This report describes part of the research effort at the Center for Automotive Structural Durability Simulation. The research described in this report deals with the proper modeling of bushing forces between suspension system components. This study was motivated by engineers at Ford Motor Company whose experience have led them to the conclusion that the modeling of bushing forces in automotive suspension systems plays an important role in predicting the dynamic behavior of the suspension system. The choice of the bushing model is especially crucial in the prediction of the loads that act on the suspension system components that are supported by bushings. The accurate prediction of the dynamic loads that act on the suspension system components is important because these loads feed directly into the fatigue life prediction of the components. Considering the modeling of bushing forces in the multibody dynamics codes, the present state-of-the-art is to model the bushing forces as a Kelvin solid which is represented by a spring in parallel with a viscous damper. The spring force is a nonlinear function of the instantaneous deformation of the bushing and viscous damping force is a linear function of the instantaneous relative velocity of the two components that are connected by the bushing. One main drawback of this approach is that the nonlinear elastic bushing model provides a

dynamic response in which the energy dissipation is a linear function of the excitation frequency. Hence, the nonlinear elastic bushing model will damp out the high frequency content of the dynamic response and will, in general, capture the correct energy dissipation characteristic of the material at only one particular excitation frequency. Experimental results indicate that bushings, which are made of elastomeric material, exhibit a deformation history-dependent behavior which can be characterized as one with'fading memory'. These results suggest that the bushing response fits nicely into the theory of nonlinear viscoelasticity. The development of the nonlinear viscoelastic force model for the bushing response is discussed in a separate report. It is sufficient to note in this report that the history-dependent, nonlinear viscoelastic force model for the bushing force is characterized by a convolution integral where the kernel is a series of exponential functions which can capture the'fading memory' characteristic of the polymeric material in the bushing. Previous studies on the dynamic analysis of structures containing viscoelastic material have been reported in the literaturell-19. Most of these studies deal with the viscoelastic behavior at the material level and use finite elements or boundary elements to build up the viscoelastic structure and subsequently determine the dynamic response either through the frequency domain or Laplace domain, or through a time domain realization of the equations of motion in the Laplace domain. Because of the use of the finite element method or the boundary element method to model the viscoelastic structure, all of the aforementioned studies require expensive calculations and are not suitable for the dynamics of multibody systems which are made up of several rigid, elastic, or viscoelastic components. A notable exception is the work by Gaul and Chen18 who used the boundary element method to develop a relation between 12 stress resultants and 12 deformation variables which represent an elastomeric mount in a multibody system. That study, however, was limited to linear viscoelastic behavior so that the correspondence principle could be used to replace the elastic isotropic material behavior with the corresponding viscoelastic constitutive equations. In contrast, the present study efficiently captures the nonlinear viscoelastic behavior by characterizing the behavior of the elastomeric bushing not at the material level, but rather, at the component level where stress resultants are expressed as nonlinear functions of the deformation variables. This approach bypasses the expensive finite element or boundary element calculations, so that the method can be readily applied to multibody system dynamics simulation (at the expense, however, of some predictive capabilities). In the remainder of this report, we describe the formulation of the bushing forces as nonlinear viscoelastic force elements in multibody systems and the numerical

implementation of the formulation into the general-purpose multibody dynamics code ADAMS20. The verification of the formulation and code implementation is performed by building a simple multibody system at the Center's test facility and comparing simulation results with experimental data. The simple multibody system consists of an automotive suspension system component that is supported by bushings and subjected to prescribed dynamic loads and boundary conditions. As an additional verification, the proposed nonlinear viscoelastic bushing model is also compared with the conventional nonlinear elastic bushing model. Comparison between the two bushing models shows that the nonlinear viscoelastic bushing model provides a more accurate prediction of the experimental results. 2. Formulation of the Equations of Motion The bushing model developed at the Center takes the form of a set of uncoupled nonlinear viscoelastic force elements. In the context of multibody dynamics, these nonlinear viscoelastic forces are treated as massless nonlinear force elements that act between two parts or bodies which are connected by a bushing. The formulation and implementation of these massless nonlinear viscoelastic force elements can be best explained by first looking at the equations of motion for a constrained set of rigid bodies, and subsequently determine how the bushing forces enter into the equations of motion. Consider a simple rigid multibody system consisting of two bodies, i and j, which are subject to mechanical constraints. The equations of motion for the rigid multibody system may be written in the following descriptor form: [mRR mROe i R\+FR1() QR + FR (1) nR MOO tQo Fo subject to the holonomic constraint equations 4(R,) = 0 (2) In Eqs. (1) and (2), the generalized coordinates (R,0) refer to the rigid body translation and rigid body orientation, respectively. The rigid body coordinates describe the position and orientation of the body-fixed frames attached to each part or body in the multibody system, as shown in Figure 1. Also in Eq. (1), mRR and meo are the mass matrices associated with the translational mass and rotational inertia of the rigid body, respectively, and mRO is the

inertial coupling between the rigid body translation and rigid body rotation. The force vector Q is the generalized force due to external forces, and the force vector F is the generalized force due to quadratic velocity terms including centrifugal forces. The matrix [R*B] is the constraint Jacobian matrix, and the vector X is the vector of Lagrange multipliers associated with the constraints. The equations of motion, Eqs. (1) and (2), form a set of differential algebraic equations (DAE's) which are, in general, numerically more difficult to solve than ordinary differential equations (ODE's). Now consider the two bodies i and j to be connected to each other by a bushing whose attachment points are located at point Pi on body i and point Pi on body j as shown in Figure 1. If the bushing forces are treated as massless force elements acting between bodies i and j, the only way that these force elements can enter the equations of motion is through the generalized force vector Q by treating the bushing forces as external forces applied to the bodies which are connected by a bushing. In what follows, we will consider only bushing forces in the radial and axial directions to simplify the discussion. A similar procedure will apply for the conical and torsional bushing moments. If we consider a bushing force vector fi acting on body i and whose components are measured with respect to the bushing coordinate system ri - si - ti, the generalized force vector associated with the translational coordinates of body i due to the bushing force vector fi is given by Qi=AiApfi (3) and the generalized force vector associated with the rotational coordinates of body i due to the bushing force vector fi is given by QO=G!u Apfi (4) where Ap is the orthogonal transformation matrix from the bushing coordinate system to the body-fixed reference frame, Ai is the orthogonal transformation matrix from the bodyfixed reference frame to the inertial reference frame, and up is the position of point Pi measured with respect to the body-fixed reference frame. The tilde symbol above the vector ub represents the skew-symmetric matrix operator, and Gi is a matrix that maps the time

derivatives of the orientation coordinates to the angular velocity of the body, and is defined by the equation co=GOi = i (5) where aj is the angular velocity of body i and whose components are measured relative to the body-fixed reference frame. In Eq. (3), the generalized force vector associated with the rigid body translational coordinates is simply the transformation of the bushing force friom the bushing coordinate system to the inertial reference frame, while in Eq. (4), the generalized force vector associated with the rigid body orientation coordinates characterizes the monment of the bushing force vector about the origin of the body-fixed reference frame. Considering the bushing force vector fi, the bushing force is modeled as a set of nonlinear viscoelastic forces which depend not only on the instantaneous bushing deformation, but also on the history of the bushing deformation. The nonlinear viscoelastic bushing forces developed at the Center are based on a modified Pipkin-Rogers superposition principle2s where the nonlinear viscoelastic force is given by itaR[Ai(S),t-S fi(t) = R[Ai(t),O] + J~ a(t-S) where fi is the bushing force vector measured with respect to the bushing coordinate system, Ai is the bushing deformation vector measured with respect to the bushing coordinate system, and R is the relaxation function which characterizes the bushing's viscoelastic response. The first term on the R.H.S. of Eq. (6) represents the instantaneous response of the bushing, while the second term on the R.H.S. of Eq. (6) represents the history-dependent response of the bushing force. The bushing deformation vector Ai whose components are measured with respect to the bushing coordinate system is related to the global bushing deformation vector by the following transformation: Ai= At A!dij (7) where dij is the bushing deformation measured with respect to the inertial reference frame, and the orthogonal matrices Ap and Ai are the bushing coordinate system-to-body reference frame and body reference frame-to-inertial reference frame transformation matrices, respectively. The global bushing deformation vector dij is in turn related to the generalized coordinates through the following relation: di =((Rj+Aj ujp)-(Ri +Aiuj)) (8)

where RI and RJ are the position vectors of the origins of the body-fixed reference frames attached to bodies i andj, respectively, Ai and Aj are the body reference frame-to-inertial reference frame transformation matrices for bodies i and j, respectively, and up and ujp are the position vectors of the bushing attachment points on bodies i and j, measured with respect to their respective body-fixed reference frames. Eqs. (8), (7), (6), (4), and (3) form the algorithmic sequence for the computation of the generalized forces due to bushing forces acting on body i. A similar procedure applies to the computation of the generalized forces due to the bushing forces acting on body j. 3. Uniaxial Nonlinear Viscoelasticity Thus far, we have not made any assumption on the mathematical description of the nonlinear viscoelastic bushing force, except for the fact that the bushing force can be decomposed into an instantaneous force response and a history-dependent force response that is characterized by a convolution integral, the kernel of which is the time derivative of the vector of relaxation functions as depicted in Eq. (6). As a first attempt, we can assume that the bushing force components are decoupled, i.e., radial forces, axial forces, torsional moments, conical moments, and their work-conjugate deformation components are independent. With this assumption, the viscoelastic bushing force can be decomposed into two force components along two orthogonal radial directions and a force component along the axial direction. Similarly, the bushing moments can be decomposed into two conical moments about the two orthogonal radial directions and a torsional moment about the axial direction. Furthermore, the bushing force or moment components depend only on their associated work-conjugate deformation component. Using the above assumption, each component of the nonlinear viscoelastic force-displacement relation of Eq. (6) simplifies to the uniaxial nonlinear viscoelastic force-displacement relation, which can be written as the following scalar equation: f(t) =RJA(t),O]+ R[a(t-s ds (9) where f(t) is the viscoelastic bushing force component along a radial or axial direction, and A(t) is the corresponding bushing deformation along the same direction. The scalar function R[Av(s),t-s] is the relaxation function of the bushing force component along the given direction. We further assume that the relaxation function can be expressed as the sum

of integer powers of the bushing deformation where each term in the series is multiplied by a time-dependent function, i.e., R[A(s),t-s]= Ak() Gk(t-S) (10) k=1 where each of the time-dependent functions Gk(t-s) is expressed as a series of exponential terms, i.e., Nk Gk(t-s) = gkO + E gkj exp[-(t-s)/kj] (11) j=1 Combining Eqs. (9) - (11), we obtain the following expression for the nonlinear viscoelastic bushing force along a particular direction: F(t) = A (t) Gk() gkj I exp[-(t-s)/kj] Ak(s) ds (12) k=l k=l j=l kj o where each of the stiffness coefficients in the instantaneous force response is given by Nk G(O) = gkj (13) j=o The parameters N, Nk, gko, gkj, and tkj (j=l,Nk and k=l,N) are bushing material parameters which are determined from experiments and processed through an extrapolation procedure. The extrapolation procedure and the bushing parameter identification process are described in a separate report. The bushing force component of Eq. (12) is calculated by applying the trapezoidal rule to compute the convolution integral, resulting in the following discrete time approximation of the bushing force: N k N Nk j pi* F(tn+l) = Y A (tn+) Gk(0) E gkij I (14) k=l k=l j=l tkj where Tkj is the numerical approximation to the convolution integral and computed from the following recursive equation:

Ij = — kj exp(-(tt ) + t) exp(-(tn+l-tn)/Tlj) + Ak(tn+l)] (tn+ tn) (15) and the recursion is started by setting kj = O (16) The recursive algorithm for the numerical computation of the nonlinear viscoelastic force obviates the need to store the entire deformation history of the bushing, thus improving the computational efficiency in the calculation of the bushing forces. It is worthwhile to note that the recursive algorithm is possible only because we have chosen to use exponential terms in the time-dependent functions Gk(t-s). The computation of the viscoelastic bushing forces is carried out in ADAMS through a user-written subroutine. Appendix A contains the source code of an ADAMS user-written subroutine that calculates the nonlinear viscoelastic forces for any bushing whose material parameters are stored in a data file. 4. Validation In order to validate the implementation of the nonlinear viscoelastic bushing force element into the multibody dynamics code, we use the code that we have developed to model a simple multibody system consisting only of two bodies, namely, an automotive suspension system component and the ground, wherein the two bodies are connected by bushings. Two bushing models are used for comparison: the conventional nonlinear elastic model and the proposed nonlinear viscoelastic bushing model, hereafter referred to as NLE and VISCO,respectively. The material parameters for both models were obtained from a suite of tests performed for every bushing. For this purpose, a special bushing testing machine was built at the Center. The bushing testing machine has two linear actuators which can provide axial and radial forces along the bushing's principal directions, and two rotary actuators which can provide torsional and conical bushing moments. The bushing testing machine was used in conducting static tests, harmonic tests, ramp tests, and simulated road load tests on each individual bushing in order to determine the timedependent response of the bushing. A parameter identification procedure was employed to process the results obtained from the bushing tests and to come up with a set of bushing

material parameters that can be used directly in Eq. (12). The experimental procedures and the parameter identification procedure are described in a separate report. In order to validate the results of the simulation, an experimental set-up of the automotive suspension system component was also built at the Center. For this purpose, a remote parameter control (RPC) test set-up was built at the Center. The RPC test set-up consists of a lower control arm supported by two bushings and subjected to prescribed dynamic loads and displacements at the ball joint attachment to the spindle. Load cells were attached to the bushings in order to measure the forces transmitted across the bushings. In addition, LVDT's were attached at the outboard end of each bushing to measure the bushing deformation. Dynamic tests were conducted for this RPC test set-up, and the experimental results are compared with the simulation results that are obtained from the VISCO and NLE bushing models. A finite element model of the lower control arm shown in Figure 2, wherein the bushing supports at the front end and at the rear end are illustrated. The ball joint at the lateral end is subjected to prescribed dynamic loads along the fore-aft and lateral directions and a prescribed displacement in the vertical direction. The prescribed forces and displacements can be seen as insets in Figure 2. The entire event takes place in 1.8 seconds. The maximum peak-to-peak applied force along the fore-aft direction is 13.4 kN and the maximum peak-to-peak applied force along the lateral direction is 13.2 kN. The maximum peak-to-peak prescribed displacement along the vertical direction is 46.7 mm. These three time-varying inputs are applied to the RPC test as well as the simulation using the VISCO and NLE bushing models. Figures 3 and 4 show the fore-aft bushing forces for the front bushing and rear bushing, respectively. These figures show the time history of the bushing forces obtained from the RPC test, VISCO simulation, and NLE simulation. The figures show that the maximum difference, in terms of peak-to-peak forces, between the nonlinear viscoelastic bushing model and the nonlinear elastic bushing model is approximately 15% of the maximum peak-to-peak applied force, while the maximum difference between the nonlinear viscoelastic bushing model and the measured test data is approximately 12% of the maximum peak-to-peak applied force. We note that the high frequency response in nonlinear elastic bushing model is damped out because the linear viscous damper in the nonlinear elastic bushing model provides energy dissipation that is linear with frequency. We also note that the RPC test can not measure the high frequency response because of limitations in the sampling rate of the load cells. Table 1 summarizes the comparison of bushing forces obtained from NLE simulation, VISCO simulation, and test results. From

this table we observe that the predicted nonlinear viscoelastic bushing forces obtained from the VISCO model are closer to the test results than those obtained from the NLE model. Figures 5 shows the rear bushing displacement along the fore-aft direction. This figure shows the time history of the bushing displacement obtained from the RPC test, VISCO simulation, and NLE simulation. As expected, the nonlinear viscoelastic bushing model shows higher peaks in the predicted bushing displacement compared to the nonlinear elastic bushing model. Again this is attributed to the fact that the conventional nonlinear elastic model damps out the high frequency response. This figure also shows that the predicted peak displacements and the measured peak displacements do not compare as well as the corresponding bushing forces. It is interesting to note that the nonlinear viscoelastic bushing model captures the frequency of the measured bushing displacement whereas the nonlinear elastic bushing model shows errors in amplitude and a phase shift in the bushing displacement with respect to the measured displacements. Table 2 summarizes the comparison of the maximum peak-to-peak displacement response obtained from NLE simulation, VISCO simulation, and test results. The table shows that the nonlinear viscoelastic bushing model gives a better estimate of the actual bushing displacement than the nonlinear elastic bushing model. Also from this table, we observe that the predicted bushing displacements obtained from the multibody system dynamics simulation do not agree well with the measured displacements. The discrepancy between the predicted displacements and measured displacements can be as high as 90% of the measured displacements. This discrepancy is due to the fact that component flexibility has been ignored in the formulation of the equations of motion. This contention is supported by calculations performed using the fully nonlinear finite element code ABAQUS22 wherein the flexibility of the control arm is considered. A comparison of dynamic responses between a rigid body model and a flexible body model of the lower control arm is shown in the time history plots of Figures 6 through 7. The finite element calculations show that there is load redistribution due to the flexibility of the component, so that the front bushing experiences larger peak-to-peak displacements and the rear bushing experiences smaller peak-to-peak displacements, which are consistent with the experimental results. The results of the finite element analysis also indicated that the predicted bushing forces can be improved if component flexibility is taken into account. Figures 8 and 9 clearly show the load redistribution due to component flexibility. These figures also show the high frequency response due to the vibration of the flexible component. It should also be noted that the finite element calculations, which required a fully nonlinear dynamic analysis due to the relatively large vertical motion imposed at the ball joint, required computation times

which are orders of magnitude higher compared to that of the multibody system dynamic analysis. 5. Extension to Flexible Multibody System Dynamics Based on the results of the validation of the code as reported in the previous section, the assumption of rigid bodies in multibody dynamics is adequate if bushing forces are the only concern in the dynamic simulation. However, if the displacements of some points of interest are required from the dynamic simulation, the assumption of rigid bodies is no longer adequate even for relatively stiff components such as the lower control arm, as the results of the code validation have shown. In this regard, we present an extension of the formulation of bushing forces to the dynamics of flexible multibody systems. Similar to the formulation of the generalized forces due to bushing forces in the equations of motion for rigid multibody systems, we first present the equations of motion for a flexible multibody system, and subsequently introduce the generalized forces due to the bushing forces that act on a deformable body in the multibody system. Consider now the two bodies i and j in Figure 1 to be deformable bodies. The equations of motion for the flexible multibody system can be written in the following descriptor form23: [ mRR RO mRf QR 1 OR FR moR moo mf + Qo + 0+ + FO (17) mfR mfR mff J j t Qf J PfJ LFf Lqf.f IL subject to the holonomic constraint equations *(R,O, qf)=0 (18) where the generalized coordinates (R,9, qf) refer to the translation of the body reference frame, the orientation of the body reference frame, and the deformation of the flexible body, respectively. The mass matrix is a function of the orientation of the body reference frame and the elastic deformation. The force vector Q refers to the generalized forces due to externally applied loads, and the force vector F includes the centrifugal forces as well as the Coriolis forces arising from the elastic deformation. The equations of motion as expressed in Eq. (17) are based on the assumption that small elastic deformations are superposed on the gross rigid body motion. The force vector Pf refers to the elastic forces within the deformable body. This force vector can be decomposed into a force vector due to the

(linear) structural stiffness matrix and a force vector due to a stress-dependent stiffness matrix which captures the stiffening effects from the rigid body rotation. Similar to the case of rigid multibody systems, the bushing forces enter the equations of motion through the generalized force vector Q by treating the bushing forces as external forces that act between two bodies or parts that are connected by a bushing. The generalized force vector due to the bushing forces will be of the same form as that of the generalized bushing force vectors formulated for the case of rigid multibody systems, with some minor modifications which account for the elastic deformation of the bodies which are connected by the bushing. Again, in the following discussion, we will consider only the contribution of radial and axial bushing forces to the generalized bushing force vector. The contribution of conical and torsional bushing moments to the generalized bushing force vector can be easily derived by applying the principle of virtual work. If we consider a bushing force vector fi acting on deformable body i and whose components are measured with respect to the bushing coordinate system ri - si ti at the bushing attachment point Pi, the generalized bushing force vector associated with the translation of the body reference frame of body i is given by QR = Ai[I + ] Apfi (19) where Ap is the rotation transformation matrix from the bushing coordinate system to the body reference frame in the undeformed configuration, fp is the skew-symmetric matrix whose associated axial vector is the elastic rotation vector at the bushing attachment point. Again, Eq. (19) is valid only for small elastic rotations. The generalized bushing force vector associated with the rotation of the body reference frame of body i is given by Q= G ib[I+ p] Apfi (20) where ub is the instantaneous position vector of the bushing attachment point measured with respect to the body reference frame. Unlike the case of rigid bodies, uj is not constant since it now depends on the elastic deformation of the deformable body. When a finite element discretization or an assumed mode formulation is used for the spatial interpolation of the flexible body's deformation, position vector of the bushing attachment point can be expressed as up =(u)r + Np q (21)

where the first term (u)r is the position vector of the bushing attachment point in the undeformed configuration and second term is the elastic deformation of the bushing attachment point, which is defined in terms of the interpolation or shape function matrix Np evaluated at point Pi, and the vector of generalized elastic coordinates qif. Finally, the generalized bushing force vector associated with the deformation coordinates of body i is given by Q= N [I + p] Apfi (22) The nonlinear viscoelastic bushing force vector fi is expressed as a function of the bushing deformation history fi(t) = Ai[(t),0] + J a (t-s) ](23) where the local bushing deformation vector is related to the global bushing deformation vector through the transformation relation Ai= A[I +p] At dij (24) and the global bushing deformation vector is determined from the generalized coordinates through the following equation: dij= (Ri + Aj +Nj [(u + Np qf])- (Ri + Ai [( + N qi)) (25) Eqs. (25), (24), (23), (19), (20) and (22) form the algorithmic sequence for the computation of the generalized force vectors due to bushing forces acting on body i. A similar procedure is needed to determine the generalized force vectors due to bushing forces acting on body j. 6. Conclusion In this report, we have presented the formulation of nonlinear viscoelastic bushing forces as massless force elements between two bodies in a multibody system. The

numerical implementation of the proposed formulation into the general-purpose multibody dynamics code ADAMS has been completed for rigid multibody systems. The validation of the resulting rigid multibody dynamics code has been performed by comparing results obtained from a nonlinear viscoelastic bushing model and a nonlinear elastic bushing model to those obtained from measured test results. The code validation and comparison with measured data revealed that the nonlinear viscoelastic bushing model gives a more accurate prediction of dynamic loads and displacements than the nonlinear elastic bushing model. However, the rigid multibody system models did not provide accurate predictions of the displacements of specific points of interest. The cause of the poor displacement prediction capabilities of the rigid multibody dynamics models was traced to the fact that component flexibility played an important role in the prediction of displacements, and also to a lesser degree, in the prediction of dynamic loads. In order to achieve better prediction capabilities in the displacements and dynamic loads, the formulation of nonlinear viscoelastic bushing forces has been extended to include deformable bodies connected by bushings in the context of general flexible multibody systems. The proposed formulation of nonlinear viscoelastic force elements in flexible multibody systems will be implemented in a future release of ADAMS which includes the capability for modeling flexible components and user-subroutines for interfacing user-defined forces that act on the flexible components. References: 1. J. Unda, J. Garcia de Jalon, F. Losantos and R. Esperantza,'A comparative study on some different formulations of the dynamic equations of constrained mechanical systems', ASME Journal of Mechanisms, Transmissions and Automation in Design, vol. 109, 466474, 1987. 2. E. Bayo and A. Avello,'Singularity-free augmented Lagrangian algorithms for constrained multibody dynamics', Nonlinear Dynamics, vol. 5, 209-231, 1994. 3. 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. 4. K. Park and J. Chiou,'Discrete momentum-conserving explicit algorithm for rigid body dynamic analysis', International Journal for Numerical Methods in Engineering, vol. 36, 1071-1083, 1993. 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 Journalfor 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. U. Holzlohner,'A finite element analysis for time-dependent problems', International Journalfor Numerical Methods in Engineering, vol. 8, 55-69, 1974. 12. Y. Yamada, H. Takabatake and T. Sato,'Effect of time-dependent material properties on dynamic response', International Journalfor Numerical Methods in Engineering, vol. 8, 403-414, 1974. 13. C. Johnson and D. Kienholz,'Finite element prediction of damping in structures with constrained viscoelastically damped structures', AIAA Journal, vol. 20, 1284-1290, 1982. 14. K. Morman and J. Nagtegaal,'Finite element analysis of sinusoidal small-amplitude vibrations in deformed viscoelastic solids. part 1: theoretical development', International Journalfor Numerical Methods in Engineering, vol. 19, 1079-1103, 1983. 15. R. Bagley and P. Torvik,'Fractional calculus in the transient analysis of viscoelastically damped structures', AIAA Journal, vol. 23, 918-925, 1985. 16. D. Golla and P. Hughes,'Dynamics of viscoelastic structures - a time domain finite element formulation', ASME Journal of Applied Mechanics, vol. 52, 897-906, 1985. 17. K. Xie, L. Roystre and R. Ciskowski,'A boundary element method formulation for fractional operator modeled viscoelastodynamic structures', Advances in Boundary Elements: Vol. 3, Stress Analysis, 55-64, Computational Mechanics Publications, 1989. 18. L. Gaul and C. Chen,'Modeling of viscoelastic elastomer mounts in multibody systems', Advanced Multibodry System Dynamics, 257-276, Kluwer Academic Publishers, 1993.

19. S. Yi and H. Hilton,'Dynamic finite element analysis of viscoelastic composite plates in the time domain', International Journalfor Numerical Methods in Engineering, vol. 37, 4081-4096, 1994. 20. ADAMS (Automatic Dynamic Analysis of Mechanical Systems) Reference Manual, Mechanical Dynamics, Inc., Ann Arbor, Michigan. 21. A. Pipkin and T. Rogers,'A non-linear integral representation for viscoelastic behavior', Journal of Mechanics and Physics of Solids, vol. 16, 59-72, 1968. 22. ABAQUS User's Manual, Hibbitt, Karlsson & Sorensen, Inc., Pawtucket, Rhode Island. 23. A. Shabana, Dynamics of Multibody Systems, Wiley, 1989.

Appendix A: ADAMS User-Subroutine for Nonlinear Viscoelastic Force Elements subroutine gfosub(id,time,par,npar,dflag,iflag,result) c c external variable definitions c integer id double precision time double precision par(*) integer npar logical dflag logical iflag double precision result(6) c c local variable definitions c integer brflag brflag=nint(par( 1)) if(brflag.eq. 1) then call visco l(id,time,par,npar,dflag,iflag,result) else call errmes(.true.,' unrecognized branch in gfosub',id,'STOP') endif return end c subroutine for implementing uncoupled nonlinear viscoelasticity c with up to 5 exponential terms for each of G1,G2,G3,G4 and G5. c subroutine viscol(id,time,par,npar,dflag,iflag,result) implicit double precision (a-h,o-z) c external variable definitions c integer id double precision time double precision par(*) integer npar logical dflag logical iflag double precision result(6) c c local variable definitions c integer bushid integer idstrn integer markl integer mark2 integer nsize integer nstates integer ipar(3) integer i,j,k integer mode integer nchars integer istat

double precision u(6),disp(3) double precision atime,sig double precision g(5),s(5) double precision fn,fnpl,dk,h logical errflg character*80 string,dummy integer n(5,6,100) double precision timen( 100),timenp 1(100) double precision lan(6, 100),lanpl (6,100) double precision kn(5,5,6,100) double precision knp 1(5,5,6, 100) double precision c0(5,6, 100),c(5,5,6, 1 00),f(5,5,6,100) save timen save timenp save lan save lanpl save kn save knpl save cO save c save f save n c assign understandable names to passed parameters bushid=nint(par(2)) markl=nint(par(3)) mark2=nint(par(4)) idstrn=nint(par(5)) c markl=marker attached to outer cylinder of bushing c mark2=marker attached to inner cylinder of bushing c —---------------------— initializationc if (iflag) then call gtstrg(idstrn,string,nchars,istat) open(unit=96,file=string,status='old') do k=1,6 read(96,*) dummy read(96,*) (cO(i,k,bushid),i=1,5) read(96,*) dummy do j=1,5 read(96,*) (c(i,j,k,bushid),i=1,5) end do read(96,*) dummy do j=1,5 read(96,*) (f(i,j,k,bushid),i=l,5) end do end do close(96) do k=1,6 do i=1,5 n(i,k,bushid)=O do j=1,5 if(dabs(c(i,j,k,bushid)).gt. 1.OD-10) n(i,k,bushid)=j end do

end do end do endif c —--------------— end initialization —---------- c find the relative displacement between the 2 markers nsize=3 ipar(1)=mark2 ipar(2)=markl ipar(3)=markl call sysary('TDISP',ipar,nsize,disp,nstates,errflg) call errmes(errflg,'error calling sysary for TDISP',id,'STOP') do i=1,3 u(i)=disp(i) end do c find the relative rotation between the 2 markers nsize=2 ipar(l)=mark2 ipar(2)=markl call sysfnc('AX',ipar,nsize,u(4),errflg) call errmes(errflg,'error calling sysfnc for AX',id,'STOP') call sysfnc('AY',ipar,nsize,u(5),errflg) call errmes(errflg,'error calling sysfnc for AY',id,'STOP') call sysfnc('AZ',ipar,nsize,u(6),errflg) call errmes(errflg,'error calling sysfnc for AZ',id,'STOP') call getmod(mode) if(mode.eq.5) then c static analysis (fully relaxed state) do k=1,6 lanpl(k,bushid)=u(k) end do timen(bushid)=O.OdO timenpl(bushid)=O.OdO do k=1,6 do i=1,5 g(i)=cO(i,k,bushid) end do do i=1,5 s(i)=g(i)*lanpl(k,bushid)**i end do sig=0.0d0 do i=1,5 sig=sig+s(i)

end do result(k)=sig c initialize static variables do i=1,5 do j=l,n(i,k,bushid) kn(i,j,k,bushid)=O.OdO knpl(i,j,k,bushid)=O.OdO end do end do lan(k,bushid)=lanpl(k,bushid) end do else c dynamic analysis call timget(atime) if(atime.lt. 1.OE-14) then do k=1,6 lanpl(k,bushid)=u(k) end do timenp 1 (bushid)=time do k=1,6 do i=1,5 g(i)=cO(i,k,bushid) if(n(i,k,bushid).gt.0) then do j=l,n(i,k,bushid) g(i)=g(i)+c(i,j,k,bushid) end do endif end do do i=1,5 s(i)=g(i)*lanp 1 (k,bushid)* *i if(n(i,k,bushid).gt.O) then do j=l,n(i,k,bushid) fn =exp(-timenpl(bushid)*f(ij,k,bushid))*lan(k,bushid)**i fnp 1 =lanpl(k,bushid)**i dk = O.5dO*timenpl(bushid)*(fn+fnpl) knpl(i,j,k,bushid)=dk s(i)=s(i)-c(i,j,k,bushid)*f(i,j,k,bushid)* * knpl(i,j,k,bushid) end do endif end do sig=O.OdO do i=1,5 sig=sig+s(i) end do result(k)=sig c store static variables timen(bushid)=time

do i=1,5 if(n(i,k,bushid).gt.O) then do j=l,n(i,k,bushid) kn(i,j,k,bushid)=knp 1 (i,j,k,bushid) end do endif end do lan(k,bushid)=lanpl(k,bushid) end do else c determine if previous trial converged if(atime.ge.timen(bushid)) then c previous trial converged, so update static variables timen(bushid)=timenp 1 (bushid) do k=1,6 do i=1,5 if(n(i,k,bushid).gt.O) then do j=l,n(i,k,bushid) kn(i,j,k,bushid)=knp1 (i,j,k,bushid) end do endif end do lan(k,bushid)=lanp 1 (k,bushid) end do endif c determine viscoelastic force for current trial do k=1,6 lanpl(k,bushid)=u(k) end do timenp 1 (bushid)=time h=timenpl(bushid) - timen(bushid) do k=1,6 do i=1,5 g(i)=cO(i,k,bushid) if(n(i,k,bushid).gt.O) then do j=l,n(i,k,bushid) g(i)=g(i)+c(i,j,k,bushid) end do endif end do do i=1,5 s(i)=g(i)*lanp 1 (k,bushid)* *i if(n(i,k,bushid).gt.O) then do j=l,n(i,k,bushid) fn =exp(-h*f(i,j,k,bushid))*lan(k,bushid)* *i fnp 1 =lanp 1 (k,bushid)**i dk = O.5dO*h*(fn + fnpl) knp 1 (i,j,k,bushid )=ex p(-h*f(i,j,k,bushid))*kn(i,j,kbushid)+dk s(i)=s(i)-c(i,j,k,bushid)*f(i,j,k,bushid) * knpl(i,j,k,bushid) end do

endif end do sig=O.OdO do i=1,5 sig=sig+s(i) end do result(k)=sig end do endif endif return end

Maximum Peak-to-Peak Bushing Forces Force Component ADAMS-NLE ADAMS-VISCO Test Results Front Bushing 8.396 kN 10.349 kN 9.503 kN Fore-Aft Direction Front Bushing 10.807 kN 12.399 kN 10.969 kN Lateral Direction Rear Bushing 6.384 kN 5.753 kN 4.056 kN Fore-Aft Direction Rear Bushing 12.077 kN 12.408 kN 12.600 kN Lateral Direction Table 1: Comparison between simulation and test results: maximum peak-to-peak bushing forces Maximum Peak-to-Peak Displacements Displacement ADAMS-NLE ADAMS-VISCO Test Results tOTJ20nent Component Front Bushing Disp. 6.796 mm 8.705 mm 11.762 mm Fore-Aft Direction Front Bushing Disp. 8.104 mm 10.439 mm 10.716 mm Lateral Direction Rear Bushing Disp. 7.460 mm 9.522 mm 5.041 mm Fore-Aft Direction Rear Bushing Disp. 0.790 mm 1.296 mm 1.861 mm Lateral Direction Table 2: Comparison between simulation and test results: maximum peak-to-peak displacements

List of Figures Fig. 1. Mulibody system and reference frames. Fig. 2. Lower control arm and prescribed forces and displacements at the ball joint. Fig. 3. Front bushing fore-aft force: comparison between test results, nonlinear viscoelastic bushing model, and nonlinear elastic bushing model. Fig. 4. Rear bushing fore-aft force: comparison between test results, nonlinear viscoelastic bushing model, and nonlinear elastic bushing model. Fig. 5. Rear bushing fore-aft displacement: comparison between test results, nonlinear viscoelastic bushing model, and nonlinear elastic bushing model. Fig. 6. Front bushing fore-aft displacement: comparison between flexible lower control arm and rigid lower control arm. Fig. 7. Rear bushing fore-aft displacement: comparison between flexible lower control arm and rigid lower control arm. Fig. 8. Front bushing fore-aft force: comparison between flexible lower control arm and rigid lower control arm. Fig. 9. Rear bushing fore-aft force: comparison between flexible lower control arm and rigid lower control arm.

N /~ \ \~ ~ ~~~~~~~~~~~~~~~ / owe ai \5 / 0. cr tax \ C oa N mm* P-. tC~~~~~~~~~~~1

~~~~~~~~~~~~~~~~~......:. —-.l...............l_-..i.....I I i - i...;._ ~"~! —: I -:...'dii N,~~~I3*gg.~.~ —. — ~~~-i ~-.. "'9 1: r....... ~i.....~~.. o ~ ~~~~~~~~ ~~~~~~~~~~~~.ll..11,~.At p~llI I.... F~.~ ~tg; gNUIII * t ~ t L L rt~~~~~~~~~] ~ ill:: DP~~ __________, I )'r ~~~~~~~~~~~~~~~~~~~~........ L........~..:............,...:~lt - -I __ __I 1~ ~~~~~_, ii;lll,,,,,I /I i'i"-,"~~~~~~~~~~_l::,:,-;,r-,!ls~~~l~i.... 4: i~~~~'i ~Il'

1~~ F-o Lfl 8 C' 0~~~~~~~~~~ I~~~~~~~~~' Q LO L.. 0)? o.........................................c.,. LL~~~~~~~~~~~~~~~~~~~~L ~~~~~~~~~~~~~~~~~~~~~~r~ 0~~~~~~~~~~~~~~~~~~~~(.................~~~~~~r.. d I ~ ~ ~ ~~~~~ j__ (0 5! (N~I) eaJo:: I, ( ~r~~~~~~~~~c r ~ ~ ~ ~ ~ ~ L........................................................... (

— qo ~' —'' LO~~~\... ~.c~..............., ~t~II. 1 I * I I~~~~~~ i~tI i~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ i~~~ I i ~ ~ ~~:,' ~-_~,r Z~~' I Il I~~~~~~~~~ --- ~ ~ ~ ~ ~ ~ ~.~"-.In.~~~~~~n-...-.-.-.-.~~~~~~~~~~~...~......-.-.-.-.-. r~~~~~~~. O~'~~ ~ ~ ~ ~~~~~~~~~~,1::~~~~~~~~~ d C, ___~.,:, — ~~r ~:...'-r.~~~~~~~~~~~~~b ~ r~~~~~~ rrrrrrr -5-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~_ I V ~~rrrrrrrrrrrrrrrrrrrrrrrrr~~~rrr~rr,!: t.o.....o.....,... c~ ~ ~~~~~c ~ I- ~~~~~~~~~~~~~~~u r~~~~~~ I~~~~~~~ ~ Oi i i t ~. c ~ ~ ~~~~~~~~~~~~~~~~~~~~:!: 03 ~ I - - - ~~~~~~~~~~ 1 IIL L- ~ ~ II LL~~~~~~~~~NI eo~ ~o_-!

I8''','':'' I';- - ---- - ---------------- ------------- -.. ~. —:.... _..~"',.cry. c/................ I...........7..................'..Y.................... -----—............I............... 4~~~~~~~~~~:::,-..,Y!::::,',L-f —- - I-~ —--- -- 2 L O - -,'.!:: m. --- -- —,.............................................. co X~~~~~~~~~~~~~~~~~~~~~~~~ —---.....,...................:::: s::::::3........,~.......................................................(.................................... I::::;:::: ~,-:, -:::.: Z!! i! _ i"!!:0: -.J.;.,....:::::: O~~~~~~~~~~~~~~~~~~~~L II~~~~......r..... _...................................... co@ @ 0 0 \ e @ @ -:D::':O:d, r - C::':J:' (ww) ]uawooBldsln

"'I:,: v /.....,...... rr'.:.:.::, i~~~~ 0I OU 1 1 ~.(.::::::. m * ---------— ) —-- ----- ------—:'1 I'''' I *, _, ~. Q.....:: 00 11 i i i i: -'"~''~ A..............":..............i........ -.....-.-....-.....-..... —.-...............-....... ~~~~~~~~~~o** Cl)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~f o............ ~::.~..........=...... _~~~~~~~~~~~~.............-.....................................:........~~~~~~~~~~~~~~~~~~~~~~~~~~.. /........ 0.,. ~~~~~~~~~~0.,, ~.::::::: (0 0 q..,f C} ~ t _- ^,, ls, 5,,, _ ~~~~~~~- C D.CL~~~~~~~~,CI) ~~~~~~~~~~~ II. I (wwL) leuawe~!ds!a... t.., \.J.-.~~~~~~~~~ \~~~~~L~~l~~~~~C L''' )''' @' I (9~~~~~~6 ClClO N OdoU ~~~.~.~.~.~.~. ~.~~~.~.~~~..1,.~~.~.~~~~-~-~I.~~~.~.~~.~.~.~~~~~.~.~.~..~.~~~.~.~~.~.~~~.~~~...~~~.~~~. 00 (Uln]) la~asldSd

........................................ C~ C...........................................i.......................... coo ml~(.'e.,_,:. ). 1(:w s 1t I.COn t _................................................................................................... ed I —----------— 1 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ c

I ~ ~,, I I.:1 x:::::: f':: az t,...,r LLI''.'''', r..::::;;: i i _ _ *~~~~~~~~~~~~~~~,.,,,.,,............~~~~~~~~~~~~~ ~~~~-~- — i —--—,-i —~.,.,...-...,...,,.,...........-~,..,..........,.~ 0 <CO1. =.~~~~~~~~~~~~~~~~~~ i.. i F 1L~~~~~~~cni I ~ ~ o LL,,..' C 3 ~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ I~ I....o..,,,,(, I.::. _...................................................................................._......I....I....... #.I....,.,,,., d...... _.......,,_c...... ~ —-- ~ — ~~~,~~~~~-~-~~,.- o.. ~~~......,~.~~ ~ ~ ~....~...~....)......rr rr rr~r rr rr rr r.~~~~~~~.... L............... f.~ i. i ii { I, r~~~~~~~ N ( U)t XCM O M C Q)~~~~~~(M O3o

II I I.. A.},. ~7.............. o. 1 0 LL 0) C /Jt...-............................................................................................ L I..., ILO [.,.. 0 (I0 LO) O ~ ~ ~ ~ ~ ~~~NI eo~o-! —---- ~~~ — ~-~-~ —1 —--— ~ —~ —-------