Design of energy conserving algorithms for frictionless dynamic contact problems

This paper proposes a formulation of dynamic contact problems which enables exact algorithmic conservation of linear momentum, angular momentum, and energy in (cid:2)nite element simulations. It is seen that a Lagrange multiplier enforcement of an appropriate contact rate constraint produces these conservation properties. A related method is presented in which a penalty regularization of the aforementioned rate constraint is utilized. This penalty method sacri(cid:2)ces the energy conservation property, but is dissipative under all conditions of changing contact so that the global algorithm remains stable. Notably, it is also shown that augmented Lagrangian iteration utilizing this penalty kernel reproduces the energy conserving (i.e. Lagrange multiplier) solution to any desired degree of accuracy. The result is a robust, stable method even in the context of large deformations, as is shown by some representative numerical examples. In particular, the ability of the formulation to produce accurate results where more traditional integration schemes fail is emphasized by the numerical simulations.


INTRODUCTION
The contact problem attracts considerable attention from the computational mechanics community, due in large part to its highly non-linear and discontinuous nature. Indeed, engineering analysts charged with solving such problems will attest that merely achieving convergence of non-linear solution schemes can be di cult under many circumstances. These di culties stem primarily from the fact that contact surfaces are unknown a priori, and must be evolved in a manner consistent with the equilibrium conditions of the two bodies and the Kuhn-Tucker conditions governing contact evolution. These physical requirements imply a variational inequality for the overall system, with admissible variations being constrained by the conditions governing the solution. 1; 2 Many ÿnite element modelling e orts have been addressed most directly to the quasi-static contact problem, where inertial e ects are essentially negligible and the contact conditions can be e ectively devised and implemented on individual conÿgurations, without much regard to the temporal variation of contact kinematic measures. This work, by contrast, concerns itself with treatment of the dynamic contact problem. In general the most prevalent technique for treating dynamics in non-linear solid mechanics has been the semidiscrete ÿnite element method. Formulations of this type can be roughly grouped into two classes: explicit schemes, which are favoured for highly transient problems; and implicit schemes, which are more appropriate for systems dominated by low-frequency behaviour. Implementations of the dynamic contact problem that have been proposed might also be grouped similarly, with explicit treatments appearing for example in References 3-6, and implicit schemes being described, e.g. in References 7 and 8. Roughly speaking, the conceptual approach to contact constraints in many of these works can be described as follows. A semidiscrete ÿnite element system is developed exactly as would be done in a problem with no contact, except that a contact force vector is also included. The contact tractions deÿning this force vector must be subject to some type of contact conditions, and the description of these conditions is critical in determining the properties of the resulting method. A common choice, for example, is to take the contact conditions to be the same as those governing quasi-static response, and to apply any required numerical integration techniques (e.g. for frictional response) exactly as would be done in the quasi-static case. An analogue to this treatment is to be found in many implementations of elastoplasticity, where the integration algorithm for the constitutive equations at each quadrature point does not in general depend on the global time integration scheme, or even on whether the problem is static or dynamic. Stretching this analogy a bit further, we might describe most contact implementations as treating the contact force vector in the global equilibrium equations essentially as an extremely non-linear (and non-smooth) internal force vector, with the 'constitutive law' for the contact being deÿned by the complementarity conditions and the friction law, if present. This general conceptual framework, while widely applied, appears to have several drawbacks. In explicit calculations, when using the penalty method for constraint enforcement, one ÿnds that large penalties cannot be used in a fully explicit contact treatment without a ecting the Courant stability criterion. 5 In softening the penalties, the accuracy of constraint is sacriÿced to some degree, and because central di erence methods possess no numerical damping, the noise generated by the contact treatment will tend to obscure the solution as the calculation proceeds. In fact, one can show 6 that the Lagrange multiplier formulation of the fully explicit contact treatment is singular, calling into question the basic validity of the fully explicit penalty approach. Accordingly, Carpenter et al. 6 have advocated an implicit treatment of the (quasi-static) contact constraints, along with an otherwise explicit time integration of the momentum equations. This approach possesses a well-deÿned Lagrange multiplier formulation which unfortunately couples the equilibrium equations in general. Carpenter and co-workers propose a Gauss-Seidel iteration scheme for constraint enforcement, while Zhong 9 uncouples the constraints by using an alternative discretization he refers to as the 'defense node' approach. Both alternatives appear to be reasonably well behaved in comparison with the fully explicit approach, but the very fact that the contact conditions are treated completely di erently than the rest of the system calls into question our basic understanding of the appropriate dynamic contact constraints.
In the implicit arena the state of the art is somewhat similar. For example, it has long been recognized that use of the second-order accurate trapezoidal rule with a fully implicit treatment of the contact constraints produces signiÿcant oscillations which can become worse as time steps and spatial discretizations are reÿned (see, e.g. Reference 6). Some authors have proposed corrections to the Newmark updates aimed at correcting such di culties; examples are a correction based on wave propagation results for linear elastic materials 7 and a more recent treatment 8 in which new contacts are e ectively forced to be persistent. While in some cases e ective, the need for such 'corrections' suggests some sort of inadequacy in our basic understanding of contact constraints in a dynamic context. This paper makes an attempt at such understanding by formulating an implicit time integration method which is fully conservative, i.e. it conserves all momenta and the total system energy for hyperelastic bodies undergoing perfectly elastic impact events. Since we concentrate here on such conservative systems, we consider only non-dissipative contact and therefore assume frictionless response. Our approach is an extension of the work of Simo and Tarnow, 10 where conservative algorithms were proposed and demonstrated for hyperelastic systems without contact constraints. In the current context, we will see that a Lagrange multiplier formulation of a particular rate constraint on interfaces is completely consistent with global conservation laws. Since the penalty method can be much simpler to apply than Lagrange multipliers in general, we propose a penalty regularization which is unconditionally dissipative, altering the energy conservation property but retaining the stability we desire. We will also demonstrate an augmented Lagrangian update scheme which uses this penalty kernel to reproduce the energy conserving (Lagrange multiplier) solution. In this manner, an algorithm is produced which is stable, and which can conserve energy to any desired degree of accuracy through the augmented Lagrangian iteration procedure. All development will be done without restricting the amount of motion or deformation that can occur, enabling application of the method in a very broad context.
The plan of the paper is as follows. Section 2 outlines the contact problem and discusses the conservation laws and their implications in a continuum context. Section 3 reviews the conservative time integration method proposed by Simo and Tarnow 10 and discusses its extension to the frictionless contact problem. Imposition and regularization of the contact constraints in the context of the conservative algorithm is discussed in Section 4. In particular, Lagrange multiplier, penalty, and augmented Lagrangian algorithms will be proposed in this section. Section 5 brie y discusses the spatial discretization of the contact problem and some associated issues involved in ÿnite element implementation. Finally, Section 6 presents some numerical examples demonstrating the performance of the method and some implications of the conservative scheme.

