# Lagrange constraints for transient finite element surface contact

Nicholas J. Carpenter, Robert L. Taylor, Michael G. Katona

► To cite this version:

Nicholas J. Carpenter, Robert L. Taylor, Michael G. Katona. Lagrange constraints for transient finite element surface contact. International Journal for Numerical Methods in Engineering, 1991, 32 (1), pp.103-128. 10.1002/nme.1620320107. hal-01389918

HAL Id: hal-01389918

https://hal.science/hal-01389918

Submitted on 30 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# LAGRANGE CONSTRAINTS FOR TRANSIENT FINITE ELEMENT SURFACE CONTACT

NICHOLAS J. CARPENTER

TRW Inc., 1 Space Park, Redondo Beach, California 90278, U.S.A.

ROBERT L. TAYLOR

Department of Civil Engineering, University of California at Berkeley, Berkeley, California, U.S.A.

MICHAEL G. KATONA

HQ AFESC/RD, Tyndall AFB, Florida, U.S.A.

## SUMMARY

A new approach to enforce surface contact conditions in transient non-linear finite element problems is developed in this paper. The method is based on the Lagrange multiplier concept and is compatible with explicit time integration operators. Compatibility with explicit operators is established by referencing Lagrange multipliers one time increment ahead of associated surface contact displacement constraints. However, the method is not purely explicit because a coupled system of equations must be solved to obtain the Lagrange multipliers. An important development herein is the formulation of a highly efficient method to solve the Lagrange multiplier equations. The equation solving strategy is a modified Gauss-Seidel method in which non-linear surface contact force conditions are enforced during iteration. The new surface contact method presented has two significant advantages over the widely accepted penalty function method: surface contact conditions are satisfied more precisely, and the method does not adversely affect the numerical stability of explicit integration. Transient finite element analysis results are presented for problems involving impact and sliding with friction. A brief review of the classical Lagrange multiplier method with implicit integration is also included.

## 1. INTRODUCTION

Surface contact kinematic conditions can be enforced by prescribing displacement constraints to prevent structural or continuum domains from overlapping and to control surface contact sliding. Surface contact also involves contact force conditions, typically consisting of a tension limit condition for normal forces and a friction limit condition for tangential forces. Because of the non-linearity associated with surface contact force conditions, an iterative strategy is generally required to obtain a precise solution.

Lagrange multiplier methods and penalty function methods are the two most common used approaches to enforce finite element surface contact displacement constraints. Lagrange multiplier methods are alternatively referred to as mixed or hybrid variational methods by some authors, and penalty methods are commonly referred to as 'contact', 'gap', or 'joint' element methods. For transient analyses by explicit integration, penalty methods have received the most attention in the literature and in commercial finite element programs. Some of the recent work involving the use of Lagrange methods and penalty methods for finite element surface contact is
found in References 1-3, 6, 8, 10, 13 and 14. And an extensive reference to the literature on this subject is given in Reference 14.

The primary focus herein is on Lagrange multiplier methods. The presentation begins in Section 2 with a preliminary discussion of the finite element equation of motion and two dimensional surface-contact. A brief review of the classical Lagrange multiplier method is presented in Section 3, where it is shown that the classical method is not compatible with explicit integration operators. In Section 4 an alternative formulation that is compatible with explicit operators is presented, which is referred to as the 'forward increment Lagrange multiplier' method. Compatibility with explicit operators is established by referencing the Lagrange multipliers one time increment ahead of the associated surface contact displacement constraints.

A one dimensional impact example involving a single contact constraint is presented in Section 5. Finite element analysis results are compared for alternative methods of enforcing the contact constraint. Enforcing the constraint by the forward increment Lagrange multiplier method leads to a highly accurate and well behaved solution. By contrast, the performance of the classical Lagrange method with an implicit integration operator is shown to be poor. Results obtained using a penalty function method to enforce the constraint are also presented.

A two dimensional finite element surface contact formulation based on the forward increment Lagrange multiplier method is developed in Sections 6, 7 and 8. Kinematic conditions and displacement constraints are considered in Section 6. In Section 7 the Gauss-Seidel method is introduced to solve the coupled forward increment Lagrange multiplier equations. The Gauss-Seidel method is then modified in Section 8 to allow for the enforcement of contact force conditions during iteration. In Section 9 the forward increment Lagrange multiplier method and the modified Gauss-Seidel method are exercised to solve a two dimensional surface impact example. A two dimensional finite element sliding problem is then presented in Section 10, followed by closing remarks in Section 11.

The forward increment Lagrange multiplier method is an extension of the ideas presented in Reference 14. The most significant contribution in the present paper is the formulation of an efficient method to solve the coupled forward increment Lagrange multiplier equations that arise in two dimensional surface contact.

## 2. CONSTRAINED EQUATION OF MOTION

The finite element semi-discretized equation of motion is expressed in general form as

$$ M\ddot{u} + F(u, \dot{u}) = R \quad (1) $$

in which M is the mass matrix, u is the vector of displacement degrees of freedom, $\dot{u}$ is velocity, $\ddot{u}$ is acceleration, F is the internal force vector and R is the external force vector. In addition to the usual prescribed boundary conditions, it is assumed that the solution of equation (1) is also subject to surface contact displacement constraints.

An illustration of two dimensional surface contact between bodies that are spatially discretized using low order continuum finite elements is shown in Figure 1. Contactor nodes are denoted by a 'C' and target nodes by a 'T'. Displacement constraints are prescribed to prevent the contactor nodes from penetrating the target domain and to control tangential sliding of contactor nodes along target surfaces. These constraints may be expressed as

$$ G\{u + X\} = 0 \quad (2) $$

where X is the material co-ordinate vector, the sum of u and X is the spatial co-ordinate vector, and G is a surface contact displacement constraint matrix.
<figure><img src="image_2.png" /><figcaption>Figure 1. Two dimensional finite element surface contact illustration</figcaption></figure>

The components of G are typically unknown *a priori* and generally change as displacement and deformation occur. Starting from a configuration in which the surfaces are separated, the motion of contactor and target nodes must be tracked so that displacement constraint components can be introduced in G as contact occurs. During contact the components of G may change with time as required to ensure that the associated contact force reactions satisfy contact force conditions. For example, if a force component normal to a target surface approaches a tension limit force condition, then the associated displacement constraint must be eliminated to allow surface contact separation. Similarly, if a tangential force component approaches a friction limit force condition, then the associated tangential displacement constraint must be relaxed to allow sliding. The components of G change with time when sliding occurs.

It is convenient to first consider the less complicated problem of treating equations (2) as known linear equality constraints that do not change during an integration time increment. This assumption simplifies the initial discussion of Lagrange multiplier methods as presented in the following two sections. A method for relaxing displacement constraints to enforce two dimensional surface contact force conditions is formulated in Section 8.

### 3. LAGRANGE MULTIPLIER METHOD

Lagrange multipliers may be introduced into the equation of motion to give

$$ M\ddot{u} + F(u, \dot{u}) + G^T\lambda = R \quad (3) $$

where the components of the Lagrange multiplier vector $\lambda$ are the surface contact forces. The Lagrange multiplier method proceeds by treating $\lambda$ as unknown and solving equations (2) and (3) simultaneously.

For an elementary small displacement problem in which internal forces are strain-rate independent and proportional to displacement, the constrained equation of motion referenced to time $t_{n+1}$ is

$$ M\ddot{u}_{n+1} + Ku_{n+1} + G_{n+1}^T \lambda_{n+1} = R_{n+1} \quad (4a) $$

$$ G_{n+1}\{u_{n+1} + X\} = 0 \quad (4b) $$

Equations (4) may be solved by direct time integration, see for example Hughes *et al.*⁷ and Bathe and Chaudhary.¹ Herein, the following second order direct time integration operator is considered,

$$ u_{n+1} = q_0 + b_0 \Delta \ddot{u}_{n+1} \quad (5a) $$

$$ \dot{u}_{n+1} = q_1 + b_1 \Delta \ddot{u}_{n+1} \quad (5b) $$

$$ \ddot{u}_{n+1} = q_2 + b_2 \Delta \ddot{u}_{n+1} \quad (5c) $$
where

$$ q_0 = u_n + h\dot{u}_n + \frac{1}{2}h^2\ddot{u}_n \quad (5d) $$

$$ q_1 = \dot{u}_n + h\ddot{u}_n \quad (5e) $$

$$ q_2 = \ddot{u}_n \quad (5f) $$

and

$$ b_0 = \frac{1}{2}h^2\beta_0 \quad (5g) $$

$$ b_1 = h\beta_1 \quad (5h) $$

$$ b_2 = 1 \quad (5i) $$

$$ h = t_{n+1} - t_n \quad (5j) $$

This temporal discretization is equivalent to the well-known Newmark method.¹² However, the operator form presented here is formally known as the Beta-2 method, which is a subset of the generalized Beta-m method developed by Katona and Zienkiewicz.⁹ Two well-known Beta-2 methods are:

