close

Вход

Забыли?

вход по аккаунту

?

290

код для вставкиСкачать
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)
Документ
Категория
Без категории
Просмотров
2
Размер файла
218 Кб
Теги
290
1/--страниц
Пожаловаться на содержимое документа