PROBLEM FORMULATION AND CONSERVATION LAWS
In the following, we brie y discuss the governing equations and contact conditions for the system of interest. The interested reader should consult Laursen and Simo 11 for more details on the continuum formulation of large-deformation contact problems.

Governing equations
We consider the open sets (i) ⊂ R n s d ; i = 1; 2, which represent the reference conÿgurations of two bodies expected to contact during a time interval of interest I = [0; T ]. For each body, we deÿne a portion of the boundary (i) ⊂ @ (i) so that all expected areas of contact are included. Adopting a Lagrangian description of the problem, we designate material points in the contact surfaces as X ∈ (1) and Y ∈ (2) . Writing the unknown conÿguration mappings at any time t ∈ I as D (i) t , i = 1; 2, we can express the spatial positions of the contact surfaces as (i) t = D (i) t ( (i) ). Typical points x ∈ (1) t and y ∈ (2) t are then given by x = D (1) t (X) and y = D (2) t (Y). We assume that the (n sd − 1) dimensional manifolds (i) are parametrized by mappings (i) 0 such that (i) = (i) 0 (A (i) ), i = 1; 2, where A (i) is a parent domain for the surface in question and the mappings (i) 0 are assumed to be su ciently smooth. In particular, considering (2) , we denote points in A (2) as^and write Y = (2) 0 (^) and y = (2) t (^), where (2) t = D (2) t • (2) 0 . Considering any point X ∈ (1) , the normal (impenetrability) contact conditions are written in terms of a gap function g(X; t): ‡ It is deÿned at any time t in terms of a closest point projection ‡ We use the term gap function throughout the manuscript to be consistent with most contact mechanics literature. Actually, because of the sign convention chosen, the term penetration function might be more appropriate since the function g is positive when interpenetration occurs and negative otherwise in the spatial conÿguration: g(X; t) = sign(g(X; t))|g(X; t)| where |g(X; t)| = min and (1) Impenetrability is enforced by the condition g(X; t)60. The point in (2) achieving the minimization in (1) is written as Y, with its counterpart in A 2 denoted as^. It is important to remember that given a point X, the identiÿcation of Y and^will both in general vary with time, so that we will often write Y(X; t) and^(X; t). Given these deÿnitions, a basis can be constructed at each contact point by deÿning (2) 0; (^); = 1; : : : ; n s d − 1 It will be convenient in the following to augment this basis with a surface normal ], which points out of body 2. In three dimensions, ] would be deÿned via where it is assumed that the parameterizations are deÿned so that ] has the proper sense. Although both B and ] are to be associated with X and vary with time, we suppress these arguments in the following to reduce notation. The (Piola) contact traction t (1) (X; t) is resolved into normal and tangential parts via where ] is the outward normal to (2) t at y (thus the inward normal to (1) t ), and P t (1) is the projection of t (1) onto the associated tangent plane. The contact pressure t N (X; t) should be positive for compressive contact. We assume no frictional tractions in the present discussion; therefore P t (1) = 0. The conditions for normal contact can now be written as g(X; t)60 t N (X; t)¿0 t N (X; t)g(X; t) = 0 t N (X; t)ġ(X; t) = 0 (5) which must hold for all X ∈ (1) and for all t ∈ I. Equations (5) 1-3 represent the classical Kuhn-Tucker complementarity conditions between gap and pressure. Equation (5) 4 is a constraint called the persistency condition, and requires that non-zero traction may only be generated during persistent contact. It will be of particular importance in subsequent developments.
With the contact conditions written, we specify the problem to be solved as follows.
Given the following boundary conditions on body force, traction, and boundary displacement: and contact conditions (5) are satisÿed on (1) .
is the prescribed boundary traction, and V (i) 0 is the initial material velocity. The constitutive relations governing P t are at this point left unprescribed. The subregions (i) , (i) ' , and (i) are assumed to be non-intersecting and invariant with time, while satisfying