(i) the constant-average-acceleration method, also known as the trapezoidal rule; which corresponds with $\beta_0 = \beta_1 = \frac{1}{2}$;

(ii) a single step version of the central difference method, which corresponds to $\beta_0 = 0$ and $\beta_1 = \frac{1}{2}$.

Moreover, a Beta-2 operator is referred to as implicit if $\beta_0 \neq 0$, and explicit if $\beta_0 = 0$.

Substituting equations (5) into equations (4) leads to the following incremental equation of motion:

$$ \begin{bmatrix} [b_2 \mathbf{M} + b_0 \mathbf{K}] & \mathbf{G}_{n+1}^{\mathrm{T}} \\ b_0 \mathbf{G}_{n+1} & \mathbf{0} \end{bmatrix} \begin{Bmatrix} \Delta \ddot{\mathbf{u}}_{n+1} \\ \lambda_{n+1} \end{Bmatrix} = \begin{Bmatrix} \mathbf{R}_{n+1} - \{\mathbf{M}\mathbf{q}_2 + \mathbf{K}\mathbf{q}_0\} \\ - \mathbf{G}_{n+1} \{\mathbf{q}_0 + \mathbf{X}\} \end{Bmatrix} \quad (6) $$

For the surface contact constraints, the rows of $\mathbf{G}_{n+1}$ are linearly independent. Therefore the above system of equations is non-singular if $[b_2\mathbf{M} + b_0\mathbf{K}]$ is non-singular and $b_0 \neq 0$. Conversely, the system is singular for $b_0 = 0$, thereby excluding admissibility of explicit integration operators.

If the non-linearity associated with surface contact force conditions is considered, then an iterative form similar to equations (6) arises. A general two dimensional non-linear surface contact formulation of this type has been developed by Bathe and Chaudhary,¹ wherein a system of equations similar to equations (6) is re-solved for each iteration. The behaviour and accuracy of the method are good for static and slow transient problems. However, the method is ill behaved when inertial forces are relatively large. This is demonstrated by a one dimensional finite element example in Section 5.

## 4. FORWARD INCREMENT LAGRANGE MULTIPLIER METHOD

It is easy to show that the contact forces $\lambda_{n+1}$ directly influence the forward differences $\Delta\ddot{\mathbf{u}}_{n+1}$, but act one time step too late to influence the history terms, $\mathbf{u}_n$, $\dot{\mathbf{u}}_n$ or $\ddot{\mathbf{u}}_n$. Furthermore, the forward difference terms $\Delta\ddot{\mathbf{u}}_{n+1}$ are absent in the second order incremental equation

$$ \mathbf{u}_{n+1} = \mathbf{u}_n + h\dot{\mathbf{u}}_n + \frac{1}{2}h^2\ddot{\mathbf{u}}_n + b_0\Delta\ddot{\mathbf{u}}_{n+1} \quad (7) $$
for explicit integration because $b_0 = 0$. Therefore, the contact forces $\lambda_{n+1}$ have no influence on displacements $\mathbf{u}_{n+1}$, which explains the explicit operator singularity of equations (6).

An alternative formulation that is compatible with explicit integration is developed by relating displacement constraints at time $t_{n+1}$ with Lagrange multipliers at time $t_n$. This approach is referred to here as the forward increment Lagrange multiplier method, for which the incremental equation of motion is

$$M\ddot{\mathbf{u}}_n + F(\mathbf{u}_n, \dot{\mathbf{u}}_n) + \mathbf{G}_{n+1}^T \lambda_n = \mathbf{R}_n \quad (8a)$$

$$\mathbf{G}_{n+1}\{\mathbf{u}_{n+1} + \mathbf{X}\} = \mathbf{0} \quad (8b)$$

For explicit integration of equations (8), consider the well-known multi-step central difference method, which is given by

$$\dot{\mathbf{u}}_n = \frac{1}{2h} \{\mathbf{u}_{n+1} - \mathbf{u}_{n-1}\} \quad (9a)$$

$$\ddot{\mathbf{u}}_n = \frac{1}{h^2} \{\mathbf{u}_{n+1} - 2\mathbf{u}_n + \mathbf{u}_{n-1}\} \quad (9b)$$

Direct substitution from equations (9) into equations (8) leads to

$$\mathbf{u}_{n+1} = \mathbf{u}_{n+1}^* + \mathbf{u}_{n+1}^c \quad (10a)$$

where

$$\dot{\mathbf{u}}_n = \frac{1}{h} \{\mathbf{u}_n - \mathbf{u}_{n-1}\} \quad (10b)$$

$$\mathbf{u}_{n+1}^* = h^2 \mathbf{M}^{-1} \{\mathbf{R}_n - F(\mathbf{u}_n, \dot{\mathbf{u}}_n)\} + 2\mathbf{u}_n - \mathbf{u}_{n-1} \quad (10c)$$

$$\lambda_n = [h^2 \mathbf{G}_{n+1} \mathbf{M}^{-1} \mathbf{G}_{n+1}^T]^{-1} \mathbf{G}_{n+1} \{\mathbf{u}_{n+1}^* + \mathbf{X}\} \quad (10d)$$

$$\mathbf{u}_{n+1}^c = -h^2 \mathbf{M}^{-1} \mathbf{G}_{n+1}^T \lambda_n \quad (10e)$$

The incremental solution of equations (10) proceeds by first calculating the usual explicit displacement update for all degrees of freedom, including surface contact nodes. These displacements $\mathbf{u}_{n+1}^*$ are denoted with a right superscript asterisk to indicate that the effect of surface contact forces $\lambda_n$ is neglected. Following the explicit update, the spatial co-ordinates $(\mathbf{x}_{n+1}^* = \mathbf{u}_{n+1}^* + \mathbf{X})$ are processed by a search algorithm to identify all surface contactor nodes that have penetrated target surfaces. Based on the penetration search information a constraint matrix $\mathbf{G}_{n+1}$ can be assembled and equations (10d) can be solved for the contact forces $\lambda_n$. The incremental displacements $\mathbf{u}_{n+1}^c$ associated with these contact forces are then calculated from equation (10e), and the total displacements $\mathbf{u}_{n+1}$ are calculated from equation (10a).

An efficient Gauss-Seidel iterative approach to combining and solving equations (10d) and (10e) for two dimensional surface contact is developed in Section 7. A modification is then introduced in Section 8 which relaxes the displacement constraints as required to ensure that $\lambda_n$ satisfies surface contact force conditions.

An incremental algorithm based on equations (10) is given in the Appendix as Algorithm 1. The incremental algorithm for Beta-2 explicit methods is similar to the multi-step central difference algorithm and is discussed briefly in Section 10. Two advantages of the Beta-2 explicit approach are that numerical damping can be prescribed by simply selecting an appropriate value of $\beta_1$ greater than $\frac{1}{2}$, and only a slight modification is needed for the more general case of non-uniform time increment integration.
## 5. ONE DIMENSIONAL IMPACT EXAMPLE

A simple one dimensional surface contact problem is solved by the classical Lagrange multiplier method of equations (6), and by the forward increment Lagrange multiplier method of equations (10). A penalty method solution is also presented for comparison.

The problem consists of two linear elastic prismatic rods moving with equal speed in opposite directions. The configuration and properties are given in Figure 2. The rods are initially undeformed and the problem is symmetric about the point at which the two rod faces impact.

The exact solution for displacement, velocity and stress at the contact surface of the left rod is given by

$$t \le t_{\text{impact}}$$

$$u_1 = v_0 t$$

$$\dot{u}_1 = v_0$$

$$\lambda = 0$$

(11a)

$$t_{\text{impact}} < t < t_{\text{release}}$$

$$u_1 = g$$

$$\dot{u}_1 = 0$$

(11b)

$$\lambda = -v_0 \sqrt{E\rho}$$

$$t \ge t_{\text{release}}$$

$$u_1 = 2v_0(g/v_0 + L\sqrt{\rho/E} - t/2)$$

$$\dot{u}_1 = -v_0$$

(11c)

$$\lambda = 0$$

where

$$t_{\text{impact}} = g/v_0$$

(11d)

and

$$t_{\text{release}} = g/v_0 + 2L\sqrt{\rho/E}$$

(11e)

which is plotted in Figures 3(a) to 3(c). The three finite element solutions also plotted in these figures are described in the following.

For the finite element analyses each rod is discretized by 20 equal length uniform strain elements and corresponding diagonal mass matrix. All three finite element solutions are obtained by direct time integration using a uniform time step increment, $h = 0.2226E-5$ sec.

The 'Lagrange-Implicit' finite element results presented in Figures 3 correspond with the classical Lagrange multiplier method given by equations (6). The integration is performed using

