INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 44, 439–459 (1999) A GENERAL CONJUGATE-GRADIENT-BASED PREDICTOR–CORRECTOR SOLVER FOR EXPLICIT FINITE-ELEMENT CONTACT EDWARD ZYWICZ ∗ AND MICHAEL A. PUSO The University of California, Lawrence Livermore National Laboratory, P.O. Box 808, Livermore, CA 94550, U.S.A. SUMMARY A Lagrange-multiplier based approach is presented for the general solution of multi-body contact within an explicit ÿnite element framework. The technique employs an explicit predictor step to permit the detection of interpenetration and then utilizes a corrector step, whose solution is obtained with a pre-conditioned matrixfree conjugate gradient projection method, to determine the Lagrange multipliers necessary to eliminate the predicted penetration. The predictor–corrector algorithm is developed for deformable bodies based upon the central dierence method, and for rigid bodies from momentum and energy conserving approaches. Both frictionless and Coulomb-based frictional contact idealizations are addressed. The technique imposes no timestep constraints and quickly mitigates velocity discontinuities across closed interfaces. Special attention is directed toward contact between rigid bodies. Algorithmic moment arms conserve the translational and angular momentums of the system in the absence of external loads. Elastic collisions are captured with a two-phase predictor–corrector approach and a geometrically approximate velocity jump criterion. The ÿrst step solves the inelastic contact problem and identiÿes inactive constraints between rigid bodies, while the second step generates the necessary velocity jump condition on the active constraints. The velocity criterion is shown to algorithmically preserve the system kinetic energy for two unconstrained rigid bodies. Copyright ? 1999 John Wiley & Sons, Ltd. This paper was produced under the auspices of the U.S. Government and it is therefore not subject to copyright in the U.S. KEY WORDS: contact; Lagrangian multiplier; ÿnite element 1. INTRODUCTION In large-scale explicit ÿnite element simulations, e.g. vehicle crashworthiness analyses, the computational cost of contact typically exceeds all other costs combined and is dominated mainly by the search algorithm. Yet, when contact is detected, restoration forces are commonly generated by the computationally inexpensive penalty method1 or approximate Lagrange-multiplier forms2 that impose global time-step limitations, may require persistent penetration, and even suer pathological deÿciencies. An alternative approach to these methods are Lagrange-multiplier-based techniques ∗ Correspondence to: Edward Zywicz, Structural=Applied Mechanics Group, Lawrence Livermore National Laboratory, University of California, P.O. Box 808, Livermore, CA 94550, U.S.A. Contract=grant sponsor: U.S. Department of Energy; Contract=grant number: W-7405-Eng-48 This paper was produced under the auspices of the U.S. Government and it is therefore not subject to copyright in the U.S. CCC 0029–5981/99/040439–21$17.50 Copyright ? 1999 John Wiley & Sons, Ltd. Received 12 September 1997 Revised 21 April 1998 440 E. ZYWICZ AND M. A. PUSO that determine the restoration forces via a predictor–corrector formulation.3 – 5 These techniques predict the nodal positions at the next time step and then, unlike penalty methods or approximate Lagrange-multiplier forms, generate corrective contact forces that eliminate all predicted penetration. Lagrange-multiplier methods are an attractive alternative from a mechanical standpoint because they obviate penalty sensitivity, eliminate persistent penetration, and improve the resolution of contact forces. On the other hand, they are slightly more expensive computationally. Nevertheless, this extra cost can be oset by incorporating simpler, more robust and cheaper search algorithms that exploit the method’s ability to eliminate penetration. In the explicit forward increment Lagrange-multiplier method,5 an initial predictor step at each time step is made in which no contact constraints are employed. In the predicted conguration, interpenetrations between candidate node and segment pairs are identied and used to form the contact matrix. In the corrector step, the discrete equations of motion involved with contact are re-solved with Lagrange multipliers included to eliminate the predicted penetrations. Along with eliminating penetrations, it is also important that the solution satisfy the Kuhn Tucker complementary conditions, i.e. no tensile contact forces or forces on open gaps are allowed, otherwise the system may be overconstrained. The corrector step involves the solution of a coupled system of equations whose dimension equals the current number of Lagrange multipliers, and it is the solution of these equations and constraints that represent the largest computational cost in this approach. For a variety of reasons, e.g. computational speed, non-linear contributions, linearly dependent equations, it is preferable to use iterative techniques to solve the system of equations.4 Recent work has demonstrated the capability of the explicit Lagrange multiplier method, with iterative solution strategies, in both two-dimensional (2D) and three-dimensional (3D) implementations.3 – 6 In Reference 5 a modied, matrix-free Gauss–Seidel iterative method is used to solve the system of equations. During the iterations, tensile contact forces are prohibited by setting the force to zero whenever a tensile value is calculated; however, no explicit enforcement of zero forces on open gaps is made. In References 3 and 6, it is shown that the discretized equations of motion can be reduced to a linear complementary problem, the solution of which is equivalent to a quadratic programming problem. In their work, a Fletcher Reeves conjugate gradient (CG) projection method is developed to solve the linear complementary equations. This algorithm precludes both tensile contact forces and forces on open gaps. Because the CG iterations are performed with constant global matrices, non-linear eects are not addressed. This work focuses on developing an ecient Lagrange multiplier implementation for general 3D utilization. It extends the method to include the non-linear kinematics of rigid bodies and, therefore, permits exible-to-exible, exible-to-rigid, and rigid-to-rigid contacts. Because of its superior convergence characteristics, a CG solution strategy is adopted to solve the linear complementary problem. In an eort to make the CG projection algorithm more general and computationally ecient, a matrix-free version is developed which includes a preconditioner. In contrast to the aforementioned works, friction forces are backcalculated from the Lagrange multipliers associated with the normal contact forces by using a penalty approximation. While this can potentially compromise the modelling of the stick-slip condition, it reduces the number of Lagrange multipliers for a 3D implementation by one-third (References 3 and 6) or one-half (Reference 5). Also, unlike Carpenter et al.,5 the contact matrix is not updated during the iteration procedure. While the update is a trivial calculation in 2D, it is considered too expensive in the present 3D implementation. The importance of this update is unclear, especially since the contact matrix denition is not necessarily unique. Finally, the entire algorithm is fully implemented in the explicit 3D nite element code DYNA3D.7 Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) EXPLICIT FINITE-ELEMENT CONTACT 441 In DYNA3D, the equations of motion for rigid bodies are integrated implicitly using energyand momentum-conserving algorithms similar to those presented by Simo and coworkers.8; 9 The implicit treatment is done so that the rigid materials do not control the time step for the entire problem. The work of Simo is extended here to include contact such that algorithmic conservation of energy and momentum is maintained for a system of contacting rigid bodies free of external forces. To accommodate the resulting rigid-body integration algorithms, as well as the abovementioned treatment of frictional contact, a nonlinear version of the linear CG projection method3 is employed. Many practical aspects of the contact solver are also discussed. As just noted, contact between rigid bodies is rigorously considered. Special attention is paid in the implementation of the iterative solver in order to minimize computational cost. Because many commonly used search algorithms employ dual pass approaches and generate multiple, linearly dependent, contact equations, the practical treatment of redundant equations is addressed and demonstrated. Also, the consequences of enforcing a previously missed contact are quantied. The paper proceeds as follows. The governing equations for the forward increment Lagrangemultiplier method are developed for deformable bodies based upon the central dierence timeintegration scheme and for rigid bodies based upon a momentum and energy-conserving algorithm. Special attention is directed to algorithmic momentum and energy conservation for contact between rigid bodies. The kinematic behaviour of the gap is explored, and the criterion needed to model elastic collisions between rigid bodies is developed. Next, the pre-conditioned, iterative solver, and matrix-free update solution strategy are presented. Finally, several academic and practical examples show the eectiveness and overall computational cost of the approach. 2. GOVERNING EQUATIONS FOR DEFORMABLE BODIES 2.1. Equations of motion In the forward increment Lagrange-multiplier method,5 the equations of motion are expressed in their semi-discrete form at time tn , with the Lagrange multipliers associated with contact explicitly segregated, as Mx n + fnInt (xn ) + Gn+1 [n + dn+1 ([n ; xn+1 ) = fnExt (1) Here M is the mass matrix, G the contact matrix in the normal direction,5; 6 [ the Lagrange multiplier (normal contact force) vector, d the frictional contact force vector operator, f Int the stress divergence or internal force vector, f Ext the external force vector, and x the nodal position vector The notation (˙) denotes dierentiation of ( ) with respect to time, and repeated indices do not imply summation. While (1) can be integrated temporally with various one-step methods, attention within this paper is focused solely upon the explicit central dierence method. For DOF (degrees of freedom) associated entirely with deformable bodies, application of the central dierence method to (1) yields the well-known update relationships fn = fnExt − fnInt (xn ) − Gn+1 [n − dn+1 ([n ; xn+1 ) ẋn+1=2 = ẋn−1=2 + t˜n M−1 fn xn+1 = xn + tn ẋn+1=2 (2) (3) (4) where d(0; x) = 0, tn = tn+1 − tn and t˜n = (tn−1 + tn )=2 Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 442 E. ZYWICZ AND M. A. PUSO 2.2. Gap function The distance between a candidate point or node and the potential contact surface is quantied via the gap function, dened as g = nT (Xs − Xm ) − (5) s Here X is the position of the candidate or ‘slave’ node, is an oset (e.g. used to account for shell element thickness), and n is the unit outward normal of the contact surface at the ‘master’ or predicted contact point Xm . (Henceforth, slave and master quantities are denoted by the superscripts s and m, respectively.) In general, the point Xm , obtained from a closest point projection of the slave node onto the master segment, does not coincide with a node, and is interpolated from nodal values associated with the contact segment. The gap vector, g, is dened in terms of x and G as g = GT x − T (6) To prevent interpenetration, ensure compressive normal tractions at the contact points, and create restoring forces only during active contact3 g¿0; [60 and gT [ = 0 (7) 2.3. Friction Translational motion between slave and master points is restricted with an approximate, historydependent slack-type Coulomb friction model.10 The frictional force vector, e, between a slave and master contact pair is given in its discrete form by en = −mn min(−n (kuk=tn ); kẽn k) (8) where mn = ẽn =kẽn k (9) ẽn = (en−1 − u) − n((en−1 − u)T n) m u = Xm n (tn ) − Xn−1 (tn ) (10) (11) (kuk=tn ) is the velocity-dependent friction coecient function, and is the slack stiness. m Prior to contacting the current master surface, e = 0. (In general, Xm n and Xn−1 do not correspond to the same material point.) The vector operator d([; x) represents the scattered nodal contributions of the individual frictional forces. When x is held constant, d([) demonstrates a pseudo-linear dependence upon [. In the case when i i ¿ − kẽi k for all i, d depends linearly upon [, i.e. d = F[; where F is a constant matrix. Conversely, when i i ¡ − kẽi k for all i, d is a constant vector independent of [. This pseudo-linear behaviour is exploited by the predictor–corrector solution procedure. 2.4. Predictor–Corrector Following Carpenter et al.5 the predictor–corrector relationships are assembled. The predicted conguration, denoted with barred quantities, is constructed at tn+1 by assuming [n = 0. It is given by (12) ẋ n+1=2 = ẋn−1=2 + t˜n M−1 [fnExt − fnInt (xn )] x n+1 = xn + tn ẋ n+1=2 Copyright ? 1999 John Wiley & Sons, Ltd. (13) Int. J. Numer. Meth. Engng. 44, 439–459 (1999) EXPLICIT FINITE-ELEMENT CONTACT 443 Contact, potential contact, and interpenetration are identied using ẋ n+1=2 , x n+1 , and traditional global and local search algorithms. Based upon this information, the operators G n+1 and dn+1 are assembled. To generate a pseudo-explicit procedure for deformable bodies, G n+1 and dn+1 are used in (1) in place of the unknown quantities Gn+1 and dn+1 , respectively. Equation (4) is now expressible entirely in terms of the known quantities and the unknown Lagrange multipliers as xn+1 = x n+1 − tn t˜n M−1 [G n+1 [n + dn+1 ([n ; x n+1 )] (14) Substitution of (14) into (6) yields the relationship gn+1 = g n+1 − Bn [n − bn ([n ) (15) where g n+1 = G Tn+1 x n+1 − T Bn = tn t˜n G Tn+1 M−1 G n+1 bn = tn t˜n G Tn+1 M−1 dn+1 ([n ; x n+1 ) (16) (17) (18) Thus, (15) needs only be solved subject to the constraints imposed by (7). The current three-dimensional predictor–corrector algorithm deviates from the others3 – 6 in several ways. In the two-dimensional work,4; 5 Gn+1 and dn+1 in (1) and (14) are, in general, based upon the contact segment and the isoparametric co-ordinates from the contacting point at tn , Xm n. Thus, forces are always applied to the tn contact segment even if the contact segments, at tn and tn+1 , for a slave node dier. In equations (16) – (18), versions of Gn+1 and dn+1 with iteratively updated nn+1 are used, where nn+1 is the surface normal at time tn+1 of the material point cor3;6 a constant G is employed, but responding to Xm n . (In the three-dimensional implementations, its precise denition is omitted.) In the current work G n+1 and dn+1 are used in lieu of Gn+1 and dn+1 . This means that a constant n n+1 is used in place of the actual nn+1 , and the normal contact force fN parallels n n+1 , not nn+1 or an intermediate normal.11 Also, f is applied on the predicted master surface at time tn but at the material point corresponding to the predicted contact location at tn+1 . Therefore, in general, a net moment is imposed on the system because the slave and master forces are applied at dierent spatial locations at tn . While it is possible to use the two-dimensional constructs in three dimensions, the signicant computational cost associated with iteratively performing the additional closest point projections and normal updates does not currently appear justied given the typically small dierences in congurations between any two successive states. The friction model and implementation also dier from the previous treatments. Unlike in Reference 3, no additional Lagrange multipliers are used to explicitly characterize the frictional forces. Rather an implicit approach, in the spirit of Reference 5, is employed to satisfy equations (8) – (11) approximately. By holding ẽ xed in d, dn+1 ([n ; x n+1 ) is a pseudo-linear function in [. Consequently, (15) retains certain linear-like characteristics even when frictional eects are present, and the inclusion of friction does not signicantly alter the form of equation (15) or the ability to obtain a solution with an iterative solver. While the current frictional treatment is approximate, equations (8) –(11) can be solved rigorously by updating them each iteration. (Aside: more precise treatment of the Coulomb relations can be readily incorporated by iteratively updating (11) and by replacing (10) with an accumulated, mass-scaled, frictional force as used in Reference 5.) Since Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 444 E. ZYWICZ AND M. A. PUSO the local updates introduce more non-linearities into (15) and increase the overall computational cost, the approach was not explored in the present work. 3. GOVERNING EQUATIONS FOR RIGID BODIES 3.1. Equations of motion Rigid-body updates are obtained using an explicit central dierence method for translational motions and an implicit integrator for the rotational variables which algorithmically conserves linear momentum, angular momentum, and energy for incremental force free motions.8 The method is essentially the same as that in Reference 8 except a change in incremental rotational variables from the material body frame to the spatial frame is made, i.e. Xn = n n . The alternative representation reduces the computational eort required to solve and perform the updates. In the following sections, it is shown that the appropriate application of torques and velocity constraints can maintain total momentum and energy for a system of contacting rigid bodies free of external forces. This novel treatment provides unconditional stability for such a system. The discrete momentum conservation laws, when modied with the forward increment Lagrange multiplier assumptions, require that for each rigid body Mẋn+1=2 − Mẋn−1=2 = t˜n fn (19) In+1=2 Ẋn+1=2 − In−1=2 Ẋn−1=2 = t˜n cn (20) where In+1=2 = n In−1=2 Tn n = Exp[X n ] (21) (22) Xn = t˜n (Tn Ẋn+1=2 + Ẋn−1=2 )=2 (23) fn = Axn [fn − Gn+1 [n − dn+1 ([n ; xn+1 )] (24) cn = An [fn (25) − Gn+1 [n − dn+1 ([n ; xn+1 )] I denotes the time-dependent spatial inertia tensor dened relative to the centre of mass, Axn and An transform the global forces into generalized forces and torques, ( ) is the skew symmetric tensor formed from the vector ( ), and f, c, x, and Ẋ represent the generalized force, torque, position, and spatial rotational velocity vectors of the centre of mass, respectively. The tensor operator Exp is the exponential map (e.g. see Reference 9). The velocity update formulae, obtained by solving for the unknown quantities, are given by ẋn+1=2 = ẋn−1=2 + t˜n M−1 fn (26) T ˜ Ẋn+1=2 = n I−1 n−1=2 n [t n cn (27) + In−1=2 Ẋn−1=2 ] (Although equation (27) is not explicit since n implicitly depends upon Ẋn+1=2 , its solution is suciently important to justify the local implicit solution procedure developed in the Appendix.) Points located on or within the body are updated via xn+1 = xn+1 + n+1=2 [xn − xn ] Copyright ? 1999 John Wiley & Sons, Ltd. (28) Int. J. Numer. Meth. Engng. 44, 439–459 (1999) EXPLICIT FINITE-ELEMENT CONTACT 445 where xn+1 = xn + tn ẋn+1=2 n+1=2 = Exp[X n+1=2 ] (29) Xn+1=2 ≈ tn Ẋn+1=2 and (30) Note, to maintain compatible rigid-body and deformable update procedures, i.e. in which positions are updated at the whole step and velocities (and inertias) are calculated at the half step, Xn+1=2 is approximated. 3.2. Momentum-conserving contact algorithm When contact occurs between rigid bodies, angular momentum conservation requires that an algorithmic An be used. Consider a system composed of q unconstrained rigid bodies free of any external forces or torques, except those induced from p frictionless contact interactions. Discrete momentum conservation requires that for each body (19) and (20) be satised and that the changes in system momentums, given by, p = pn+1=2 − pn−1=2 = q P i=1 Mi (ẋin+1=2 − ẋin−1=2 ) (31) and, with respect to the origin, by h = hn+1=2 − hn−1=2 q P i i = [xn+1=2 × Mi ẋin+1=2 − xn−1=2 × Mi ẋin−1=2 + Iin+1=2 Ẋin+1=2 − Iin−1=2 Ẋin−1=2 ] i=1 (32) be identically 0. Since only contact forces are present, fni and cin (i.e. (24) and (25)) are expressible as sums of the ‘s(i)’ slave and ‘m(i)’ master contact forces acting on the ith rigid-body surface as fni = m(i) P j=1 fnmj + s(i) P j=1 fnsj (33) sj (csj n × fn ) (34) and cin = m(i) P j=1 mj (cmj n × fn ) + s(i) P j=1 where c is the nominal moment arm, i.e. the vector from the centre of mass to the postulated contact point. Substitution of (33) into (19) permits (31) to be expressed solely in terms of the p contact forces as p = t˜n p P i=1 (fim + fis ) (35) With fim and fis dened by fim ≡ i ni Copyright ? 1999 John Wiley & Sons, Ltd. and fis ≡ −i ni (36) Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 446 E. ZYWICZ AND M. A. PUSO p = 0, and linear momentum is automatically conserved. From equations (20), (26), (29), and (32) – (34), and by use of xn−1=2 = xn − tn−1 ẋn−1=2 =2 and xn+1=2 = xn + tn ẋn+1=2 =2, it can be shown that q P h = [xin × Mi (ẋin+1=2 − ẋin−1=2 ) + t˜n cin ] i=1 q P [xin × fni + cin ] " # m(i) s(i) q P P mj P mj sj sj rn × f n + rn × fn = t˜n = t˜n i=1 i=1 j=1 j=1 (37) where, corresponding to the ith rigid body, rnj = xni + cnj : (38) Since each contact only generates one master and slave contribution, h = t˜n p P j=1 mj sj sj [rmj n × fn + rn × fn ] (39) Equation (39) represents the net rotational impulse imposed on the system by the contact forces mj with respect to the origin. The terms rsj n and rn dene the slave and master contact locations at tn , respectively. In general, the locations are not coincidental, and thus h 6= 0. In order to conserve angular momentum for rigid-body contact, algorithmic moment arms, ĉs and ĉm , are used in An and cn in place of cs and cm , and are constructed so that h = 0, i.e. rsn = rm n . To minimize errors in contact-induced torques, the algorithmic moment arms are chosen as s m sn+1 ) − (xm m ĉm n =c n+1 + [(xn + c n +c n+1 + n)](1 − ) m ĉsn = c sn+1 − [(xsn + c sn+1 ) − (xm n +c n+1 + n)] (40) (41) with = 1=2. Equations (40) and (41) render cs and cm independent of Ẋ when Ẋ n = Ẋn . Unfortunately, these denitions also generate errors in cs and cm that are linear in t and of order (ĉ− c ) × f. While the torque error in either body can be minimized or eliminated by adjusting , in general, it is not possible to eliminate all torque errors and conserve the system angular momentum T m s m via this approach. However, when ẋsn+1=2 −n(nT ẋsn+1=2 ) = ẋm n+1=2 −n(n ẋn+1=2 ) and n × (xn −xn ) = 0, s m use of (40) and (41) render c = c = 0. Thus, two spinning, homogeneous, rigid balls can linearly collide without incurring any algorithmic torque errors. 3.3. Predictor–Corrector The predictor–corrector algorithm for contact between rigid bodies is similar to the one formalized for deformable bodies. The predicted conguration, constructed at tn+1 by assuming [n = 0, is given by n+1=2 [xn − xn ] x n+1 = xn+1 + (42) xn+1 = xn + tn ẋn+1=2 (43) Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) EXPLICIT FINITE-ELEMENT CONTACT 447 ẋn+1=2 = ẋn−1=2 + t˜n M−1 fn (44) n I−1 T ˜ Ẋ n+1=2 = n−1=2 n [t n cn + In−1=2 Ẋn−1=2 ] (45) n is calculated using Ẋ n+1=2 in (23). Again, to evaluate g , the operators Gn+1 and dn+1 Here n+1 x n+1 ) and dn+1 . are replaced with G( The actual position at tn+1 , (28), is expressible in terms of [n and known quantities as xn+1 = xn+1 − tn t˜n M−1 Axn [G n+1 [n + dn+1 ([n ; x n+1 )]− n+1=2 [xn − xn ] (46) (Note, because is a non-linear matrix function, no advantage is gained by using the compound n+1=2 .) Insertion of (46) into (6) gives the non-linear relationship rotation involving n+1=2 )[xn − xn ] gn+1 = g n+1 − Bn [n − bn ([n ) − G Tn+1 (n+1=2 − (47) As with (15), interpenetration is eliminated by solving (47) for [n subject to (7). 4. GAP RESPONSE 4.1. Kinematic response of contacting points When persistent contact occurs between two points, enforcement of g = 0 generates identical kinematic behaviour in the direction parallel to n after only three time steps.4 By use of (3) – (5), the update formulae for the gap function are expressible as gn+1 = gn + tn ġn+1=2 (48) ġn+1=2 = ġn−1=2 + t˜n g n (49) where ġn−1=2 = nT (Ẋsn−1=2 − Ẋm n−1=2 ) and sn − X m g n = nT (X n) (50) If gn = gn+1 = gn+2 = 0, then, via (48) ġn+1=2 = ġn+3=2 = 0, and, from (49), g n+2 = 0. Furthermore, if s gn+3 = 0, then the sequence extends and ġn+5=2 = g n+3 = 0. Consequently, nT (Ẋn+i−1=2 −Ẋm n+i−1=2 )= 0 sn+i − X m ) = 0 for i¿2, and means that the normal velocity component remains the same and nT (X n+i for both the slave and master contact points until separation (or translational movement) occurs. For unconstrained rigid bodies free of any external forces or torques, except those generated by the contact, enforcement of g = 0 generates inelastic collisions. For example, consider a onedimensional problem of two rigid bodies impacting. Assuming the initial contact occurs between tn and tn+1 , satisfying (48) requires that ġn+1=2 60. Physically, this means that the bodies are not separating, and the bodies will be in contact at tn+2 . Continued enforcement in the next and all subsequent time steps yields g = ġ = 0, and therefore the two bodies remain in intimate contact until perturbed by an external force. Had the bodies been elastic, the initial behaviour would be identical. However, the initial contact force would eventually manifest a rarefaction wave that could separate the two bodies. An undesirable result may arise if interpenetration previously existed and is subsequently resolved. Assume gn ¡0 and gn+1 = 0. This condition implies that ġn+1=2 = −gn =tn ¿0 and means that the gap is opening at a rate dictated by gn . If insucient forces exist in subsequent time steps Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 448 E. ZYWICZ AND M. A. PUSO to overcome the nodal momentum imparted by the contact force at tn , e.g. the attached elements reach their elastic limits, non-physical and unbounded gap growth and nodal motion occurs. (Sometimes the colloquial term ‘shooting nodes’ is used to describe this undesirable contact-induced behaviour.) 4.2. Elastic collisions between rigid bodies As just demonstrated, the enforcement g = 0 between rigid bodies results in inelastic collisions and energy dissipation. An algorithmic criterion is now developed to generate elastic collisions between rigid bodies similar to those developed for deformable bodies, e.g. see Reference 11. Consider two unconstrained rigid bodies free of any external forces or torques, except those induced from contact. Assume that the bodies collide at a single point between time tn + and tn+1 . The dierence in the system kinetic energy between time tn−1=2 and tn+1=2 , T , is dened as T = Tn+1=2 − Tn−1=2 2 P i i i i iT iT i iT i = ẋiT n+1=2 M ẋn+1=2 − ẋn−1=2 M ẋn−1=2 + Ẋn+1=2 In+1=2 Ẋn+1=2 − Ẋn−1=2 In−1=2 Ẋn−1=2 i=1 (51) Via (21), (51) reduces to T = 2 P i=1 i i i i i i iT i i T i (ẋiT n+1=2 − ẋn−1=2 )M (ẋn+1=2 + ẋn−1=2 ) + (Ẋn+1=2 − n Ẋn−1=2 ) In+1=2 (Ẋn+1=2 + n Ẋn−1=2 ) (52) By use of (21) and after some algebraic manipulation, (20) can be rewritten as In+1=2 (Ẋn+1=2 − n Ẋn−1=2 ) = t˜n n cn + (1 − n )In+1=2 Ẋn+1=2 (53) where 1 is the second-order identity tensor. Substituting both (19) and the transpose of equation (53) into (52) yields T = = 2 P i=1 2 P i=1 [t˜n fni T (ẋin+1=2 + ẋin−1=2 ) + (t˜n in cin + (1 − in )Iin+1=2 Ẋin+1=2 )T (Ẋin+1=2 + in Ẋin−1=2 )] i iT i i i [t˜n fni T (ẋin+1=2 + ẋin−1=2 ) + (t˜n cinT + ẊiT n+1=2 In+1=2 (n − 1))(n Ẋn+1=2 + Ẋn−1=2 )] (54) = f (Reference 9) and equation (23), (54) simplies to From the identity Exp[f]f T = t˜n 2 P i=1 iT i i [ fniT (ẋin+1=2 + ẋin−1=2 ) + ciT n (n Ẋn+1=2 + Ẋn−1=2 )] (55) When fn and cn are expressed in terms of n, ĉ, and , the algorithmic T is given by T = t˜n nT (ṽm − ṽs ) (56) ṽ = (ẋn+1=2 + Tn Ẋn+1=2 × ĉ) + (ẋn−1=2 + Ẋn−1=2 × ĉ) (57) where Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) EXPLICIT FINITE-ELEMENT CONTACT 449 ˙ and the dynamic gap, f, by Next, dene the algorithmic gap velocity, f, f˙ = nT (ṽs − ṽm ) (58) f = t˜n f˙ (59) and respectively. When is determined such that f = 0, the system kinetic energy, T , is algorithmically conserved. To represent elastic contact between multiple rigid bodies, equation (59) is evaluated at each contact site. An elastic collision occurs between the contacting pair when is found such that f¿0; 60 and f = 0 (60) Note that in the limit as t → 0, (58)–(60) converge to the continuum denition of elastic contact presented by Kane and Levinson12 for the two-body case. Several important dierences exist between the continuum description and the algorithmic representation given by (58)–(60). The proposed method implicitly assumes that the collision occurs between the pseudo-points xns + ĉs and xnm + ĉm , and the collision happens at time tn and in the tn geometric conguration. Only in the limit as t → 0 is the continuum geometric description fully recovered. Unlike (7), (60) can be enforced only between contact pairs that actually collide during the current time step. Otherwise, non-physical interactions are introduced. In light of this, the following procedure is developed for the solution of elastic contact between rigid bodies. To model simultaneous elastic collisions between an arbitrary number of rigid bodies, a two-step predictor–corrector process is employed. The predicted conguration is constructed and potential contact pairs are identied. In the ‘inelastic’ phase, equation (47) is solved for [ subject to equation (7) and used to predict which potential pairs actually collide. (Rigid-on-rigid contact is inelastic in this phase.) It is assumed that only pairs with ¡0 actually make contact during this time step, and therefore, rigid-body pairs whose = 0 are removed from the list of potential contact pairs. In the second or ‘elastic’ step, the elastic contact criterion, (59), is solved subject to the constraint (60) for [n using the reduced pair list. This two-step scheme is consistent with explicit, predictor– corrector, nite element methods, and converges, in the limit as t → 0, with more traditional rigid-body approaches (e.g. References 13 and 14). 5. MATRIX-FREE ITERATIVE SOLVER After identication of all interpenetrating and=or possibly interpenetrating contact pairs, a preconditioned, Fletcher–Reeves CG method is used to nd [n in a newly developed two-phase process. In the traditional inelastic phase, (15) and (47) are combined symbolically into a single vector, given by gn+1 = V([), which is solved subject to (7). Iterations are performed until the convergence criteria are satised. If no rigid-on-rigid pairs exist or if all of their corresponding = 0, then [n = [, and the contact forces are scattered to the appropriate DOF. However, if the inelastic phase terminates with rigid-on-rigid pairs having 6= 0, then rigid-on-rigid pairs whose = 0 are eliminated from the system of equations and an elastic phase is performed. Equations (15) and (59) are combined symbolically into a new vector, given by ĝn+1 = V̂([), and iteratively solved subject to (7), (60), and the convergence criteria. The resultant [ is then scattered to Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 450 E. ZYWICZ AND M. A. PUSO the appropriate DOF. For rigid-on-rigid contact, the inelastic phase determines which pairs are actually in contact, and only these pairs are considered in the elastic phase where the elastic gap idealization is enforced. CG convergence is achieved when the maximum normalized gap, max , is smaller than the tolerance parameter Tol, and k[k, the norm of the change in force, is less than k[k Tol. The normalized gap is dened as = |g | (61) where the characteristic length, 1= , is based upon the element thickness or the segment diagonal. Current tests indicate that Tol = 0·001 yields rapid convergence and acceptable results. Since the resultant gap errors do not appear to inuence the solution or propagate in an unbounded manner, a tighter tolerance appears unnecessary for deformable bodies. However, for rigid-on-rigid contact, energy dissipation is directly related to the accuracy of the iterative solution, and therefore, for rigid-on-rigid pairs, 1= is chosen so that machine precision accuracy is achieved when 6Tol. The following sections describe the CG solver and matrix-free implementation. The same CG algorithm is used in both the inelastic and elastic solution phases. The only dierences are that during the elastic contact phase, V̂ and ĝn+1 are used in place of V and gn+1 , and [ at the start of the elastic phase is initialized with nal value of [ from the inelastic phase. 5.1. Conjugate gradient solver The CG based projection method3 is augmented with an untransformed, linear, Jacobi acceleration, pre-conditioner,15; 16 i.e. diagonal scaling. The diagonal pre-conditioning matrix K is assembled before the iteration process begins and is derived from the denition of Bn . In general, the components are given by Kii = 1 masss + P4 1 j=1 1 h2j mass tt˜ m (62) j except for rigid-on-rigid contact, where they are given by Kii = masss massm tt˜ (masss + massm ) (63) Here h is the local nodal weight and ‘mass’ pertains to the nodal mass, except for rigid-on-rigid contact, where it is the total rigid-body mass. The CG solver is now schematically presented. The iteration process starts with i = 1 and an initial guess for [ given by [(1) = K g , with the restriction that j(1) 60 for all j. Step 1: Update gap, gradient, search direction, and error parameter (i) = V([(i) ) (skip for linear systems when i¿1) gn+1 ( (i) 0; j ∈ (gn+1j ¿0; j(i) ¿0) (i) (i) z → zj = (i) ; otherwise gn+1j q(i) = Kz(i) (i) max = max(j(i) ) Copyright ? 1999 John Wiley & Sons, Ltd. (64) (65) (66) (i) j ∈ (gn+1j ¿0; j(i) ¿0) (67) Int. J. Numer. Meth. Engng. 44, 439–459 (1999) EXPLICIT FINITE-ELEMENT CONTACT Step 2: Adjust search direction (i) zT(i) q(i) = zT(i−1) q(i−1) 0 when i¿1 and otherwise (i−1) = (i−1) or (i−1) (i−1) = max and (i−1) (i−1) ¿0; sj ¡0)} 6= {0} { j: (j (i) + s(i−1) s(i) = −gn+1 451 (68) (69) Step 3: Calculate current step size (i) (i) = min(n+1j =sj(i) ) j ∈ (sj(i) ¡0; j(i) ¡0) max (i) = −zT(i) s(i) − g n+1 ) (70) (71) sT(i) (V(s(i) ) (i) (i) = min(max ; (i) ) (72) Step 4: Update force, iteration counter, and gap (as necessary) [(i+1) = [(i) − (i) s(i) √ k[k(i+1) = |(i) | s(i)T s(i) (73) (74) i=i + 1 (i+1) gn+1 (i) = gn+1 (75) (i) (i) + (V(s ) − g n+1 ) (linear systems only) (76) Steps 1– 4 are repeated until after step 4, max 6Tol and k[k(i+1) 6Tolk[(i) k, or i exceeds ≈ p, where [ ∈ Rp . When the contact system is linear, i.e. friction and rigid bodies are absent, (76) provides a computationally cheaper update for gn+1 than (64). Under these circumstances, equations (76) and (64) are mathematically equivalent because V([) = g n+1 − B[. However, periodic use of (64) in lieu of (76) helps minimize numerical drift. (Note, such reductions cannot be achieved with the Gauss–Seidel iterative strategy used in Reference 5.) If an elastic phase fails to converge within approximately p iterations, a search for redundant or nearly redundant equations is performed. When excessive equations exist, they are eliminated and the CG solver is re-initialized. (In our DYNA3D implementation, redundant equations arise because the search procedure employs a dual pass algorithm and generates a pair of equations when node-on-node contacts occur. Also, global searches sometimes have overlapping domains.) Since convergence often occurs even with multiple sets of redundant equations and the search for redundant equations is rather expensive, it appears more ecient to search for redundant equations only after the rst set of iterations. Convergence of the present CG solver is inuenced by many factors. If the eective B is positive-denite and constant, exact convergence should theoretically occur in less than p + 1 iterations provided that [i ¡0 for all i.17 For [i 60, Sha et al.3 proved that the projection-based CG procedure is convergent, but did bound the number of iterations required to achieve convergence. Unfortunately, in the general case, the eective B may be extremely ill-conditioned or even singular, and no formal convergence proofs exist. To date, the CG solver has demonstrated excellent convergence behavior even when the initial eective B was poorly conditioned. Because the CG projection method eliminates equations Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 452 E. ZYWICZ AND M. A. PUSO Figure 1. Bar impact geometry that do not satisfy (7), excessive constraints are apparently removed during the iteration process. This not only reduces the active number of equations, but also improves the matrix conditioning. When equations are removed, the convergence rate, as measured by the change in max , often increases dramatically. This is especially true when the set of equations in which n ¡0 has been established. Nonetheless, sometimes as many as 6p iterations are required to reach convergence. 5.2. Matrix-free gap update The gap update operations, V([) and V̂([), can be eectively performed in a matrix-free manner. First, the nodal positions resulting from the application of [, denoted by x̂n+1 , are linearly updated via an alternative form of (14), given symbolically by x̂n+1 = x n+1 + j P i=1 xi (i ) (77) where j denotes the number of contact pairs. Although not directly obvious, xi depends only upon xed global quantities and variables associated with the i-th contact pair. Next, when rigid bodies are present, x̂n+1 − x n+1 is used to calculate contributions to t˜n c and, via mass scaling, also to xn+1 . This eliminates the need to segregate deformable and rigid nodes in (77). For each rigid body, (29) and (27) are solved, subject to any imposed kinematic boundary conditions, and the position of rigid-body nodes actively involved in contact are updated in x̂n+1 via (28). ˆ While it is advantageous to Finally, gin+1 or ĝin+1 is updated using x̂n+1 ([) and, as required, ḟ([). calculate and temporarily store some vectors, none of the contact matrices, e.g. B, G, or A, need to be explicitly formed or stored. 6. EXAMPLES The rst example investigates the normal impact of two identical slender bars.18 Both 50 m long bars have a square cross section 1 m wide and deep, and their Young’s modulus, Poisson’s ratio, and density are 100 Pa, 0·0, and 0·01 kg=m3 , respectively. As shown in Figure 1, both bars are meshed with 50 uniform eight-node continuum brick elements. The initial x-component velocity of the left bar is 1 m=s, and the right bar is at rest. The impacting bar surfaces are initially separated by 0·01 m. Figure 2 shows the x-component of the interface velocities for both bars using the proposed Lagrange-multiplier based method and DYNA3D’s default penalty method. The x-component of displacement at the common interfaces is displayed in Figure 3 for both contact idealizations. The penalty method plots have been oset by 0·2 s for visualization purposes. Clearly, the Lagrangian approach represents contact more accurately, i.e. no interpenetration, but slight overshoots are observed in the velocity and, not shown, the normal traction. Furthermore, dierences in the total displacements are observed between the two idealizations. The second example is similar to the rst example, except rigid material idealizations are used and the bars are fully unconstrained. The previous time-step size is employed, and the contact is modelled as fully elastic. The x-component of velocity is shown in Figure 4 for both contact Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) EXPLICIT FINITE-ELEMENT CONTACT 453 Figure 2. Normal interfacial velocity versus time Figure 3. Normal interfacial displacement versus relative time Figure 4. Rigid-body translational velocity versus time idealizations. With the Lagrangian approach, unlike the penalty method, the velocities jump in a single step because the duration of contact is independent of the initial velocity, contact stiness, and time-step size. Because of the dual pass search algorithm used, eight separate node-on-node contacts are identied. Although the resulting system of equations contains four pairs of linearly dependent equations, both Lagrangian solver phases converge in only a few steps and do not Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 454 E. ZYWICZ AND M. A. PUSO Figure 5. Initial and deformed bumper assembly geometries require the duplicate equations to be removed. Although not shown, results generated with the Lagrangian approach are completely independent of the cross sectional discretization even when the two bars are meshed dierently. The third example examines the impact of a 9·47 m=s bumper assembly with a stationary, elastic pole. The assembly consists of an elastic–plastic box-beam structure covered with a thin, elastic– plastic outer shield. The 1766 node model contains 496 eight-node continuum brick elements in the pole and 841 four-node quadrilateral shell elements in the bumper. Contact is permitted between the exterior surfaces of all brick and shell elements (1493 four-node contact segments) via an automatic contact algorithm.19 The initial geometry and deformed geometry at 0·1 s, calculated with frictionless Lagrange-multiplier based contact, are shown in Figure 5. For the frictionless idealization, the number of unknown Lagrange multipliers increases from zero to a peak value of 149 as the shield impacts the pole and the box beam collapses. The average number of multipliers is 93, and, on average, 0·165 CG iterations are required for each multiplier to achieve convergence, i.e. on average approximately sixteen CG iterations are performed per step. The raw computational times and number of time steps necessary to solve the problem with the Lagrange-multiplier and penalty contact methods are summarized in Table I. The CPU time, for simulations with and without friction, is separated into the core solution time (Solution), time spent in contact routines (Contact), and time consumed in initialization and input=output routines (Other). When scaled by the number of steps used to solve the problem with penalty-based contact, the dierences in contact costs between the Lagrange-multiplier method and penalty method are 237 and 343 s for frictionless and frictional contact, respectively. The ratio of core solution cost to Lagrange-multiplier cost is approximately 29 per cent for frictionless contact and 45 per cent for frictional contact. The primary reasons for the higher frictional cost are the gap updates must be calculated explicitly, instead of via (76), and more CG iterations are required to reach convergence. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 455 EXPLICIT FINITE-ELEMENT CONTACT Table I. Computational expense and step count for Lagrange-multiplier and penalty contact methods Frictionless CPU Solution Contact Other Steps Friction ( = 0:1) Lagrangian Penalty Lagrangian Penalty 793 2768 7 23 455 821 2605 7 24 079 832 2942 7 23 938 791 2512 7 23 233 Overall, the additional cost of Lagrange-multiplier contact is 3 and 13 per cent for the frictionless and frictional cases, respectively. Two, traditionally dicult, contact problems concurrently arise in this example. Toward the end of the impact event, the outer shield wraps around the now attened box-beam structure as well as the curved pole. Here four shell layers and one brick layer, each with essentially a different discretization, are simultaneously in intimate contact with each other. The diculties are that ‘shell-to-shell-to-shell’ contact conditions arise and, due to the local curvature and discretization, a set of over constrained and redundant contact equations is generated. Without the proper treatment and enforcement of the constraint inequalities, i.e. (7), the contact equations are unsolvable. Such problems are common in complex applications such as crashworthiness analyses, and, as demonstrated by this example, can be easily and cost eectively solved with the proposed method. In the fourth and nal example, a moving 1 m long rigid rod is conned within a hollow rigid cube. The unconstrained cube, initially stationary, measures 2 m by 2 m by 2 m, each side is 0·1 m thick, and is discretized with one shell element per side. The rod, initially centered within the cube and coincidental with the x-axis, is made from ve collinear beam elements, and its cross section measures 0·1 m by 0·1 m. It has an initial velocity of ẋT = (50; 2; 0) m=s and ẊT = (0; 40; 0) rad=s. Both materials have a density of 7·0 × 10−4 kg=m3 , and the contact is modelled as fully elastic. The analysis is run for 10 s with a xed time-step size of 1 × 10−6 s. The mid-surfaces of the shell elements dene the contact segments, and the six beam nodes are used as the slave nodes. The thickness of both the beams and shells is ignored by the contact algorithm. During the 107 steps, the rod impacts the cube 447 separate times. The initial and nal values of the system linear momentum and angular momentum, with respect to the origin, are listed in Table II. The normalized system kinetic energy, dened as (T=T0 ) − 1, reaches a peak value of 9·46 × 10−11 . Neglecting numerical noise, clearly both the momentums and kinetic energy of the system are conserved. 7. SUMMARY AND CONCLUDING REMARKS A generalized forward increment Lagrange-multiplier based formulations has been constructed along with a matrix-free iterative solution strategy to determine the contact forces between multiple bodies in three dimensions. The predictor–corrector approach, developed here for both deformable and rigid bodies, imposes no Courant stability limitations, demonstrates no processing dependent behaviour (like some kinematic algorithms), and is more accurate than penalty methods. Although Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 456 E. ZYWICZ AND M. A. PUSO Table II. Initial and nal system momentums x-component y-component z-component Linear momentum Initial Final 3·50000000000000 × 10−4 3·50000000000002 × 10−4 1·40000000000000 × 10−5 1·40000000000003 × 10−5 −5·60000002031396 × 10−13 −5·60000015583924 × 10−13 −2·52000000560000 × 10−5 −2·52000000487943 × 10−5 −4·23516473627150 × 10−22 −4·38342938335889 × 10−16 Angular momentum Initial Final 0 −5·03119952383752 × 10−15 problem dependent, the additional computational cost of Lagrange-multiplier based contact over alternative methods does not appear substantial, i.e. less than 15 per cent, and is justied by the improved force resolution, elimination of interpenetration, reduced contact coecient dependency, and its ability to generate elastic contact between rigid components. Several new derivations regarding the algorithmic conservation of momentum and energy for contacting rigid bodies were presented. These include algorithmic moment arms to preserve angular momentum and a velocity criterion to generate elastic collisions. To correctly impose the latter criterion, a two-step solution procedure, compatible with the current iterative strategy, was constructed and demonstrated. While the present paper has addressed many aspects, several issues regarding the current Lagrange-multiplier strategy remain unresolved. To minimize CG computational costs, the proposed CG tolerance for deformable bodies is set rather large, and friction, when present, is implicitly and continuously incorporated in the CG loop. Does the large CG tolerance pose a stability issue, and do alternative treatments of friction oer computational savings or other advantages? Furthermore, how implicit must Lagrange-multiplier methods be to generate suciently accurate and stable results within an explicit nite element framework? For example, Lagrange multipliers may generate deformation that conicts with global constraints and boundary conditions. Can the resultant forces be adjusted in the usual explicit manner or must the constraints and boundary conditions be implicitly incorporated into the contact solution algorithm? Clearly, these issues remain topics for future research. APPENDIX: PRACTICAL COMPUTATION OF THE ANGULAR VELOCITY The solution of equation (27) for Ẋn+1=2 is easily and accurately achieved using successive substitution and, as necessary, a Newton–Raphson iterative scheme. In many practical problems only two to four successive substitution iterations are necessary to converge to machine precision. Because of the structure of the equations, it is computationally faster to iterate for the variable vn+1=2 , dened as vn+1=2 ≡ Tn Ẋn+1=2 than for Ẋn+1=2 directly. Step 0: Initialize variables and estimate new velocities: h = t˜n cn + In−1=2 Ẋn−1=2 Ẋ(0) n+1=2 = 2Ẋn−1=2 − Ẋn−3=2 Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) EXPLICIT FINITE-ELEMENT CONTACT 457 (0) (0) n = Exp[t˜n Ẋn+1=2 ] (0) vn+1=2 = Tn (0) Ẋ(0) n+1=2 Step 1: Successive substitution iterations: i=i+1 Xn(i) = t˜n (i−1) (vn+1=2 + Ẋn−1=2 ) 2 (i−1) n(i) = Exp[X n ] (i) T (i) h vn+1=2 = I−1 n−1=2 n (i) (i−1) (i) E(i) = kvn+1=2 − vn+1=2 k − tolkvn+1=2 k Step 2: Convergence check: If E(i) ¿ 0 and i68; goto Step 1 If E(i) ¡ 0; goto Step 5 Step 3: Newton–Raphson iteration: i=i+1 Xn(i) = t˜n (i−1) (vn+1=2 + Ẋn−1=2 ) 2 (i−1) n(i) = Exp[X n ] (i−1) u(i) = In−1=2 vn+1=2 a(i) = n(i) u(i) − h F(i) = (i) @n(i) u(i) (i) @vn+1=2 + n(i) In−1=2 (i−1) vn+1=2 = vn+1=2 − F−1 (i) a(i) (i) (i−1) (i) E(i) = kvn+1=2 − vn+1=2 k − tolkvn+1=2 k Step 4: Convergence check: If E(i) ¿0; goto Step 3 Step 5: Recover n + 1=2 velocity: (i) Ẋn+1=2 = n(i) vn+1=2 Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) 458 E. ZYWICZ AND M. A. PUSO Remark: The derivative @u=@v can be conveniently calculated by expressing the operator Exp in terms of quaterions and then dierentiating. Skipping several intermediate steps, the nal result is given by ( T ) X u q @q @q @u = t˜ − sin(kXk=2) q0 u − ⊗ u − (q0 u − (qT u)1) +q⊗ @v 2 kXk @X @X where q0 sin(kXk=2) X ⊗ X sin(kXk=2) 1 − + 2 kXk kXk2 kXk 1 sin(kXk=2) X q= 2 kXk @q = @X and q0 = cos(kXk=2) ACKNOWLEDGEMENTS This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract W-7405-Eng-48. The authors beneted greatly from discussions with Prof. R. L. Taylor of the University of California, Berkeley. REFERENCES 1. J. O. Hallquist, G. L. Goudreau and D. J. Benson, ‘Sliding interfaces with contact-impact in large-scale Lagrangian computations’, Comput. Meth. Appl. Mech. Engng., 51, 107–137 (1985). 2. T. Belytschko and M. O. Neal, ‘Contact-impact by the pinball algorithm with penalty and Lagrangian methods’, Int. J. Numer. Meth. Engng., 31, 547– 572 (1991). 3. D. Sha, K. K. Tamma and M. Li, ‘Robust explicit computational developments and solution strategies for impact problems involving friction’, Int. J. Numer. Meth. Engng., 39, 721–739 (1996). 4. M. W. Salveson and R. L. Taylor, ‘Solution of dynamic contact problems by implicit=explicit methods’, Rept. UCRLCR-125780, Lawrence Livermore National Laboratory, Livermore, CA, 1996. 5. N. J. Carpenter, R. L. Taylor and M. G. Katona, ‘Lagrange constraints for transient nite element surface contact’, Int. J. Numer. Meth. Engng., 32, 103 –128 (1991). 6. K. K. Tamma, M. Li and D. Sha, ‘Conjugate gradient based projection: a new explicit computational methodology for frictional contact problems’, Comm. Num. Meth. Engng., 10, 663 – 648 (1994). 7. R. G. Whirley and B. E. Engelmann, ‘DYNA3D: user manual’, Rept. UCRL-MA-107254, Lawrence Livermore National Laboratory, Livermore, CA, 1993. 8. J. C. Simo and K. K. Wong, ‘Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentum’, Int. J. Numer. Meth. Engng., 31, 19 –52 (1991). 9. J. C. Simo and L. Vu-Quoc, ‘On the dynamics in space of rods undergoing large motions—a geometrically exact approach’, Comput. Meth. Appl. Mech. Engng., 66, 125 –161 (1988). 10. J. O. Hallquist, ‘DYNA3D course notes’, Rept. UCID-19899; Rev. 2, Lawrence Livermore National Laboratory, Livermore, CA, 1987, p. 226. 11. T. A. Laursen and V. Chawla, ‘Design of energy conserving algorithms for frictionless dynamic contact problems’, Int. J. Numer. Meth. Engng., 40, 863 – 886 (1997). 12. T. R. Kane and D. A. Levinson, Dynamics: Theory and Applications, McGraw-Hill, San Francisco, 1985, pp. 231–235. 13. J. E. Huag, S. C. Wu and S. M. Yang, ‘Dynamics of mechanical systems with Coulomb friction, stiction, impact and constraint addition-deletion—I’, Mech. Machine Theory, 21, 401– 406 (1986). 14. S. C. Wu, S. M. Yang and J. E. Huag, ‘Dynamics of mechanical systems with Coulomb friction, stiction, impact and constraint addition-deletion—II’, Mech. Machine Theory, 21, 407– 416 (1986). Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999) EXPLICIT FINITE-ELEMENT CONTACT 459 15. J. R. Shewchuk, ‘An introduction to the conjugate gradient method without the agonizing pain’, School of Computer Science Rept., Carnegie Mellon University, Pittsburgh, PA, 1994. 16. T. J. R. Hughes and R. M. Ferencz, ‘Implicit solution of large-scale contact and impact problems employing an EBE preconditioned iterative solver’, Proc. IMPACT 87: Int. Conf. on Eects of Fast Transient Loading in the Context of Structural Mechanics, Lausanne, Switzerland, 1987, pp. 1–17. 17. G. H. Golub and C. F. Van Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, 1989, pp. 520 – 527. 18. S. C. Lovejoy and R. G. Whirley, ‘DYNA3D example problem manual’, Rept. UCRL-MA-105259, Lawrence Livermore National Laboratory, Livermore, CA, 1990, pp. 32–39. 19. R. G. Whirley and B. E. Engelmann, ‘An automatic contact algorithm in DYNA3D for impact problems’, Rept. UCRL-JC-114820, Lawrence Livermore National Laboratory, Livermore, CA, 1993. Copyright ? 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 439–459 (1999)

1/--страниц