Variational principle
At any time t ∈ I, one can introduce admissible variations D * (i) on each body and construct a variational principle. Following the development in Laursen and Simo, 11 which should be consulted for more details, this weak form of the equations can be stated as such that for all D * (i) ∈ V (i) ; i = 1; 2: where A (i) t is the material acceleration of the body, and arguments X and t are dropped from t (i) (X; t) for convenience of notation.
The time-dependent solution spaces C (i) t and time-independent variational space V (i) are deÿned such that and One may add the virtual work expressions implied by (9) to deÿne a global variational principle where D t is understood to be the collection of mappings D (i) t ; i = 1; 2 (similarly for D * ). In (12) we also use the notations G (i) (· ; ·) to indicate the sum of the internal virtual work and those of the applied forces. The right-hand side of (12), representing the contact virtual work, can be expressed as a single integral over (1) by requiring the contact forces on either side of the interface to be equal and opposite, which leads to Equation (14) can be simpliÿed further by considering linearized variations of the kinematic quantities, denoted here by the symbol (·): Direct calculation will verify that Using the fact that the frictional traction is zero, (16) can be substituted into (14) to produce

Conservation laws
Before we begin development of the algorithm, it is instructive to consider the global conservation laws in the context of the problem at hand. Speciÿcally, we wish to verify that the linear momentum L t , the angular momentum J t , and the total energy E tot t are globally conserved by the formulation we propose. The paper of Simo and Tarnow 10 should be consulted for more detail on the general approach. Using the current notation, we deÿne the total system linear momentum L t and total system angular momentum J t for any time t as L t := (1) (1) The theorem of power expended states the global energy balance in rate form where the total kinetic energy is deÿned as the stress power is written as the expended power of the external loading is given by and the total power input of the contact stresses is given by is the second Piola-Kirchho stress tensor and is related to Introducing the notations D int t for the internal dissipation function and E int t for the total stored internal energy, we can write the following form of the second law of thermodynamics (i.e. a reduced dissipation inequality): Combining (20) and (25) yields where E tot t is the sum of the kinetic and internal energies of the two bodies in question. In the study considered in this paper, we will make the following assumptions about the system at hand: 1. The bodies are subject to no body forces, so that f (i) t = 0 on (i) × I for i = 1; 2. 2. There are no Dirichlet (displacement) boundary conditions, and the tractions are zero on the Neumann boundaries, so that (i) There is no internal dissipation in the bodies under consideration, so that D int t = 0 for all time t.
Under this set of assumptions, we can examine the conservation properties of our system by making appropriate substitutions into the virtual work expression (13). For example, substituting (13), where W is a constant arbitrary vector on (1) ∪ (2) , gives which in turn implies that dL t =dt = 0 (W is arbitrary). In a similar manner, considering the variation whereŴ is the skew-symmetric tensor whose axial vector is W. In the last step of (28), the ÿrst term disappears because the contraction of a skew symmetric tensor (Ŵ) with a symmetric tensor is zero, while the last term, the contact contribution, disappears because in the frictionless case the traction t (1) (X; t) is collinear with (D (1) t (X) − D (2) t (Y(X; t))). Again, since W is arbitrary, we can conclude that dJ t =dt = 0, so that angular momentum is conserved.
To examine the energy, we consider use of the material velocity ÿeld V t as the variation. Accordingly, we compute In view of our assumptions above, P ext = 0. Using this fact and comparing (29) with (20) leads to the identiÿcation P con (1) . By similar reasoning used to calculate g in (16), the material time derivativeġ(X; t) is given bẏ so that the expression for P con t becomes P con If we now examine equation (5) 4 , the persistency condition for frictionless contact, we are led to conclude which, in view of (26) and the fact that P ext In other words, total energy is conserved for the system at hand as a direct result of persistency condition (5) 4 . Another way to say this is that if we wish to ensure that all contacts are perfectly elastic, so that the net contact power input to the system is zero, the persistency condition should be satisÿed. This observation will be key to the ensuing algorithmic development, where a counterpart of this condition in the temporally discrete framework must be found.