<figure><img src="image_3.png" /><figcaption>Figure 2. Problem description for one dimensional impact example</figcaption></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3. One dimensional impact example: (a) contact surface displacement $u_1$ versus time; (b) contact surface velocity $\dot{u}_1$ versus time; (c) contact force $\lambda$ versus time; (d) centre of mass velocity versus time</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 3 (continued)</figcaption></figure>

the implicit operator recommended by Chaudhary and Bathe,⁵ for which $β₀ = 1$ and $β₁ = ½$. The solution commences without a contact constraint and continues until the following condition is satisfied:

$$u_{n+1} + X_1 > u_{n+1} + X_2 \qquad (12a)$$

which indicates penetration. When this occurs the constraint

$$u_{n+1} + X_1 = u_{n+1} + X_2 \qquad (12b)$$

is introduced and $u_{n+1}$ is re-solved while concurrently solving for the contact force $λ_{n+1}$. The constraint remains active until the contact force becomes tensile, at which time the constraint is inactivated and $u_{n+1}$ is re-solved to allow for contact surface separation. Fortunately these two conditions for activating and inactivating the displacement constraint are never mutually exclusive. It is observed in Figures 3 that the Lagrange-Implicit results are generally poor with regard to velocity and contact force. And the results worsen if a smaller time step increment is used. Impact and release conditions can be used to improve the Lagrange-Implicit method, as proposed by Hughes et al.⁷ The more recent theoretical developments by Moreau¹¹ might also serve as a basis to formulate impact and release conditions; however, such efforts are not easily extended to the more general problem of two and three dimensional finite element surface contact. As a final comment regarding the Lagrange-Implicit method it is noted that use of the trapezoidal rule method leads to much greater inaccuracy than that observed with the implicit operator employed here.

The 'Penalty-Explicit' finite element results are obtained using a linear penalty function constraint that is active whenever inequality (12a) is satisfied. The contact force is calculated by multiplying the depth of penetration by a penalty parameter. The parameter used is $7.5E + 6$ lb/in. The transient solution is obtained using the central difference method given by equations (9). The ratio of the time integration increment to the critical time step size is 90 per cent for the unconstrained mesh. When the penalty constraint is active the numerical stability ratio increases to 99 per cent. A larger penalty parameter would reduce the contact surface penetration observed in Figure 3(a) but would also further stiffen the system and therefore require use of
a smaller time step increment to maintain numerical stability. Overall the Penalty-Explicit results are a good approximation to the exact surface contact solution.

The 'Lagrange-Explicit' finite element results are obtained by the forward increment Lagrange multiplier method with central difference integration as given in equations (10). The results obtained with this method are the most accurate in predicting the contact force and virtually exact for displacement and velocity at the contact surface during contact. It is noted that following surface contact release (i.e. $t > t_{\text{release}}$) the Lagrange-Explicit solution exhibits a noticeable error in centre of mass velocity, as observed in Figure 3(d). A closer study of the forward increment Lagrange multiplier method reveals that, during the time step in which impact occurs, the collision of surface nodal masses corresponds with the mechanics of a plastic collision and thus an associated loss of kinetic energy. However, in the mesh refinement limit this kinetic energy loss vanishes because surface nodal mass vanishes. Overall the Lagrange-Explicit results are an excellent approximation to the exact surface contact solution, and the accuracy of the method is found to improve with mesh refinement. Furthermore, the forward increment Lagrange multiplier method has no adverse affect on the numerical stability of explicit integration; the stability ratio is actually less than 90 per cent during contact.

The forward increment Lagrange multiplier method of equations (10) is outlined in the following for the one dimensional impact problem. Surface contact penetration is

$$p_{n+1}^* = G_{n+1} \{u_{n+1}^* + X\} = (u_{n+1}^* + X_2) - (u_{n+1}^* + X_1) \quad (13a)$$

$$\text{in which} \qquad G_{n+1} = [-1, 1] \qquad (13b)$$

If $p_{n+1}^*$ is negative, then penetration exists and the contact force is calculated

$$\lambda_n = [h^2 G_{n+1} M^{-1} G_{n+1}^T]^{-1} p_{n+1}^* = \frac{p_{n+1}^*}{h^2(1/m_1 + 1/m_2)} \quad (13c)$$

in which

$$M = \begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix} \qquad (13d)$$

The displacements due to this contact force are then

$$\left\{ \begin{matrix} u_{n+1}^c \\ u_{n+1}^2 \end{matrix} \right\} = -h^2 M^{-1} G_{n+1}^T \lambda_n = h^2 \lambda_n \begin{Bmatrix} \frac{1}{m_1} \\ \frac{-1}{m_2} \end{Bmatrix} \quad (13e)$$

Given the symmetries, $m_2 = m_1$, $X_2 = X_1 + 2g$, and $u_{n+1}^* = -u_{n+1}^*$, leads to

$$p_{n+1}^* = 2(-u_{n+1}^* + g) \quad (14a)$$

and if $p_{n+1}^*$ is negative, then

$$\lambda_n = \frac{m_1 p_{n+1}^*}{2h^2} \quad (14b)$$

$$u_{n+1}^c = \frac{1}{2} p_{n+1}^* \quad (14c)$$

$$u_{n+1}^1 = g \quad (14d)$$

Note that the displacement in equation (14d) corresponds with the exact displacement solution of equation (11b).
## 6. TWO DIMENSIONAL SURFACE CONTACT KINEMATIC CONDITIONS

The explicit displacement update given by equation (10c) allows contactor nodes to penetrate target surfaces. This is illustrated in Figure 4 for a two dimensional target surface spanning between target nodes $A_I$ and $B_I$ which is penetrated by the contactor node $C_I$ during the time step from $t_n$ to $t_{n+1}^*$. The penetration is eliminated by requiring nodes $A_I$, $B_I$ and $C_I$ to lie on a common line at time $t_{n+1}$. This condition may be expressed as

$$ \mathbf{x}_{n+1} = (1 - \alpha_{n+1}) \mathbf{x}_{n+1} + \alpha_{n+1} \mathbf{x}_{n+1} \quad (15) $$

where $\alpha_{n+1}^*$ is a dimensionless constraint parameter and $\mathbf{x}_{n+1}^*$, $\mathbf{x}_{n+1}^*$ and $\mathbf{x}_{n+1}^*$ are the spatial co-ordinates of nodes $A_I$, $B_I$ and $C_I$ at time $t_{n+1}$. The constraint element matrix corresponding with equation (15) is

$$ \bar{\mathbf{G}}_{n+1} = \bar{\mathbf{x}}_{n+1} = \mathbf{0} \quad (16a) $$

where

$$ \bar{\mathbf{G}}_{n+1} = \begin{bmatrix} (1 - \alpha_{n+1}^*) & 0 & \alpha_{n+1}^* & 0 & -1 & 0 \\ 0 & (1 - \alpha_{n+1}^*) & 0 & \alpha_{n+1}^* & 0 & -1 \end{bmatrix} \quad (16b) $$

and

$$ \bar{\mathbf{x}}_{n+1} = \{ \mathbf{x}_{n+1}^*, \mathbf{y}_{n+1}^*, \mathbf{x}_{n+1}^*, \mathbf{y}_{n+1}^*, \mathbf{x}_{n+1}^*, \mathbf{y}_{n+1}^* \}^T \quad (16c) $$

The relative position of node $C_I$ along the target surface is associated with the value of $\alpha_{n+1}^*$. In the following $\alpha_{n+1}^*$ is prescribed to maintain the relative position of node $C_I$, thereby imposing a tangential constraint to control sliding. If the node $C_I$ is in contact with the target surface at time $t_n$, then $\alpha_{n+1}^*$ is simply

$$ \alpha_{n+1}^* = \frac{\Delta_{TI}}{|\mathbf{V}_{ABI}|} \quad (17a) $$

where

$$ \Delta_{TI} = \frac{\mathbf{V}_{ACI} \cdot \mathbf{V}_{ABI}}{|\mathbf{V}_{ABI}|} \quad (17b) $$

$$ \mathbf{V}_{ABI} = \mathbf{x}_{n+1}^* - \mathbf{x}_{n}^* \quad (17c) $$

$$ \mathbf{V}_{ACI} = \mathbf{x}_{n+1}^* - \mathbf{x}_{n}^* \quad (17d) $$

The above geometric relations are illustrated in Figure 5.

<figure><img src="image_8.png" /><figcaption>Figure 4. Illustration of 1th constraint element: (a) configuration at time $t_n$; (b) configuration at time $t_{n+1}^*$</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 5. Relative contactor positions: (a) configuration at time $t_n$; (b) configuration at time $t_{n+1}$</figcaption></figure>

The distance between the node $C_I$ and the target surface at time $t_n$, measured in the direction normal to the target surface, is

$$ \Delta_{NI} = \frac{V_{ACIx} V_{ABIy} - V_{ACIy} V_{ABIx}}{|V_{ABI}|} \quad (18) $$