CONSERVATIVE DISCRETIZATION SCHEMES
In this section we extend the energy-momentum conserving scheme of Simo and Tarnow 10 to the contact problem posed in the last section. The reader is referred to that reference for extensive details on the approach, which will be only brie y outlined here. In the following, we pay particular attention to the temporal discretization of the contact conditions, with the aim of algorithmically reproducing the conservation properties outlined for the continuous case in the last section. We subdivide the time interval of interest I into intervals [t n ; t n+1 ], where n is an index on time steps. Given a time step n, we sometimes write t = t n+1 − t n , and note that in general t need not be uniform throughout a problem. Focusing our attention on a typical time interval [t n ; t n+1 ] we will use the notation (·) n to mean the algorithmic (i.e. time discrete) approximation to the continuum variable (·)(t n ). The conservation properties we wish to maintain are where all quantities are as deÿned in the last section.
Simo and Tarnow described their algorithm by making the following deÿnitions based on convex combinations of variables at n + 1 and variables at n: where ∈ [0; 1] is an algorithmic parameter. Considering ÿrst a problem with only one body, and hence no contact constraints, one could summarize their algorithm in the unforced case as Given all data at n, ÿnd D n+1 ∈ C n+1 , such that for all D This algorithm has the following properties: (i) Algorithmic conservation of linear momentum (i.e. satisfaction of (34) 1 ) for any ∈ [0; 1] and for an arbitrary constitutive relation describing the (symmetric) second Piola-Kirchho tensor S; (ii) Algorithmic conservation of angular momentum (i.e. satisfaction of (34) 2 ) for the case where = 1 2 and arbitrary prescription of S; (iii) Algorithmic energy conservation (i.e. satisfaction of (34) 3 ) for the case where = 1 2 , and where S is deÿned according to a gradient of a generic stored energy functionê(C): C n+ÿ is an algorithmic right Cauchy-Green tensor, deÿned via and ÿ ∈ [0; 1] is selected to satisfŷ Such a ÿ must always exist as a direct consequence of the mean value theorem. (iv) Second-order accuracy under the conditions stated in item (iii).
For more details on the proof of these properties the paper by Simo and Tarnow is recommended. In particular, properties (i)-(iii) can be veriÿed by substitution of W, W × D n+ and V n+1=2 , respectively, into (36) 1 , where W is again an arbitrary constant vector. These arguments are directly analogous to those given in the last section for the time continuous case.
Returning once more to the case at hand, where contact constraints are active between two bodies, we propose the following algorithm, which extends that in (36) by including G c (D n+ ; D * ): Given all data at n, ÿnd D n+1 ∈ C n+1 , such that for all D * ∈ V: where and Y n+ (X) is that point of (2) which minimizes D (1) n+ (X)−D (2) n+ (Y) ; i.e. the closest point projection used to deÿne the contact basis is done in the n + conÿguration. One can also infer from this deÿnition that force equilibrium on the interface is enforced in the n + conÿguration.
It is readily veriÿed that properties (i) and (ii) of the non-contact algorithm are retained by this algorithm. These veriÿcations are directly analogous to those for the continuum case and will not be repeated here. In examination of the energy conservation property, we set = 1 2 and examine the equation If both bodies are described by the constitutive law outlined in equations (37)-(39), one can show by the same arguments used without contact that In view of (43), it is clear that the following must be satisÿed for energy to be conserved: Since we consider frictionless response, and have already made the condition that the basis for t (1) will be deÿned in the n + 1 2 conÿguration, we can simplify G c (D n+1=2 ; V n+1=2 ) as whereg n+1=2 ( X), an algorithmic gap rate, is deÿned as Finally, we can conclude from (46) that if the following algorithmic persistency condition is satisÿed pointwise on (1) , then global energy will be conserved for the two-body contact system: We now have the algorithmic counterpart of (5) 4 appropriate for energy conservation. Note that the remainder of the contact conditions, i.e. equations (5) 1-3 , are to this point unenforced by our algorithmic formulation. We discuss some alternatives for carrying out this step in the next section.
Remark. Examination of equations (47) and (48) shows that the pointwise contact condition necessary to obtain conservation properties only involves quantities associated with X and Y n+1=2 (X). In particular, no history terms associated with the point X need be stored and utilized in calculation ofg n+1=2 , with only algorithmic material velocities and the surface normal at the n + 1 2 conÿguration being involved. This is particularly advantageous in large deformation, large slip problems, which are characterized by frequent changes in the elements contacted by individual points X. The fact that the conservation properties are una ected by such events is a crucial feature of the algorithm we propose.