If $\Delta_{NI}$ is greater than zero, then node $C_I$ is not in contact with the target surface at time $t_n$. In this case the relative position of node $C_I$ and the associated value for $\alpha_{I,n+1}$ should correspond with the point at which contact occurs during the time increment. However, this is quite difficult to determine precisely. Therefore an approximation for the point of contact is based on the assumption that nodal displacements are linear with respect to time during the increment $t_n$ to $t_{n+1}^*$. This leads to

$$ \alpha_{I,n+1} = \frac{a\Delta_{TI} + a^*\Delta_{TI}^*}{a|V_{ABI}| + a^*|V_{ABI}^*|} \quad (19a) $$

where

$$ a = \frac{\Delta_{NI}^*}{\Delta_{NI}^* - \Delta_{NI}} \quad (19b) $$

$$ a^* = 1 - a \quad (19c) $$

and

$$ \Delta_{TI}^* = \frac{V_{ACI}^* \cdot V_{ABI}^*}{|V_{ABI}^*|} \quad (19d) $$

$$ \Delta_{NI}^* = \frac{V_{ACIx}^* V_{ABIy}^* - V_{ACIy}^* V_{ABIx}^*}{|V_{ABI}^*|} \quad (19e) $$

$$ V_{ABI}^* = x_{n+1}^* - x_{n+1}^* \quad (19f) $$

$$ V_{ACI}^* = x_{n+1}^* - x_{n+1}^* \quad (19g) $$

The surface contact displacement constraint form given by equations (16) is typical and can be found for example in References 1, 4 and 8. The usual implementations are with low order finite elements that employ linear spatial and displacement interpolation functions along target and contactor bodies. With regard to the evaluation of $\alpha_{I,n+1}$, it might be argued that equations (19) provide a sufficiently accurate approximation for the point of contact even when the contactor node is not in contact with the target surface at time $t_n$. Nevertheless, equations (19) provide
a more accurate approximation of the point of contact without significant computational cost or complexity. At the other extreme, attempting to calculate the exact point of contact occurring between time $t_n$ and $t_{n+1}$ is generally quite complicated, and may provide no real improvement in accuracy considering the limitations of the finite element and time integration discretizations.

## 7. GAUSS-SEIDEL ITERATION

An equation solving strategy based on the Gauss-Seidel method is developed to solve equations (10d) and (10e). In the following it is assumed that the 2 × 6 constraint element matrices, given by equation (16b), are assembled to form the global matrix $G_{n+1}$. The two rows of $G_{n+1}$ corresponding with $\bar{G}_{n+1}$ are given by $G_I$. Equation (10d) is then expressed as

for $I = 1, 2, \dots, N$

$$ \sum_{J=1}^{N} h^2 G_{n+1} M^{-1} G_{n+1}^{T} \lambda_{J} = G_{n+1} \{\mathbf{u}_{n+1}^{*} + \mathbf{X}\} \quad (20) $$

where $M$ is the global mass matrix, $N$ is the number of contactor nodes, and $\lambda_J$ is the 2 × 1 vector consisting of surface contact force components $\lambda_{xJ}$ and $\lambda_{yJ}$ associated with the $J$th contactor node. The $i$ to $i+1$ Gauss-Seidel iteration of equation (20) is,

for $I = 1, 2, \dots, N$

$$ \begin{aligned} {}^{i, i+1} \mathbf{p}_I = & G_{n+1} \{\mathbf{u}_{n+1}^* + \mathbf{X}\} - \sum_{J=1}^{N} h^2 G_{n+1} M^{-1} G_{n+1}^{T} {}^i \lambda_J \\ & - \sum_{J=1}^{I-1} h^2 G_{n+1} M^{-1} G_{n+1}^{T} {}^{i+1} \lambda_J \end{aligned} \quad (21a) $$

$$ \Delta \lambda_I = [h^2 G_{n+1} M^{-1} G_{n+1}^T]^{-1} {}^{i, i+1} \mathbf{p}_I \quad (21b) $$

$$ {}^{i+1} \lambda_{n} = {}^{i} \lambda_{I} + \Delta \lambda_{I} \quad (21c) $$

Equations (21b) and (21c) may be efficiently evaluated by element level calculations. However, equation (21a) is complicated because coupling generally occurs between the $I$th and $J$th displacement constraints. Fortunately a more efficient expression for calculating ${}^{i, i+1} \mathbf{p}_I$ can be obtained by combining equations (10d) and (10e). From equation (10e), the value of $\mathbf{u}_{n+1}^c$ at the beginning of iteration $i$ to $i+1$ is

$$ {}^{i} \mathbf{u}_{n+1}^{c} = - \sum_{J=1}^{N} h^2 M^{-1} G_{n+1}^{T} {}^{i} \lambda_{J} \quad (22a) $$

and at the end of the iteration

$$ {}^{i+1} \mathbf{u}_{n+1}^{c} = - \sum_{J=1}^{N} h^2 M^{-1} G_{n+1}^{T} {}^{i+1} \lambda_{J} \quad (22b) $$

Defining the difference between ${}^{i+1} \mathbf{u}_{n+1}^c$ and ${}^i \mathbf{u}_{n+1}^c$ as

$$ \Delta \mathbf{u}^c = {}^{i+1} \mathbf{u}_{n+1}^c - {}^i \mathbf{u}_{n+1}^c \quad (23) $$

and substituting from equations (21c) and (22) leads to

$$ \Delta \mathbf{u}^c = \sum_{J=1}^{N} \Delta \mathbf{u}_J^c \quad (24a) $$
where

$$ \Delta \mathbf{u}_j^c = - h^2 \mathbf{M}^{-1} \mathbf{G}_{n+1}^T \Delta \lambda_j \quad (24b) $$

During iteration $i$ to $i + 1$, just prior to solving for $\Delta\lambda_I$, the current estimate of $\mathbf{u}_{n+1}^c$, denoted as $i, i+1 \mathbf{u}_{n+1}^c$, is

$$ i, i+1 \mathbf{u}_{n+1}^c = i \mathbf{u}_{n+1}^c + \sum_{J=1}^{I-1} \Delta \mathbf{u}_J^c \quad (25) $$

Substituting from equations (22a) and (24) into equation (25), and also substituting the difference, $(i+1\lambda_n - i\lambda_n)$, in place of $\Delta\lambda_J$, leads to

$$ i, i+1 \mathbf{u}_{n+1}^c = - \sum_{J=1}^{N} h^2 \mathbf{M}^{-1} \mathbf{G}_{n+1}^T \lambda_n - \sum_{J=1}^{I-1} h^2 \mathbf{M}^{-1} \mathbf{G}_{n+1}^T \lambda_n \quad (26) $$

Substituting from equation (26) into equation (21a), and combining with (21b) yields

$$ \Delta \lambda_I = [h^2 \mathbf{G}_{n+1} \mathbf{M}^{-1} \mathbf{G}_{n+1}^T]^{-1} \mathbf{G}_{n+1} \{i, i+1 \mathbf{u}_{n+1}^c + \mathbf{u}_{n+1}^* + \mathbf{X}\} \quad (27) $$

Following calculation of $\Delta\lambda_I$ from equation (27), the change in $i, i+1 \mathbf{u}_{n+1}^c$, according to equation (24b), is given by

$$ \Delta \mathbf{u}_I^c = - h^2 \mathbf{M}^{-1} \mathbf{G}_{n+1}^T \Delta \lambda_I \quad (28) $$

Based on equations (27) and (28), the $i$ to $i + 1$ Gauss-Seidel iteration can be expressed as, for $I = 1, 2, \dots, N$

$$ i, i+1 \mathbf{p}_I = \mathbf{G}_{n+1} \{i, i+1 \mathbf{u}_{n+1}^c + \mathbf{u}_{n+1}^* + \mathbf{X}\} \quad (29a) $$

$$ \Delta \lambda_I = [h^2 \mathbf{G}_{n+1} \mathbf{M}^{-1} \mathbf{G}_{n+1}^T]^{-1} i, i+1 \mathbf{p}_I \quad (29b) $$

$$ i+1 \lambda_n = i \lambda_n + \Delta \lambda_I \quad (29c) $$

$$ \Delta \mathbf{u}_I^c = - h^2 \mathbf{M}^{-1} \mathbf{G}_{n+1}^T \Delta \lambda_I \quad (29d) $$

$$ i, i+1 \mathbf{u}_{n+1}^c \leftarrow \{i, i+1 \mathbf{u}_{n+1}^c + \Delta \mathbf{u}_I^c\} \quad (29e) $$

where the arrow in equation (29e) means: replace $i, i+1 \mathbf{u}_{n+1}^c$ by $\{i, i+1 \mathbf{u}_{n+1}^c + \Delta \mathbf{u}_I^c\}$. In equations (29) all references are to the $I$th constraint element only. Therefore the Gauss-Seidel iteration can proceed by simple element level calculations, without assembling any global matrices. The coupling between the $I$th and $J$th constraint elements that appears in equations (21a) is accounted for very efficiently by the global vector update in equation (29e). Also, the matrix in equation (29b) is a simple $2 \times 2$ diagonal matrix given by

$$ \mathbf{G}_{n+1} \mathbf{M}^{-1} \mathbf{G}_{n+1}^T = \left( (1 - \alpha_{n+1})^2 / m_{AI} + \alpha_{n+1}^2 / m_{BI} + 1 / m_{CI} \right) \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \quad (30) $$

where $m_{AI}$, $m_{BI}$ and $m_{CI}$ are the nodal masses of nodes $A_I$, $B_I$ and $C_I$ for the $I$th constraint element. The Gauss-Seidel algorithm for a two dimensional surface contact based on equations (29) is given in the Appendix as Algorithm 2. This algorithm corresponds with step 2 of Algorithm 1.

In addition to the steps outlined above, an iterative convergence check is necessary. The following is recommended and appears in step 2.3 of Algorithm 2

$$ |\Delta \lambda| < \text{tol}|^{i+1} \lambda_n| \quad (31) $$
where $|\Delta\lambda|$ is the mean-square norm of contact force changes, and $|{}^{i+1}\lambda_n|$ is the mean-square norm of the $i+1$ estimate of $\lambda_n$.

## 8. TWO DIMENSIONAL SURFACE CONTACT FORCE CONDITIONS

The procedure developed in Section 7 for solving equations (10d) and (10e) leads to surface contact force reactions $\lambda_n$ that satisfy the prescribed displacement constraints. However, satisfying the displacement constraints may result in contact force reactions that fail to satisfy prescribed surface contact force conditions. When this occurs the associated displacement constraints are relaxed to obtain a modified solution that satisfies the contact force conditions.

In the global co-ordinate system the components of $\lambda_i$ are $\lambda_{xi}$ and $\lambda_{yi}$. These components may be transformed to obtain the components $\lambda_{TI}$ and $\lambda_{NI}$, which correspond with the directions tangent and normal to the $I$th target surface at time $t_{n+1}$. An elementary set of surface contact force conditions applicable to $\lambda_{TI}$ and $\lambda_{NI}$ are,

for $I = 1, 2, \dots, N$

$$ \lambda_{NI} \le 0 \qquad (32a) $$

$$ |\lambda_{TI}| \le \mu |\lambda_{NI}| \qquad (32b) $$

Inequality (32a) is a tension limit force condition, while inequality (32b) is a friction limit force condition in which $\mu$ is a friction coefficient.

The Gauss-Seidel strategy is easily modified to account for inequalities (32) during iteration. The necessary modifications are given in the Appendix in Algorithm 3, which is a sequence of steps to follow step 2.2.3 of Algorithm 2. The procedure in Algorithm 3 begins in steps 2.2.3a and 2.2.3b by calculating the $i+1$ iterative estimate of the $\lambda_{TI}$ and $\lambda_{NI}$. In step 2.2.3c ${}^{i+1}\lambda_{NI}$ is set to zero if it is tensile, and ${}^{i+1}\lambda_{TI}$ is reduced in magnitude if it exceeds the friction limit. Following the modification of ${}^{i+1}\lambda_{TI}$ and ${}^{i+1}\lambda_{NI}$ in step 2.2.3c, the corresponding modification to $\Delta\lambda_{xI}$ and $\Delta\lambda_{yI}$ is calculated in step 2.2.3d before proceeding to step 2.2.4 of Algorithm 2.

In the above strategy the displacement constraint element matrices $\tilde{G}_{I}$ are not altered during iteration. However, enforcement of these constraints is relaxed by the effects of step 2.2.3c. For example, if $\lambda_{NI}$ is zero because of step 2.2.3c, then the $I$th contactor node may separate from the target at time $t_{n+1}$. Similarly, sliding occurs when $\lambda_{TI}$ is reduced in magnitude as a result of the friction limit condition in step 2.2.3c. Since the displacement constraints are re-evaluated by equations (17) and (19) at the beginning of each time increment, the constraint parameters $\alpha_{I}$ will change at each time step when sliding occurs.

## 9. TWO DIMENSIONAL IMPACT EXAMPLE

A two dimensional surface contact example is presented to demonstrate the application of Algorithms 2 and 3. The example consists of a contactor surface impacting a stationary target surface during an increment of motion from time $t_n$ to $t_{n+1}$.

The target and contactor surfaces are illustrated in Figure 6. The target surface includes three segments spanning from nodes 1 to 4, which are motionless at $t_n$. Nodes 5 to 16 are contactor surface nodes, which are in motion at $t_n$. The contactor nodes penetrate the target segments during the increment from $t_n$ to $t_{n+1}^*$, as shown in Figure 6(a).

The spatial co-ordinate positions of target and contactor nodes at $t_n$ and $t_{n+1}^*$ are given in Table I. The information in Table I corresponds with completion of step 1 in Algorithm 1. Step 2 in Algorithm 1 is replaced by Algorithms 2 and 3 to solve for $\lambda_n$ and $u_{n+1}$.
<figure><img src="image_10.png" /><figcaption>Figure 6. Two dimensional impact example: (a) contactor node positions denoted by circles at $t_n$ and asterisks at $t_{n+1}^*$; (b) node positions at $t_{n+1}$ for $\mu > 0.88$; no sliding; (c) node positions at $t_{n+1}$ for $\mu = 0$; no resistance to sliding</figcaption></figure>

<table><caption>Table I. Surface node spatial co-ordinates at time $t_n$ and $t_{n+1}^*$</caption><thead><tr><th rowspan="2">Node</th><th colspan="2">Co-ordinates at t<sub>n</sub></th><th colspan="2">Co-ordinates at t<sub>n+1</sub><sup>*</sup></th><th rowspan="2">Nodal mass</th></tr><tr><th>x<sub>n</sub></th><th>y<sub>n</sub></th><th>x<sub>n+1</sub><sup>*</sup></th><th>y<sub>n+1</sub><sup>*</sup></th></tr></thead><tbody><tr><td>1</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>2.000</td></tr><tr><td>2</td><td>10.000</td><td>0.000</td><td>10.000</td><td>0.000</td><td>4.000</td></tr><tr><td>3</td><td>20.000</td><td>0.000</td><td>20.000</td><td>0.000</td><td>4.000</td></tr><tr><td>4</td><td>30.000</td><td>0.000</td><td>30.000</td><td>0.000</td><td>2.000</td></tr><tr><td>5</td><td>0.000</td><td>-2.000</td><td>2.000</td><td>2.000</td><td>1.000</td></tr><tr><td>6</td><td>2.000</td><td>-2.000</td><td>3.000</td><td>0.500</td><td>1.000</td></tr><tr><td>7</td><td>4.000</td><td>-2.000</td><td>5.000</td><td>3.000</td><td>1.000</td></tr><tr><td>8</td><td>6.000</td><td>-2.000</td><td>8.000</td><td>2.000</td><td>1.000</td></tr><tr><td>9</td><td>12.000</td><td>-2.000</td><td>11.000</td><td>2.000</td><td>1.000</td></tr><tr><td>10</td><td>13.000</td><td>-2.000</td><td>15.000</td><td>3.000</td><td>1.000</td></tr><tr><td>11</td><td>16.000</td><td>-2.000</td><td>17.000</td><td>3.000</td><td>1.000</td></tr><tr><td>12</td><td>19.000</td><td>-2.000</td><td>19.000</td><td>0.500</td><td>1.000</td></tr><tr><td>13</td><td>22.000</td><td>-2.000</td><td>21.000</td><td>3.000</td><td>1.000</td></tr><tr><td>14</td><td>24.000</td><td>-2.000</td><td>25.000</td><td>3.000</td><td>1.000</td></tr><tr><td>15</td><td>26.000</td><td>-2.000</td><td>28.000</td><td>3.000</td><td>1.000</td></tr><tr><td>16</td><td>27.000</td><td>-2.000</td><td>30.000</td><td>3.000</td><td>1.000</td></tr></tbody></table>

A constraint element is activated for each contactor node that penetrates a target segment during the increment from $t_n$ to $t_{n+1}^*$. In this example 12 constraint elements are activated, corresponding with contactor nodes 5 to 16. The constraint element connectivity and constraint parameters $\alpha_{n+1}^*$, evaluated using equations (19), are given in Table II.

Two different friction coefficient cases are considered. In the first case the friction coefficient $\mu$ is large enough to prevent all contactor nodes from sliding on the target segments during the time
<table><caption>Table II. Constraint element connectivity and parameter values</caption><thead><tr><th>Element-I</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><th>Node-A<sub>I</sub></th><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>2</td><td>3</td><td>3</td><td>3</td><td>3</td></tr><tr><th>Node-B<sub>I</sub></th><td>2</td><td>2</td><td>2</td><td>2</td><td>3</td><td>3</td><td>3</td><td>3</td><td>4</td><td>4</td><td>4</td><td>4</td></tr><tr><th>Node-C<sub>I</sub></th><td>5</td><td>6</td><td>7</td><td>8</td><td>9</td><td>10</td><td>11</td><td>12</td><td>13</td><td>14</td><td>15</td><td>16</td></tr><tr><th>α<sub>I</sub><sub>n+1</sub></th><td>0.100</td><td>0.280</td><td>0.440</td><td>0.700</td><td>0.150</td><td>0.380</td><td>0.640</td><td>0.900</td><td>0.160</td><td>0.440</td><td>0.680</td><td>0.820</td></tr></tbody></table>