DISCRETE CONTACT CONSTRAINTS
In the following we discuss three alternative formulations of the contact constraints appropriate for the above conservative framework. In view of the role of the normal pressure-gap rate complementarity condition in energy conservation, we propose the following expression of the contact conditions in the continuum setting. In contrast to (5), they emphasize complementarity conditions between the gap rate and the pressure: We note that g is precluded from becoming positive by equation (50) 2 , which simply states that if g = 0 the only change that can occur is a negative one. In the continuum case, one can verify that conditions (49) and (50) place the same physical restrictions as equations (5). However, these new conditions appear to be more readily accommodated into a conservation context, as we demonstrate below.

Lagrange multiplier formulation
In the time discrete setting we now consider, we let t N = N , where N is the Lagrange multiplier ÿeld on (1) satisfying: whereg n+1=2 is as deÿned in (47) and g n is deÿned via We note that the impenetrability condition g60 is only enforced in the limit as t → 0; in particular, examination of (51) and (52) will show that computation of a non-zero contact traction t N will lag one time step behind the detection of an interpenetration. Conditions (52) will then act to preclude further violation of the constraint by directly enforcing the gap rate to be nonpositive. Note that this prescription, when used in conjunction with the global algorithm of the last section, preserves all the conservation properties discussed. In particular, energy is conserved because Ng n+1=2 = 0 at all contact points and for all conditions of contact.