increment. A value of $\mu$ greater than 0.88 leads to this result. In the second case the friction coefficient is zero. The converged results for these two cases are illustrated in Figures 6(b) and 6(c), with numerical results for surface contact forces and spatial co-ordinates given in Tables III and IV. The time step increment for the calculations is $h = 1 \cdot 0$.

An overview of the convergence behaviour of Algorithms 2 and 3 is given in Figure 7. Mean-square norms are used to measure relative contact force error. The actual error is

$$ \text{actual error} = \frac{|^{i+1}\lambda_n - \lambda_n|}{|\lambda_n|} \quad (33) $$

where $^{i+1}\lambda_n$ is the contact force estimate at completion of the current iteration, and $\lambda_n$ is the converged solution. The estimated error is

$$ \text{estimated error} = \frac{|^{i+1}\lambda_n - {}^{i}\lambda_n|}{|^{i+1}\lambda_n|} \quad (34) $$

where ${}^i\lambda_n$ is the estimate from the previous iteration. Equation (34) is similar to the expression used for the convergence check in step 2.3 of Algorithm 2.

It is observed in Figure 7 that the order of convergence is approximately linear. Nevertheless, an accurate solution is obtained quickly. For example, the actual contact force error is reduced from 100 to 1 per cent in only 3 iterations. Furthermore, the relaxation of tangential displacement constraints in step 2.2.3c of Algorithm 3, which occurs in the frictionless sliding case, has only a slight effect on the convergence behaviour.

Note that contactor nodes 6 and 12 are not in contact with the target at time $t_{n+1}$. Enforcing the displacement constraints initially assumed for these contactor nodes results in tensile normal force reactions. This feature of the example serves to exercise the relaxation of normal displacement constraints in step 2.2.3c of Algorithm 3.

In this example the incremental displacements are large and the pattern of contactor node motions is unusually chaotic. For smaller incremental displacements and a more orderly pattern of motion, Algorithms 2 and 3 convergence more quickly. Also, increasing the number of target segments and contactor nodes has been found to have little effect on convergence behaviour.

A significant factor affecting convergence is the relative magnitude of contactor and target nodal masses. If target node masses are significantly greater than contactor node masses, the Gauss-Seidel method will converge more rapidly. Conversely, if the target nodes are relatively less massive, the convergence rate is slower. This occurs because equation (10d) becomes more diagonally dominant as the magnitude of target masses increases. This effect should be considered when assigning contactor and target surface designations.

## 10. TWO DIMENSIONAL EXAMPLE

As a further demonstration of the forward increment Lagrange multiplier method formulated in Algorithms 2 and 3, two finite element example problems are presented. The finite element
<table><caption>Table III. Forward increment Lagrange multiplier solution for surface contact forces $\lambda_n$. No sliding for $\mu > 0.88$, no resistance to sliding for $\mu = 0$</caption><thead><tr><th rowspan="2">Contact element I</th><th rowspan="2">Contact node</th><th colspan="2">μ > 0.88</th><th colspan="2">μ = 0</th></tr><tr><th>λ<sub>TI</sub><sub>n</sub></th><th>λ<sub>NI</sub><sub>n</sub></th><th>λ<sub>TI</sub><sub>n</sub></th><th>λ<sub>NI</sub><sub>n</sub></th></tr></thead><tbody><tr><td>1</td><td>5</td><td>-0.5773</td><td>-0.9155</td><td>0.0000</td><td>-0.9155</td></tr><tr><td>2</td><td>6</td><td>0.0000</td><td>0.0000</td><td>0.0000</td><td>0.0000</td></tr><tr><td>3</td><td>7</td><td>-0.2486</td><td>-1.9289</td><td>0.0000</td><td>-1.9298</td></tr><tr><td>4</td><td>8</td><td>-0.7111</td><td>-0.9423</td><td>0.0000</td><td>-0.9441</td></tr><tr><td>5</td><td>9</td><td>0.6817</td><td>-0.9349</td><td>0.0000</td><td>-0.9316</td></tr><tr><td>6</td><td>10</td><td>-1.0628</td><td>-1.8470</td><td>0.0000</td><td>-1.8454</td></tr><tr><td>7</td><td>11</td><td>-0.4890</td><td>-1.8032</td><td>0.0000</td><td>-1.8025</td></tr><tr><td>8</td><td>12</td><td>0.0000</td><td>0.0000</td><td>0.0000</td><td>0.0000</td></tr><tr><td>9</td><td>13</td><td>0.7551</td><td>-1.7067</td><td>0.0000</td><td>-1.7060</td></tr><tr><td>10</td><td>14</td><td>-0.2748</td><td>-1.5813</td><td>0.0000</td><td>-1.5809</td></tr><tr><td>11</td><td>15</td><td>-0.7292</td><td>-1.4872</td><td>0.0000</td><td>-1.4871</td></tr><tr><td>12</td><td>16</td><td>-1.2441</td><td>-1.4244</td><td>0.0000</td><td>-1.4246</td></tr></tbody></table>

<table><caption>Table IV. Forward increment Lagrange multiplier solution for spatial co-ordinates at time $t_{n+1}$. No sliding for $\mu > 0$, no resistance to sliding for $\mu = 0$</caption><thead><tr><th rowspan="2">Node</th><th colspan="2">μ > 0.88</th><th colspan="2">μ = 0</th></tr><tr><th>x<sub>n+1</sub></th><th>y<sub>n+1</sub></th><th>x<sub>n+1</sub></th><th>y<sub>n+1</sub></th></tr></thead><tbody><tr><td>1</td><td>0.4410</td><td>1.0915</td><td>0.0052</td><td>1.0939</td></tr><tr><td>2</td><td>10.2177</td><td>1.0477</td><td>9.9880</td><td>1.0465</td></tr><tr><td>3</td><td>20.1130</td><td>1.2646</td><td>19.9654</td><td>1.2608</td></tr><tr><td>4</td><td>30.7083</td><td>1.5971</td><td>29.9508</td><td>1.5732</td></tr><tr><td>5</td><td>1.4186</td><td>1.0871</td><td>1.9957</td><td>1.0845</td></tr><tr><td>6</td><td>3.0000</td><td>0.5000</td><td>3.0000</td><td>0.5000</td></tr><tr><td>7</td><td>4.7427</td><td>1.0722</td><td>4.9908</td><td>1.0702</td></tr><tr><td>8</td><td>7.2847</td><td>1.0608</td><td>7.9955</td><td>1.0559</td></tr><tr><td>9</td><td>11.7020</td><td>1.0803</td><td>11.0200</td><td>1.0686</td></tr><tr><td>10</td><td>13.9779</td><td>1.1301</td><td>15.0396</td><td>1.1550</td></tr><tr><td>11</td><td>16.5507</td><td>1.1865</td><td>17.0387</td><td>1.1979</td></tr><tr><td>12</td><td>19.0000</td><td>0.5000</td><td>19.0000</td><td>0.5000</td></tr><tr><td>13</td><td>21.8082</td><td>1.3178</td><td>21.0533</td><td>1.2949</td></tr><tr><td>14</td><td>24.7749</td><td>1.4109</td><td>25.0494</td><td>1.4199</td></tr><tr><td>15</td><td>27.3178</td><td>1.4907</td><td>28.0465</td><td>1.5136</td></tr><tr><td>16</td><td>28.8012</td><td>1.5372</td><td>30.0445</td><td>1.5761</td></tr></tbody></table>

problems are compared with an elementary single degree of freedom (SDOF) sliding contact problem, which is illustrated in Figure 8(a). The SDOF problem consists of a rigid object that slides horizontally on a rigid base, while connected to a fixed point by a linear elastic spring. The mass of the object is specified as 30, the spring stiffness is 10 000 and the undeformed spring length is 10. At time $t_0$ the object has an initial horizontal velocity of 75, and is located directly above the fixed point such that the spring is undeformed.
<figure><img src="image_11.png" /><figcaption>Figure 7. Convergence behaviour for two dimensional impact example</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 8. Two dimensional sliding contact problem: (a) single degree of freedom object sliding on rigid base; (b) finite element object sliding on rigid base; (c) finite element object sliding on finite element base</figcaption></figure>