A penalty regularization
Equations (51) and (52) have the advantage that they are readily regularized, enabling a penalty solution which in some cases may be more convenient. Consider the following prescription for t N : where N ¿ 0 is a penalty parameter (physically speaking, a viscosity) and · denotes the MacAuley bracket, or positive part of the operand. One notes that in the limit as N → ∞ conditions (51) and (52) are reproduced. However, it is no longer true that t Ng n+1=2 = 0, so we cannot expect exact energy conservation in this case. This makes sense; in fact, penalty methods rely on the ÿnite, non-zero energy associated with the penalization to enforce the constraint as the penalty is increased. We would like to ensure, however, that use of (54) with the otherwise conservative framework does not increase energy. This is readily veriÿed by using equations (43), (44) and (46): which shows that the algorithm is unconditionally dissipative. From an energy standpoint the integrator thus remains stable. The dissipated energy can also be restored via augmented Lagrangian iteration, as discussed below.

Augmented Lagrangian iteration
Consider the following expression for t N , based on an augmented Lagrangian augmentation of the penalization in (54): where (k) N is a ÿxed iterate for the Lagrange multiplier satisfying (49) and (50). One may therefore consider solving the following iterative sequence of problems within a time step, in which we begin with some initial estimate (0) N of the multipliers on (1) , and proceed with iterations (k) until some tolerance is satisÿed: where G(· ; ·) is as deÿned in (41), D (k) n+1=2 is as deÿned in (35), and G c is deÿned via n+1=2 (X))) d (1) 2. Update the multipliers on (1) : 3. IF convergence achieved EXIT, ELSE return to 1.
As we will show in the numerical examples, this algorithm can reproduce conditions (51) and (52) arbitrarily closely as iterations (k) proceed. Several convergence criteria might be proposed to deÿne step 3; for example, the error in E tot n+1 − E tot n might be a reasonable choice for deÿning a convergence tolerance for these iterations.

SPATIAL DISCRETIZATION AND IMPLEMENTATION
In this section we brie y present an overview of the spatial discretization process which, when applied to the conserving algorithm of the last section, produces the same conservation properties in the fully discrete setting.

Finite element discretization
In general terms, the ÿnite element discretization is achieved by introducing D h t and D * h , the ÿnite dimensional approximations of D t and D * . These lie in the discrete spaces C h t and V h such that and where d A t is the vector-valued nodal value of the conÿguration mapping, N A is a ÿnite element shape function with domain (1) ∪ (2) , and c A are nodal constants. Substitution of these ÿnite dimensional approximations into the time discrete weak form and enforcing it for arbitrary combinations of c A yields the following fully discrete equations to be solved in each time step [t n ; t n+1 ]: where M is the global mass matrix, d is the global solution vector, F c is the contact force vector, and F int is the internal force vector. We retain the discrete counterparts of D n+1=2 and V n+1=2 used previously, e.g.: and remark that the nodal force F int n+1=2 A associated with node A is given by where S h is as deÿned in (37)-(39). In giving a general expression for F c we summarize the main results, referring the reader to Laursen and Simo 11 for more details on the general procedure. The integral G c can be approximated via where n int is the total number of contact quadrature points on (1) , W l is a weight of integration for quadrature point l, and j l is the jacobian resulting from the local to global transformation used to describe (1) . F c can then be expressed as an assembly of individual quadrature point contributions as follows: where A is the standard ÿnite element assembly operator. In the case of interest, frictionless contact, f l c n+1=2 takes the following form if nodal quadrature is used in (65): where n l = [] l n+1=2 ; −N 1 (^l n+1=2 )] l n+1=2 ; · · · ; −N nel (^l n+1=2 )] l n+1=2 ] T with N a ; a = 1; · · · ; nel representing the shape functions interpolating the element surface containing the projection Y n+1=2 (X l ). The contact pressure t N can be described by any of the representations given in the last section. Calculation of the contact sti ness requires exact linearization of the contact force vector. We omit this calculation here, and refer the interested reader to Laursen and Simo 11 for elaboration on the general procedure.

NUMERICAL EXAMPLES
In presenting some results obtained with the proposed algorithm, we give comparisons with results obtained via two other prevalent strategies: the Newmark method 12 and the Hilber-Hughes-Taylor, or HHT method. 13 As frequently implemented for contact problems, these integrators could be summarized as follows: HHT where , ÿ and represent algorithmic parameters throughout. As discussed in the introduction, the contact constraints are commonly imposed in both strategies without explicit consideration of rate conditions, so that F c is assembled from tractions t N satisfying where g in equations (70) is evaluated from the n + 1 conÿguration in the case of Newmark and from the n + conÿguration in the case of HHT.

Impact between identical bars
In the ÿrst example we consider the axial impact of two identical straight elastic bars, which provides a simple yet illustrative demonstration of the algorithm's performance. In the initial state the two bars are collinear but out of contact as shown in Figure 1. Bar A is initially moving at a uniform velocity of 1 unit while bar B is at rest initially. The properties are: density = 1 unit, area of cross section a c = 1 unit, length l = 10 unit, Young's modulus E = 1 unit and Poisson's ratio = 0. The problem is driven only by the initial conditions with di erent amounts of initial separation considered between the bars.
The solution for displacements and velocities of the impacting ends of the two bars based on physical observations is plotted in Figure 2 for the case of small strains. The numerical solution is obtained using small strain bilinear elastic elements, with 100 elements in each bar and with a time step of t = 0·1 unit. We consider four algorithms for the temporal integration: the proposed conservative scheme, the trapezoidal rule (second-order accurate Newmark with ÿ = 0·25, = 0·5), Newmark with maximum high-frequency dissipation (ÿ = 0·49, = 0·9) and HHT with optimal high-frequency damping (ÿ = 0·3025, = 0·6, = 0·9). Solutions for these cases are shown in Figures 3-6. Unless otherwise stated, the method of augmented Lagrangians is used with iterations until the change in the multiplier value is less than 1 per cent.
Examination of the results shows that dissipative Newmark ( Figure 4) and HHT ( Figure 5) produce oscillations shortly after the initial contact event, characterized by the two bars coming into and out of contact. The high-frequency dissipation of these two integrators eventually damps out this behaviour, so that a reasonable solution is produced. The trapezoidal rule possesses no such dissipation, so that these numerically induced contact oscillations persist throughout the solution ( Figure 3). Notably, the conservative algorithm possesses no dissipation either, and yet does not excite these spurious modes due to the conservative treatment of the contact conditions ( Figure 6). In fact, one can show that up to the treatment of the contact conditions, the trapezoidal rule and the conservative algorithm are identical for small-strain elasticity problems. We therefore conclude that the key issue in this problem is the stable treatment of the contact constraints. The slight interpenetration of the bar ends evident in Figure 6 also tells us that insistence on complete impenetrability is unnecessary for accurate prediction of energy and momentum transfer across the interface. Further insight can be gained by looking at Figure 7, which shows the total system energy for the four cases. Clearly, use of constraints (70) with the trapezoidal rule results in a net energy gain with each new impact event. This energy remains in the system and continues to grow until the impacts stop, after which the system energy is constant. The dissipative integrators damp out the contact oscillations as discussed above, but lose system energy in doing so. The conservative algorithm obtains results no more noisy than those provided by the dissipative integrators, while conserving system energy. Importantly, the conservative algorithm in no way forces the contact to be persistent, but instead allows the conditions to evolve as driven by the momentum equations. It would therefore appear that formulations requiring such persistency of new contacts, as in Taylor and Papadopoulos, 8 may pose needless constraints which could in some problems be non-physical.
Finally, all integrators display oscillations after the separation of bars. This results from the numerical dispersion of the solution near discontinuities causing the higher-frequency oscillations to travel at lower speeds. A common way to deal with such numerical dispersion is to damp out these oscillations with artiÿcial viscosity (see for example References 14 and 15). An interesting issue we intend to pursue in future work is the incorporation of conserving contact treatments into dissipative global schemes, creating algorithms with high-frequency dissipation which integrate the contact conditions in a stable manner.