The first of the two finite element problems is illustrated in Figure 8(b), which features a two dimensional linear elastic object. The object is spatially discretized by 12 bilinear quadrilateral plane stress elements, which are initially square at time $t_0$. The modulus of elasticity of the object is $10^6$, Poisson's ratio is 0.25 and the mass density is 2. The spring is modelled by a two node uniform strain truss element having a modulus of elasticity of $10^5$, a cross section area of 1 and density of 1·2. Large rotation of the truss is accounted for by incrementally updating spatial co-ordinates during the finite element internal force calculations.
The second finite element problem is illustrated in Figure 8(c), where an elastic base is introduced and discretized by 15 bilinear quadrilateral plane stress elements. The base elements are initially square at time $t_0$. The material properties and mass density of the base are the same as given above for the elastic object.

For the SDOF system, the equation of motion at time $t_n$ is

$$m\ddot{u}_n - F_n + \lambda_n = 0 \quad (35a)$$

where $m$ is the mass of the sliding object, $u$ is the horizontal displacement of the object and $\lambda_n$ is the contact friction force. The horizontal and vertical components of the spring force are

$$F_n = F_n\left(\frac{u_n}{R_n}\right) \quad (35b)$$

$$F_n = F_n\left(\frac{L}{R_n}\right) \quad (35c)$$

where

$$R_n = \sqrt{L^2 + u_n^2} \quad (35d)$$

and

$$F_n = k(R_n - L) \quad (35e)$$

in which the undeformed spring length is $L$ and the spring stiffness is $k$. Since the SDOF object is constrained from moving vertically, the vertical contact force is given by $-F_y$. Equations (35) constitute a non-linear initial value problem that is readily solved by direct time integration using the central difference operators given in equations (9).

An incremental relation for the friction force $\lambda_x$ is obtained by substituting equation (9b) into (35a) and rearranging so that

$$\lambda_n = -\frac{m}{h^2}(u_{n+1} - 2u_n + u_{n-1}) - F_n \quad (36)$$

This force is defined as a 'stick' force if it prevents sliding (i.e. $u_{n+1} = u_n$) and is given by

$$\lambda_{\text{stick}} = \frac{m}{h^2}(u_n - u_{n+1}) - F_n \quad (37a)$$

Observing that the vertical force given by equation (35c) is non-negative and introducing a coefficient of friction, $\mu$, leads to the following expression for the friction force:

$$\lambda_n = \begin{cases} \lambda_{\text{stick}}, & |\lambda_{\text{stick}}| \le \mu F_n \\ \mu F_n \text{sign}(\lambda_{\text{stick}}), & |\lambda_{\text{stick}}| > \mu F_n \end{cases} \quad (37b)$$

The central difference integration solution to equations (35) and (37) is designated as 'SDOF' in Figures 9 and 10. The horizontal displacement and vertical contact force results presented in Figures 9(a) and 9(b) correspond with frictionless contact, $\mu = 0$. Results for a friction coefficient of $\mu = 0.5$ are given in Figures 10. The SDOF solution is calculated using a uniform time step increment of $h = 0.000850$, which is two orders of magnitude less than the critical time step size.

For comparison to the SDOF solutions, finite element results are also presented in Figures 9 and 10. The problem description of Figure 8(b) is referred to as 'FE object', while that of Figure 8(c) is referred to as 'FE object/base'. The finite element solutions are obtained using the Beta-2 explicit integration method with Algorithms 2 and 3.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 9. Two dimensional sliding example. Frictionless contact, $\mu = 0$: (a) centre of mass horizontal displacement of sliding object versus time; (b) vertical contact force resultant versus time</figcaption></figure>
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 10. Two dimensional sliding example. Friction contact, $\mu = 0.5$: (a) centre of mass horizontal displacement of sliding object versus time; (b) vertical contact force resultant versus time; (c) horizontal contact force resultant versus time</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 10 (continued)</figcaption></figure>

The forward increment Lagrange multiplier method with Beta-2 explicit integration is similar to the multi-step central difference integration in Algorithm 1. The modification of Algorithm 2 for Beta-2 explicit integration is accomplished by simply replacing $h^2$, by $h^2(\beta_1 + 1/2)$, in steps 2.2.3 and 2.2.5. The Beta-2 method is preferred here because it allows for high frequency numerical damping, which improves the finite element results. The Beta-2 integration parameters used are: $\beta_0 = 0$ and $\beta_1 = 1$. The FE object and FE object/base problems are solved by the Beta-2 method using a uniform time step increment of $h = 0.000850$, which is 85 per cent of critical for both finite element problems.

The contact surface of the finite element object is designated as the contactor surface, which consists of seven contactor nodes. The vertical contact force resultant given in Figures 9(b) and 10(b) corresponds with

$$ \text{vertical contact force resultant} = \sum_{I=1}^{7} \lambda_{yI} \qquad (38a) $$

where $\lambda_{yI}$ are calculated by Algorithms 2 and 3. Similarly the horizontal contact force resultant for the friction case in Figure 10(c) is

$$ \text{horizontal contact force resultant} = \sum_{I=1}^{7} \lambda_{xI} \qquad (38b) $$

The finite element object sliding on the rigid base is a special case in which target surface mass is assumed to be infinite. Algorithm 2 is easily modified to accommodate this assumption by setting the terms involving the reciprocal of target nodal masses to zero in step 2.2.3. In this
special case the matrix in equations (10d) is diagonal. Consequently, the exact solution for $\lambda_n$ is obtained by one iteration, at each time step.

By comparison, the problem in Figure 8(c) with the finite element base is solved by the more general methodology of Algorithms 2 and 3. The convergence tolerance in step 2.3 of Algorithm 2 is specified as tol = 0.001. Satisfying this rather strict tolerance requires not more than four iterations at any time step in the transient analyses, and the number of iterations required is nearly independent of the friction coefficient.

The horizontal centre of mass displacement results for the FE object and FE object/base frictionless contact are quite similar, hence the dotted line in Figure 9(a) is visually obscured by the dashed line. It is observed in Figures 9(a) and 10(a) that the horizontal centre of mass motion of the elastic finite element object is nearly identical to the motion of the SDOF object in all cases. This is true even when friction, $\mu = 0.5$, is prescribed. The high frequency oscillations in the contact force resultants are expected as a consequence of the elastic compliance of the finite element object and the rapid changes in friction forces acting on the object. This is observed in Figure 10(c), where the horizontal contact force resultant is large and suddenly changes direction.

The sliding contact problem presented here is a demanding test of the forward increment Lagrange multiplier method and the modified Gauss-Seidel solution strategy in Algorithms 2 and 3. The forward increment Lagrange multiplier method proves to be highly accurate for enforcing the coupled normal and tangential contact constraints and friction limit force condition. Furthermore, the modified Gauss-Seidel strategy converges rapidly.

## 11. CONCLUSION

It should be emphasized that the coupled system of equations involved in the forward increment Lagrange multiplier method is typically small. For two dimensional surface contact, the number of unknowns is simply equal to two times the number of active contactor nodes. Moreover, the incremental displacements of interior nodes and surface nodes away from contact regions are calculated by the usual explicit operator update without loss of efficiency. The computational cost per increment for these degrees of freedom is not affected by the forward increment Lagrange multiplier contact forces.

The modified Gauss-Seidel equation solving strategy formulated in this paper requires little computational effort per iteration, and converges quickly for the forward increment Lagrange multiplier two dimensional surface contact problems. Also, it is quite convenient that non-linear surface contact force conditions can be enforced concurrently in the modified Gauss-Seidel formulation. This combined methodology generally proves to be a very accurate and efficient alternative to the penalty function method. Since the forward increment Lagrange multiplier method does not adversely affect numerical stability, a slightly greater computational cost per increment with the method may be more than compensated in many cases by the advantage of a larger allowable time step size.

For the sake of simplicity, no distinction between static and dynamic friction coefficients is made in this paper. A full range of friction law and contact condition refinements remains a topic for further investigation.

## ACKNOWLEDGEMENT

This work was carried out under a TRW Independent Research and Development Program, project number 89268533, in cooperation with research investigator John E. Crawford.
Algorithm 1. Forward increment Lagrange multiplier method using multi-step central difference integration

Step 1 Given $\mathbf{u}_{n-1}$ and $\mathbf{u}_n$, predict velocities and displacements

$$ \dot{\mathbf{u}}_n = \frac{1}{h} \{\mathbf{u}_n - \mathbf{u}_{n-1}\} $$

$$ \mathbf{u}_{n+1}^* = h^2 \mathbf{M}^{-1} \{\mathbf{R}_n - \mathbf{F}(\mathbf{u}_n, \dot{\mathbf{u}}_n)\} + 2\mathbf{u}_n - \mathbf{u}_{n-1} $$

Step 2 Calculate contact forces and displacement corrections

$$ \lambda_n = [h^2 \mathbf{G}_{n+1} \mathbf{M}^{-1} \mathbf{G}_{n+1}^T]^{-1} \mathbf{G}_{n+1} \{\mathbf{u}_{n+1}^* + \mathbf{X}\} $$