Impact of a ring against rigid surface
This problem is motivated by a system originally discussed in Wriggers et al. 16 An elastic ring is thrown with an initial velocity of 2 units at a 45 • angle to a at rigid surface, as depicted in Figure 8. The material properties and dimensions of the ring are as follows: Young's modulus E = 10 2 units, density = 0·01 unit, Poisson's ratio = 0·0001, outer radius r o = 10 units and inner radius r i = 9 units. Bilinear large strain elements are used to discretize the ring, with the mesh consisting of 64 elements. Results are obtained for the trapezoidal rule, the dissipative Newmark method and the proposed conservative strategy, using a time step t = 0· 2 in all cases. The ring conÿgurations at various time intervals and the total energy are plotted in Figures 8 and 9, respectively.
The conÿguration plot shows large deformations during and after impact and an unstable blow-up of the simulation for the trapezoidal rule. The dissipative Newmark method provides reasonable results for small times, but damps out the low-mode deformation of the ring for larger times, preserving only the rigid-body motion with any accuracy. The conservative scheme is seen to be clearly superior to either of the above two approaches. Examination of Figure 9, a plot of the system energy for all three simulations, reveals clearly the energy instability in the trapezoidal rule simulation, the loss of energy for the dissipative scheme, and the exact conservation for the current method.

The carrom problem
This problem involves perfectly elastic and frictionless impact between rigid bodies. The system can be visualized as a 2-D version of a pool table, where a plastic disk called the striker is Figure 11. Total energy vs. time for the conservative scheme, the trapezoidal rule and dissipative Newmark; carrom problem manually struck with the aim of hitting other disks and causing them to go into holes in the four corners of the carrom board. The model simulates the motion of the striker in the absence of other disks as it hits the sides of the carrom board assuming no friction (see Figure 10). The sides of the board are simulated with four elements and the striker is placed within the closed space and subjected to an initial velocity. The material properties (contrived so that deformations will be negligible) and dimensions are as follows: striker-Young's modulus E = 10 2 units, Poisson's ratio = 0·0001, density = 0·1 unit and radius r = 1 unit; carrom board-Young's modulus E = 10 4 units, Poisson's ratio = 0·0001, density = 10 2 units, length of inner side l i = 13 units and length of outer side l o = 15 units. A time step of t = 4 units was utilized in each of three simulations, performed using the trapezoidal rule, dissipative Newmark and the conservative scheme.
The initial condition causes the striker to hit the lower side in the middle and at an angle of 45 • , with a velocity of 0·1 units. The results in Figure 10 show that in case of the trapezoidal rule the rebound is too sharp due to gain of energy; with an eventual blow-up of the solution. The dissipative Newmark result has the opposite e ect, with the angle of rebound being too at and the system energy being reduced. Finally, the conservative scheme correctly predicts all rebounds to be at 45 • after the initial impact, giving the expected diamond-shaped trajectory of the striker within the carrom board. The plot of total system energy in Figure 11 is consistent with these observations, and again veriÿes the energy conservation property.

SUMMARY AND CONCLUSIONS
This work extends the idea of algorithmic energy and momentum conservation to encompass systems featuring mechanical contacts. Importantly, these conservation properties hold for all changing conditions of contact: new contact, release and persistent. The result is an algorithm which is stable in the complete absence of any dissipation. Although the framework we have developed is by deÿnition only applicable to conservative systems featuring frictionless contacts, the work has revealed that careful consideration of algorithmic contact work input is crucial to the construction of stable integration procedures.
It is fair to observe that most physical systems of interest are not conservative, and that few impact events are reasonably modelled as perfectly elastic. E ective algorithms are needed for problems featuring physical dissipation arising from inelasticity and /or interface friction. Furthermore, discretization of continuous systems introduces non-physical modes whose energy content should ordinarily be damped out in transient simulations. When viewed in this light, perhaps the most important product of this work is its identiÿcation of algorithmic persistency conditions that enable stable numerical integration through conditions of changing contact. We intend to extend this work so that interface dissipation, consistent with the frictional model used, can be accurately produced by the algorithm in each time step of a transient simulation. Future work will also explore the use of such energetically consistent contact formulations in conjunction with otherwise dissipative global integrators (such as HHT). In this manner both physically and numerically dissipative equations could be integrated without spurious energy input from the contact conditions.