$$ \mathbf{u}_{n+1}^c = -h^2 \mathbf{M}^{-1} \mathbf{G}_{n+1}^T \lambda_n $$

$$ \mathbf{u}_{n+1} = \mathbf{u}_{n+1}^* + \mathbf{u}_{n+1}^c $$

Step 3 Calculate velocities and accelerations

$$ \dot{\mathbf{u}}_n = \frac{1}{2h} \{\mathbf{u}_{n+1} - \mathbf{u}_{n-1}\} $$

$$ \ddot{\mathbf{u}}_n = \frac{1}{h^2} \{\mathbf{u}_{n+1} - 2\mathbf{u}_n + \mathbf{u}_{n-1}\} $$

Step 4 Next increment: $n \leftarrow n + 1$, return to step 1.

Algorithm 2. Gauss-Seidel iteration algorithm for two dimensional surface contact. This algorithm corresponds with step 2 of Algorithm 1

Step 2.1 Initialization

2.1.1 evaluate $\alpha_{n+1}^i$, by equations (17), (19) for $I = 1, 2, \dots, N$

2.1.2 initialize $^0\lambda_n^i = ^0\lambda_n^i = 0$, for $I = 1, 2, \dots, N$

2.1.3 initialize $^0, ^1\mathbf{u}_{n+1}^c = 0$

Step 2.2 Gauss-Seidel iteration $i$ to $i+1$, do steps 2.2.1 to 2.2.5 for $I = 1, 2, \dots, N$

2.2.1 $\{x_{AI}, y_{AI}, x_{BI}, y_{BI}, x_{CI}, y_{CI}\}^T = \{^{i,i+1}\mathbf{u}_{n+1}^c + \mathbf{u}_{n+1}^* + \mathbf{X}_I\}$

2.2.2 $\begin{vmatrix} ^{i,i+1}p_{xI} = (1 - \alpha_{n+1}^i)x_{AI} + \alpha_{n+1}^i x_{BI} - x_{CI} \\ ^{i,i+1}p_{yI} = (1 - \alpha_{n+1}^i)y_{AI} + \alpha_{n+1}^i y_{BI} - y_{CI} \end{vmatrix}$

2.2.3 $\begin{vmatrix} d = h^2((1 - \alpha_{n+1}^i)^2/m_{AI} + \alpha_{n+1}^i^2/m_{BI} + 1/m_{CI}) \\ \Delta\lambda_{xI} = ^{i,i+1}p_{xI}/d \\ \Delta\lambda_{yI} = ^{i,i+1}p_{yI}/d \end{vmatrix}$

2.2.4 $\begin{vmatrix} ^{i+1}\lambda_{xI}^n = ^i\lambda_{xI}^n + \Delta\lambda_{xI}^n \\ ^{i+1}\lambda_{yI}^n = ^i\lambda_{yI}^n + \Delta\lambda_{yI}^n \end{vmatrix}$
$$2.2.5 \quad i, i+1 \in \mathbf{u}_{n+1}^c \leftarrow i, i+1 \in \mathbf{u}_{n+1}^c - h^2 \left\{ \begin{array}{l} \Delta\lambda_{xI}(1 - \alpha_{n+1})/m_{AI} \\ \Delta\lambda_{yI}(1 - \alpha_{n+1})/m_{AI} \\ \Delta\lambda_{xI}\alpha_{n+1}/m_{BI} \\ \Delta\lambda_{yI}\alpha_{n+1}/m_{BI} \\ - \Delta\lambda_{xI}/m_{CI} \\ - \Delta\lambda_{yI}/m_{CI} \end{array} \right\}$$

**Step 2.3** Convergence test

$$\text{if: } \sqrt{\sum_{I=1}^{N} (\Delta\lambda_{xI}^2 + \Delta\lambda_{yI}^2)} < \text{tol} \sqrt{\sum_{I=1}^{N} (i+1)\lambda_{n}^2 + (i+1)\lambda_{yI}^2}$$

then: $\mathbf{u}_{n+1} = \mathbf{u}_{n+1}^* + i, i+1 \in \mathbf{u}_{n+1}^c$,

go to step 3 of Algorithm 1.

else $i \leftarrow i + 1$,

go to step 2.2 above.

**Algorithm 3.** Surface contact force conditions. Following step 2.2.3 of Algorithm 2, do steps 2.2.3a to 2.2.3d below, then proceed to step 2.2.4 of Algorithm 2

$$2.2.3a \quad \left| \frac{\sin = (y_{BI} - y_{AI})}{\sqrt{(x_{BI} - x_{AI})^2 + (y_{BI} - y_{AI})^2}} \right| \quad \left| \frac{\cos = (x_{BI} - x_{AI})}{\sqrt{(x_{BI} - x_{AI})^2 + (y_{BI} - y_{AI})^2}} \right|$$

$$2.2.3b \quad \left| \begin{array}{l} i+1 \lambda_{nI} = (i \lambda_{xI} + \Delta\lambda_{xI})\cos + (i \lambda_{yI} + \Delta\lambda_{yI})\sin \\ i+1 \lambda_{nI} = -(i \lambda_{xI} + \Delta\lambda_{xI})\sin + (i \lambda_{yI} + \Delta\lambda_{yI})\cos \end{array} \right|$$

$$2.2.3c \quad \left| \begin{array}{l} \text{if: } i+1 \lambda_{nI} > 0, \text{ then: } i+1 \lambda_{nI} = 0 \\ \text{if: } |i+1 \lambda_{nI}| > \mu|i+1 \lambda_{nI}|, \\ \text{then: } i+1 \lambda_{nI} = \mu|i+1 \lambda_{nI}|\text{sign}(i+1 \lambda_{nI}) \end{array} \right|$$

$$2.2.3d \quad \left| \begin{array}{l} \Delta\lambda_{xI} = i+1 \lambda_{nI} \cos -i+1 \lambda_{nI} \sin -i \lambda_{xI} \\ \Delta\lambda_{yI} = i+1 \lambda_{nI} \sin +i+1 \lambda_{nI} \cos -i \lambda_{yI} \end{array} \right|$$

## REFERENCES

1. K. J. Bathe and A. B. Chaudhary, 'A solution method for planar and axisymmetric contact problems', *Int. j. numer. methods eng.*, **21**, 65-88 (1985).

2. T. B. Belytschko and J. I. Lin, 'A three-dimensional impact-penetration algorithm with erosion', *Comp. Struct.*, **25**, 95-104 (1987).

3. J. Bohm, 'A comparison of different contact algorithms with applications', *Comp. Struct.*, **26**, 207-221 (1987).

4. S. K. Chan and I. S. Tuba, 'A finite element method for contact problems of solid bodies—Part I. Theory and validation', *Int. J. Mech. Sci.*, **13**, 615-625 (1971).

5. A. B. Chaudhary and K. J. Bathe, 'A solution method for static and dynamic analysis of three-dimensional contact problems with friction', *Comp. Struct.*, **24**, 855-873 (1986).

6. J. O. Hallquist, G. L. Goudreau and D. J. Benson, 'Sliding interfaces with contact-impact in large-scale Lagrangian computations', *Comp. Methods Appl. Mech. Eng.*, **33**, 107-137 (1985).

7. T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier and W. Kanoknukulchai, 'A finite element method for a class of contact-impact problems', *Comp. Methods Appl. Mech. Eng.*, **8**, 249-276 (1976).
8. J. W. Ju and R. L. Taylor, 'A perturbed Lagrangian formulation for finite element solution of nonlinear frictional contact problems', *J. Theor. Appl. Mech.*, Special Issue, Supplement No. 1 to 7, 1-14 (1988).

9. M. G. Katona and O. C. Zienkiewicz, 'A unified set of single step algorithms, Part 3: The Beta-m method, a generalization of the Newmark scheme', *Int. j. numer. methods eng.*, **21**, 1345-1359 (1985).

10. J. A. Landers and R. L. Taylor, 'An augmented Lagrangian formulation for the finite element solution of contact problems', *SESM Report No. 85/09*, Department of Civil Engineering, University of California, Berkeley, 1985.

11. J. J. Moreau, 'Unilateral contact and dry friction in finite freedom dynamics', in J. J. Moreau and P. D. Panagiotopoulos (eds.), *Nonsmooth Mechanics and Applications, CISM Courses and Lectures No. 302*, Springer-Verlag, Berlin, 1988.

12. N. M. Newmark, 'A method of computation for structural dynamics', *J. Eng. Mech. Div.*, ASCE, **85**, 67-94 (1959).

13. J. W. Simons and P. G. Bergan, 'A finite element formulation of three dimensional contact problems with slip and friction', *Comp. Mech.*, **1**, 21-35 (1986).

14. R. L. Taylor, N. J. Carpenter and M. G. Katona, 'Analysis of non-linear transient problems with constraints', *Int. j. numer. methods eng.*, to appear.