Computer methods

In applkd
moohanks and
snglnrerlng
EISEVIER Comput. Methods Appl. Mech. Engrg. 143 (1997) 229-247

A method for contact problems using virtual elements

W.R.C. Underhill”‘“, M.A. Dokainish”, G.AE. Oravasb
“Department of Mechanical Engineering, McMaster University, Hamilton, Ontario, Canada
bDepartment of Civil Engineering and Engineering Mechanics, McMaster University, Hamilton, Ontario, Canada

Received 14 October 1992

Abstract

A direct optimization method for contact problems is shown that uses separate corrective force vectors for equilibrium and
compatability. The geometric mediation across the contact region is accomplished by defining virtual nodes and elements.
Integrations required on each contact surface are carried on a virtual mesh that reflects the mesh of the body for which forces are
required. Examples demonstrate that the method is accurate and robust.

1. Introduction

Contact problems are so common that it is surprising that they are so little understood. Examples
include stones resting on another, every bolt, screw, nail or rivet holding parts together, the bearing of
skeletal joints, the cutting or grinding of tools on workpieces, and many more than can be listed. There
are many times when there are advantages to a calculational approach to finding information needed
about contact problems. These include situations where a physical experiment must be avoided. This
could happen for reasons of expense, or of safety, or the facilities to do the experiment may not be
available. Nonetheless, detailed knowledge of conditions at a contact interface may be desired. So,
there is plenty of motivation to have a way to simulate contact.
Results that are exact to the extent of the mathematical theories of deformable bodies are almost
unknown for contact problems. Approximate results are what we can currently provide. Even these can
be found in more or less accuracy and detail. Here we consider only the contact of two or more solid
bodies.
The work of Heinrich Hertz in the nineteenth century remains one of the fundamental reference
points in the study of contact problems [1,2]. Hertz provided some analytical solutions to the problem
of linearized contact over very small regions, for elastic materials, without friction. In particular, he
approximated any surface by an ellipsoid with the same principal curvatures at the point of contact. The
case of contact at a concavity is represented by one or both radii of curvature being negative. For
techniques to solve more complicated problems this is usually considered a simplified version and so
later works attempt to demonstrate their validity and accuracy by presenting solutions to similar
problems.
Hertz showed that the contact region was an ellipse lying in a plane and that the pressure distribution
was an ellipsoid. The total pressure integrated over the contact ellipse equilibrates the forces applied to

* Corresponding author.

00457825/97/$17.00 0 1997 Published by Elsevier Science S.A. All rights reserved

SSDI 0045-7825(94)00764-O
the distant parts of the body that cause such motion as to result in contact. The distribution of stresses
can be described in closed form. The highest compressive stress, which appear at the center of the
contact, is 3/2 of the average. The maximum tensile stress appears at the edge of the contact ellipse in
the outward direction. The maximum shear appears on the normal axis through the center of contact.
The details of magnitudes and distribution of pressure and stress vary with the relative magnitudes of
the radii of curvature and the angle between the normal planes containing the principal curvatures.
Many tables have been constructed to provide calculational assistance to the analyst who would use this
approach.
For cases where these assumptions are not valid, we use the finite element method. Unfortunately,
there is a change from the usual situation in the nature of the exterior information. This leads to strong
non-linearities in the solution. We no longer necessarily have foreknowledge of either the stresses or the
displacements. The style of information on which we have relied to establish our substitute algebraic
problem is simply not available. Instead of having one or the other of the two kinds of information that
we use outside the areas of contact, we have a combination of two pieces of information that are
appropriate inside the region of contact: the Principle of Impenetrability and the Principle of Equal and
Opposite Reaction. These two principles contain all of the information needed to solve the contact
problem. It remains to find a way to express them mathematically and to put that to use.

1.1. Numerical approaches

Attempts to express these principles and use them can be lumped into four main branches. These are
those based on the use of Lagrange multipliers, those based on the penalty methods, those based on
augmented Lagrange multipliers and the direct methods. The minimization problem can be stated as:
Find the value of the nodal displacement, U. such that the potential energy is minimized while not
violating the Principles of Impenetrability of Equal and Opposite Reaction.
The first part is easily expressed as

Find Min n(u) , n+“.‘K”-u’F (1.1)

The second part is more difficult. How does one express the Principle of Impenetrability? There are two
common approaches to this. In one the separation of the bodies is measured at some time and it is
required that the displacements not exceed this gap function. The gap function may be established at
the outset or re-evaluated at every iteration. The other common approach is to detect penetration by
inspection and then do something about it.

1.2. Lagrange multipliers

These typically employ a gap function, g. To show the flavour of this approach consider first the case
of frictionless contact between a linearly elastic deformable body and a fixed, rigid barrier. As any point
in the body approaches the barrier, the value of g for that point decreases and reaches zero as the body
point just touches the barrier. So, it is required that

g(X) 2 0 (1.2)

where g is defined only for X lying on the surface of the deformable body. To enforce this condition
multiply g by a free variable, A, the Lagrange multiplier. Add this to the original potential. The
modified potential, 17*. is

(1.3)

It is also convenient to express g as a Taylor series in U, the displacements, up to linear terms as

g(U)=gl,+$U or g(U)=g,,+GU (1.4)

 W.R.C. Underhill et al. I Comput. Methods Appl. Mech. Engrg. 143 (1997) 229-247 231

where G is the derivative of the gap function with respect to the nodal displacements. Now
minimization of II* with respect to both the nodal displacements, U, and the Lagrange multipliers, A,
gives

(1.5)
Clearly, there are some disadvantages to this style of imposing constraints. Now for every node to be
constrained there is an extra variable to find in passing through the linear equation solver. There is also
a block of zeroes on the main diagonal which will disturb many available programs.
What is to be gained by solving for all these extra variables ? As is often the case in Lagrange
multiplier techniques, the multiplier is itself a quantity that can be interpreted physically. In this case it
is the contact pressure. This gives us further insight into how to reform the above for the case of two
deformable bodies in contact. We can define a contact potential, I&, which represents the work done by
the contact pressures.

n, = gT(UA,U”)A (1.6)
where superscript A and IEBlabel the bodies. The potential to be minimized is the sum of the mechanical
and contact potential. This gives

(1.7)

There are many variations on this. As presented so far, this is only suitable for analyses involving
small displacements and linear elastic materials. If the displacements become large, or if there are
material non-linearities or if there is friction involved it becomes necessary to change to an incremental
approach. The basic idea is much the same. Among the Lagrange multiplier methods there can be a
lumping into those which concentrate on the gap function and derive values for G by geometric
arguments and those which concentrate on the contact potential and derive G as an integral over the
contact area of derivatives of 17c with respect to the nodal displacements.
Among those that use the geometric approach the discussions range mostly on the topic of the
difficulties associated with sliding motion with friction. Ordinarily, this leads to extra terms being added
to the contact forces to represent the effects of friction. The result is that the system matrix becomes
asymmetric. This is unfortunate since this doubles the computer storage requirements. Pascoe and
Mottershead [3] address the problem of asymmetry of storage. The advantage in reducing the
requirement for computer storage would be important for many users. It is proposed that, for small
incremental steps, it is acceptable to replace the displacement constraints with the transpose of the force
constraints. Since the force constraints include extra terms for friction, the contact forces are not
normal to the contact surface. Likewise, this attempts to make the nodes slide along a surface that is
not tangent to the physical surface. The tangent of the angle between the two surfaces is the coefficient
of friction. It is pointed out that for surfaces with little curvature this is easily corrected in later
iterations and for surfaces with great curvature the steps have to be so small that the errors of this
method are not a problem.
The group of techniques that integrate the contact potential explicitly have much the same
experience. The integrals are usually taken over the faces of the elements in contact rather than over
some independent discretization but may use different shape functions. These are often referred to as
contact elements and often called mixed because of the explicit representation of stresses. The results
tend to oscillate for Hertzian problems.

1.3. Penalty methods

The penalty methods take a different approach. Essentially, these are attempts to change the
problem from kinematically non-linear to materially non-linear. Instead of trying to build the equations
232 W.R.C. Underhill et al. ! Compuf. Methods Appl. Mech. Engrg. 143 (1997) 229-247

so that only admissible solutions will appear, they attempt to detect and correct errors as they occur.
This is done by associating a large penalty of energy with committing the error. The penalty is set up in
such a way that nodal motions toward an admissible solution decrease the penalty to be paid and so
minimizing the energy of the problem tends to lead to acceptable results. For contact there are two
main approaches to applying these penalties: the filled gap and the penetration approaches. They differ
in where they allow compatibility errors to be committed. The filled gap techniques fill the space
between the bodies with a fictitious material. After displacement a thin sheet of this material remains
between the bodies. The penetration detection methods apply corrective forces proportional to the
depth of penetration. After displacement only small penetrations are left because large forces prevent
greater violation of the compatibility constraints.
In the simplest version of the filled gap method, Osmont [4], watch is kept for pairs of nodes about to
touch. When the nodes are very close together a sort stiff spring is inserted between them. In other
versions the gap material can be quite sophisticated. Zolti [5] points out that while the mechanical
properties are fictitious, the chemical, thermal and other properties may represent real conditions in the
gap. The constitutive properties in the direction across the gap must be highly non-linear. As the
normal strain across the contact elements gets to be very close to -1 (the gap has nearly closed) the
modulus of the fictitious material rapidly becomes very large. Constitutive properties in the directions
tangent to contact can be arranged to imposed sticking or sliding friction, or no friction at all. This
method of coupling is only suitable for small relative displacement unless some scheme is included for
rediscretizing the gap material into new elements that reflect the current arrangement of the nodes and
elements of the real bodies.
The penetration methods add to the mechanical potential a penetration potential, JZ,, proportional to
the square of the penetration depth, P. The penetration depth is treated similarly to the gap function
used in the Lagrange multipliers. It is useful to express it as a function of the nodal displacements. So

(1.8)
So the system of equations to be solved is
K”+K;
BA (1.9)
K,
The penalty terms, K, and Fp, are activated for any node in contact and are zero for the rest. Because
the penalty terms are large, the system matrix has similar mathematical difficulties as problems that
involve non-linear materials with rapidly varying constitutive relations. Such analyses also occasionally
have problems of ill-conditioning. For the simpler schemes of Cheng et al. [6], the penetration function
is found by a node matching method much like that used in the simple filled gap methods. Yagawa and
Hirayama [7] fill the space between the bodies with eight-noded brick elements. These are not filled
with contact material, but are used to establish the value and derivatives of p.

1.4. Augmented Lagrange multiplier methods

There are other approaches that are combinations of the Lagrange multipliers and the penalty
method. These hybrids are called the augmented Lagrangian methods. So, in one version the potentials
included are the mechanical potential, n, the penetration potential, II,, and the contact potential, II,.
The modified potential is
n*=n+n,+nc (1.10)
The resultant system of equations is
KA+ KtA K;” GA’ U’
BA
K, KB+$FB GfT][ yB] = [ :<!I (1.11)
GA
Another mixture, also called the augmented Lagrangian method, attempts to improve on the
 W.R.C. Underhill et al. I Comput. Methods Appl. Mech. Engrg. 143 (1997) 229-247 233

shortcomings of the method of Lagrange multipliers, so that the multipliers have their own inequality
constraint. A penalty method is used on the Lagrange multipliers, themselves. This helps overcome
numerical imprecision in the Lagrange multiplier method and helps damp the oscillatory behaviour that
often occurs when these methods are used for contact. These include a separation potential, Us, and
minimize as before. In assembled form this is

(1.12)

The augmented Lagrange penalty matrix, CX~,is diagonal. Each entry is zero for those multipliers that
are negative and a large number for those that are positive. This is the form most often used in the
literature on contact.
One may imagine combining both forms of augmented Lagrangian methods. This might improve the
performance of the penalty method by restricting the solutions and tend to keep the Lagrange
multipliers negative. No examples of such an approach were found in the literature for contact
problems.

1.5. Direct optimization methods

These attempt to solve the contact problem by finding appropriate increments to the forces acting on
the contact surfaces. It is assumed that the modelling of forces has been inadequate because of errors in
the terms for forces due to contact. Correction of these forces is to give correction to displacements and
reach the minimum potential energy.
A simple version is found in the pioneering work of Francavilla and Zienkiewicz [8]. An error in the
displacements is calculated as the depth of penetration. The error in forces due to missing contact terms
is found as the force required to cause a displacement increment to correct the displacement error.
AF=KAU (1.13)

In fact, many of these methods have set the equations up in terms of the flexibility matrix, K-‘. This
approach is extended to sliding with Coulomb friction by Sachdeva and Ramakrishnan [9]. This is
extended to a more general approach in the work by Torstenfelt [lo].
The inclusion of friction and the manner of choosing AF lead to the differences in the various direct
optimization methods. The majority of methods explicitly identify node pairs. This has the advantage of
reducing the number of equations to solve and a built-in guarantee of nodal equilibrium between
bodies. However, such methods are inherently restricted to problems where relative displacement of
the nodes is small. If two nodes identified as a pair slide apart by a distance comparable to the size of an
element, then the problem must be halted and rediscretized. It also has the disadvantage that the bodies
in contact must be kept in storage simultaneously.

1.6. Friction laws

The effects of friction are of great importance in some contact calculations. The simplest friction law
is none at all. This is also a common condition to use in the development of analyses that do use
friction. It can be used to check the qualitative behaviour of such an analysis. It is also used for
quantitative checks, using the closed form solution for Hertz contact.
The next simplest friction law is sticking contact. This can be easily implemented since under these
conditions the nodal displacements of one body can be eliminated in terms of the nodal displacements
of the other.
Sliding friction is a source of a great deal of effort in numerical simulation of contact. It raises a
number of issues. Now there are two different friction regimes to consider. It is necessary to explicitly
state what kind of friction law is to be used. Almost without exception, the law proposed by de
Coulomb is used in the literature. This was originally meant as a phenomenological description of the
234 W.K.C‘. Underhill et ul. i C’omput. Methods Appl. Mech. Engrg. 143 (1497) 229-247

averaged properties of the sliding of blocks of material with macroscopically flat surfaces; not as a
detailed, microscopically correct model. Nonetheless, Coulomb friction has been used in many
calculations as if it were an accurate description of the microscopic conditions at a point. Or, to use the
words of the study of continuum mechanics, Coulomb friction is used as if it were true for every
differentially small neighbourhood. It leads to the use of various terms at the contact surface that are
not differentiable or continuous. For terms involving friction the mathematical inconvenience is similar
to problems that involve phase change boundaries. The positions of the boundaries between sticking
and sliding friction are not necessarily at the boundaries of the elements. This has led to efforts by some
workers to arrange that over any particular element only one friction regime is present. The most
successful and elegant are those done with the moving finite element methods of Haber [ll], of Haber
and Harandria [12] and of Gu [ 131.
Proponents of ‘non-local’ friction laws have several approaches. The simplest and commonest way is
to concern themselves only with the resultant forces at the nodes. These forces are necessarily some
kind of average of the frictional conditions between the nodes. Another way is used by Shyu et al. [14],
Chang et al. [15] and Simo et al. [16]. These workers have used an augmented Lagrangian method
based on the Lagrange multipliers to restrict the set of possible solutions to the set of admissible
solutions. Instead of using nodal Lagrange multipliers they use elemental ones. This gives the averaged
contact pressure over the entire element and is only C ’ continuous. It has interesting side effects on
convergence and accuracy. Contact equilibrium and compatibility are satisfied exactly only in the limit
of fine meshes. In practice it may not be necessary for the mesh to be extremely fine. This allows very
easy integration of the stresses.
One of the most interesting versions of the Coulomb friction law arises from the observation that in
many ways this law behaves very much like the constitutive relation of an elasto-plastic material. The
work done by slippage against friction cannot be recovered; the slippage only occurs if a certain
criterion of stress is met; the frictional stresses never exceed a particular bound. Padovan et al. [17],
Zhong and Sun [18] and Plesha et al. [19] use such an ‘elasto-plastic constitutive relation’ for friction
and slippage.

2. Virtual elements

A flow chart for a typical non-linear algorithm (materially or geometrically non-linear) is shown in
Fig. 1. It is assumed that the solution is known up to a certain point. A small increment of the known
loads is applied and the corresponding increment in the unknown displacements is found. Since the
internal information (material properties) may be dependent on the new conditions it is necessary to
check whether the internal information has changed. It it is found that it is necessary to update, then
that is done. This implies that the external information (loads and displacements) may no longer be
correct, or least not sufficiently close. A simple and popular iterative method is to estimate the load
increment from the stiffness and the displacement increment. The difference between the estimated
loads and the actual increment of the loads is substituted and solved again. After convergence is
achieved the algorithm is started again with a new increment of external information. For mechanically
stable structures this process will always converge.
The key point is the test for equilibrium. This is done after the imposition of current conditions on all
known quantities; usually stiffness, displacement and forces.
Fig. 2 shows a prototype flow chart for a contact algorithm. The key difference between Figs. 2 and 1
is in the nature of the test which triggers another iteration. In the contact algorithm the test is for two
possible triggers: compatibility and equilibrium. Different algorithms distinguish themselves in how they
detect and generate corrections.
The method developed here is a direct optimization method. The only nodal degrees of freedom are
the displacements. The stiffness matrices are those that represent actual material characteristics. There
are no extra terms in the system matrices to represent attempts to arithmetically coerce the behaviour
of the solution. So, this can be added to any displacement-based finite element program if the coding is
adequately modular.
 W.R.C. Underhill et al. I Comput. Methods Appl. Mech. Engrg. 143 (1997) 229-247 235

Fh F+AF

Use Iterative force Calculate force

increment for next

Fig. 1. Flow chart for a non-linear problem.

The treatment that each body receives is exactly the same as every other. There is no nomination of
one body as target or contractor or master or slave. There is no explicit coupling between the
displacement degrees of freedom of the bodies. This allows the information for each body to be stored
separately, instead of all in one huge matrix. This alleviates the restrictions inherent in some node
pairing schemes and allows any amount of sliding of the touching surfaces without requiring
rediscretization.
The algorithm is objective. All corrections are applied as force terms rather than displacement
increments. These do not have to be applied before other non-linearities but are applied at the global
level. So, it does not have to be nested inside other non-linear coding.
Each kind of error, compatibility and equilibrium, is treated separately. To the first order the
corrective force increment for compatibility errors does not affect equilibrium. Likewise, to the first
order the corrective force increment for equilibrium errors does not affect compatibility. This is
different from other techniques that require a single force increment to correct both kinds of errors.

2.1. One-dimensional analogy

In a one-dimensional setting the measure of a penetration error is the length of overlap. Forces are to
be applied to both bodies to correct this situation. So, the compatibility condition on corrective
displacements is
-AP + AU” = P (2.1)
where P is the penetration error and AU is a corrective displacement. If the displacements are small,
then to the first order the change in the end reactions must be
236 W.R.C. Underhill et al. I C‘omput. Methods Appl. Mech. Engrg. 143 (1997) 229-247

Solve [K,“,] AU = A+---

update Kc-- K(U+AU) Stiffness

FM F+AF Loads
U+- U+AU Displacements

Check equilibrium
AND
Check compatability

NO Yes
Are both satisfactory?

1 1

Find some kind of Calculate force

increment for next

Fig. 2. Flow chart for a contact problem.

AF” = KA AU’
(2.2)
AF” = K” AU’
where AF is a change in end reaction and K is a stiffness. These displacements are to correct
compatibility without disturbing equilibrium. So, the condition is put on the reactions that

AF”+AF”=O (2.3)
These conditions can be combined to give the required increments of force to correct a compatibility
error without disturbing equilibrium.

AUA = -K”Pl(K’ + K”)

AU’ = KAPI(KA + KB)

(2.4)
AFA = - K”K”PI(KA i- KB)

AF@ = KAK”PI(KA + K”)

An error in equilibrium is simply defined as well. It is just the difference between the net force on the
end nodes of each body. The force increments to be applied to correct an equilibrium error are required
not to upset compatibility. This requires that

AU”= AU” (2.5)

 W.R.C. Underhill et al. I Comput. Methods Appl. Mech. Engrg. 143 (1997) 229-247 231

The effect of the condition in Eq. (2.5) is that the end reactions are tied as well through Eq. (2.2).
These force increments are to correct an equilibrium error. So we require of them that
AFW + AF” = -E, (2.6)

where E, is the equilibrium violation. These conditions may be combined to give

AF’ = -KWE,I(KW + K”)

(2.7)
AF” = -K”EFI(KW + K”)

These are applied at the global level along with whatever other force correction terms may be used for
the sake of other non-linearities.

2.2. Complications of two- or three-dimensional contact

There are complications that arise in contact in two or three dimensions. They include:
(a) There is more than one direction for compatibility and equilibrium considerations.
(b) There may be more than one node of each body in each contact area.
(c) Stiffness is now a matrix quantity rather than a scalar.
(d) The nodes in contact may not be touching nodes in the other body.

(a) More than one direction

In two-dimensional problems one has to consider equilibrium in directions normal and tangent to the
contact surface. In three-dimensional contact there are two tangent directions. The chief difficulty is in
deciding which is the normal direction to the surface. This is affected by the assumptions which are
made on where the true contact surface lies. In this work the difficulty is handled by the assumption
that the normal to the exterior of the elements in the current configuration is an adequate
approximation to the normal to the contact surface. As the elements are deformed during iterations for
contact or any other cause, the estimate of the normal direction is updated. As the iterations converge,
the element normals improve as estimates of the normal to the contact surface.
It is also necessary to use a depth of penetration. This is simple in one dimension but must be defined
in some way for spaces that have more than one dimension. This choice of definition alone could be the
subject of much argument. Here, the choice is made that the depth of penetration is measured in the
same direction as the normal forces in the equilibrium check.

(b) More than one node in contact

In the one-dimensional case there was no doubt as to which nodes were involved in contact. Only the
nodes at the end of a body were ever considered. In two or three dimensions the exterior of a body is
no longer just one node at either end. Each node on the exterior is a possible contact node. Of
necessity each one must be checked for contact. A node is considered to participate in contact if
(1) it is an exterior node AND
(2) is embedded in any other body OR
(3) retains compressive contact forces acting on it from previous calculation

(C) Stiffness is a matrix

For purposes of computation there are many approaches to determining which nodes are in contact
and which are not. An important feature of calculation that is implicit in the one-dimensional case
becomes explicit in the two- and three-dimensional cases. Any node that is not in contact (all interior
nodes and some exterior nodes) may be eliminated statically from the calculation of contact forces. This
greatly reduces the number of degrees of freedom involved and so reduces the computational difficulty.
The important conceptual step gained is that the entire body has been represented by only the contact
surfaces. All deformations due to contact will exhibit the correct global strain energy and work.
In the one-dimensional case the stiffness is a scalar. This is coincidental. One may correctly claim that
the stiffness is a 1 x 1 matrix, rather than something that is inherently a scalar. In two- and
 23x W.R.C. Underhill er al. I Compur. Methods Appl. Mech. Engrg. 143 (1997) 229-247

three-dimensional problems the stiffness that remains after eliminating the non-contact nodes is a
matrix. The number of degrees of freedom is the number of degrees of freedom per node times the
number of participating nodes. The stiffness is treated in two stages:
(1) Degrees of-freedom not associated with contact nodes are statically eliminated.
(2) The remaining stiffness matrix is reduced statically to each contact degree of freedom. In this
stage of reduction all remaining degrees of freedom are considered as free variables. Any
displacement prescriptions are ignored. The remaining stiffness is a measure of the strain energy
required throughout the body for a motion to occur in that dof, as such this resembles a
Castigliano stiffness. The actual calculation of the individual reduced stiffnesses is done
recursively to minimize the computational effort.

(d) Nodes do not align

In realistic problems the meshes of opposing bodies will not necessarily match over the region of
contact. So each node does not always find an opposing node, the forces and stiffnesses of which can be
used to estimate the corrective forces for compatibility and equilibrium. Different quantities are
mapped across the contact surface in different ways. Stiffness is among those found by interpolation
between opposing nodal values. Forces are among the quantities that are redistributed and reintegrated
over the mesh of node’s own body.

2.3. Redistribution of nodal quantities

Some nodal quantities are found as the integrals of distributed quantities. For example, nodal forces
can be found as the integral of stress over the surface of the element. In more detail, there are several
steps in calculating a nodal force from surface traction. First, the traction is weighted by the shape
functions. Then the product is integrated over the element to give elemental contributions to the nodal
forces. Finally, in the assembly of elements, the individual elemental contributions are accumulated to
give the nodal force.
This process may be illustrated here for an unspecified distributed quantity q. q is shown as a scalar.
The argument is more long winded if q is of higher vector order, but is essentially the same. It is
necessary in that case to expand the calculation to treat each component of q separately. First form q as

9 = 44 (2.8)
where 4 is a column of shape functions, q is the unspecified scalar field and q is a column of the same
length as 4. We define a global equivalent quantity, Q, concentrated at the nodes. The contribution
from an element is calculated as

Qel=j- O’=l,l +qdV (2.9)

%I rl
where Q,, is the elemental contribution to the concentrated nodal quantity and O,_, is the elemental
domain.
Finally, the global concentrated nodal quantity, Q, corresponding to the distributed quantity, q, is
found by assembly. That is

Q=CQ,, (2.10)

where C is the assembly operator.

There is a certain amount of information lost in this process. The exact details of the distribution of q
are not important to the result, Q. Any other distribution, say q’, which has the same integral after
weighting by the shape functions will produce the same result. In particular, one may imagine replacing
q by some distribution, q’, of form
qZq’=+Tq’ (2.11)
where q’ is an approximation for q and q’ is a column of pseudo-nodal values. Let Q’ be the
 W.R.C. Underhill et al. I Comput. Methods Appl. Mech. Engrg. 143 (1997) 229-247 239

concentrated nodal value found by using q’ instead of q. Q’ is then an approximation for Q. If the
values of q’ are chosen correctly, then it may be possible to match Q’ to Q exactly.
For the purposes of Eq. (2.9) it would be allowable for each element to maintain its own vector, q:,,
of pseudo nodal values. Many finite element codes use this form for the input for distributed quantities
such as pressures and thicknesses. But, to ease our calculation, let us use only a global list, with the
result that every element which involves any given node uses the same value for q’ at that node. This
also gives that the approximation, q’, is continuous even though there was no requirement that the
original function, q, is continuous. We can now express Q6, as

(2.12)

where qk, are those entries of q’ for this element. Since the q’ are constants, rather than functions of
position they may be taken outside the integral

(2.13)

We may define the elemental distribution matrix, N,,, for each element as the integral in Eq. (2.13).

Ne,= I R,, MT dV (2.14)

So Eq. (2.13) may be written as

Q:, = Ned, (2.15)

The assembly operation may be applied to the Qd, to give Q’

Q’ = x Q:, = c Wed,) (2.16)

The assembly of the qil is just the global q’.This leads to a definition of a global distribution matrix, N,
as the assembly of elemental distribution matrices, N,,. So

N=h’,, (2.17)

This allows the global approximation to the concentrated nodal quantity to be expressed concisely as
Q’=Nq’ (2.18)

The approximation q’ can be said to adequately represent the original function q if both give the same
global concentrated nodal quantity. Since the only degrees of freedom in q’ are the q’, all adequacy
conditions rest on them. The condition put above can be fulfilled by choosing q for the q’, where q is
the solution of
W=Q (2.19)

This requires that the global distribution matrix, N, not be singular. For many situations this will be
true. It may be noted that N,, has the same form as the consistent mass matrix for many classes of
elements. For any such class that has a non-singular consistent mass matrix the global distribution
matrix is also non-singular.
This manner of approximating q not only loses the local detail of the distribution, but also may lose
larger features. For example, suppose that some region (element or set of elements) had q of exactly
zero over it. q’ would not necessarily also show this but would tend to average the value among those
elements and their neighbours. In the interior of such a region the value of q’ may be quite close to zero
or even exactly zero but near the periphery the results change smoothly to the surrounding non-zero
conditions. Care can be taken to allow for different results. In the above example it would be possible
to suppress q’ over an element by imposing values of zero for q’ for all nodes in the affected elements
before solving Eq. (2.19). This would have the effect of forcing the distribution to zero at the edges of
240 W.R.C. lJnderhiil et al. I C‘omput. Methods Appi. Mech. Engrg. 143 (fY97) 229-247

adjacent elements as well. An alternative possibility would be to introduce a deliberate discontinuity.

This is done by omitting the affected elements from the assembly stage to give a modified global
distribution matrix N*. Substitute N* into Eq. (2.19) for N. The pseudo nodal values 9’ so obtained are
not restricted to be zero around the edges of the unaffected elements. It must be remembered in such
an approximation that a separate approximation of zero has been made for the omitted region.
In the example above any function that is known a priori can be used. The part of the distribution
that is already known need not be zero. The point to be made is that in some circumstances care may be
needed to achieve desired results or else a greater than minimal magnitude of approximation may
occur.

2.4. The virtual mesh

The equations for the corrective force increments described in the one-dimensional explanation
assume that each node in contact is touching a node of the opposing body. This is achieved in a trivial
fashion for one-dimensional problems. This was achieved in early two-dimensional finite element codes
by restricting the allowed mesh of nodes and elements so that each node in one body was paired with a
node in the touching body. Only very small misalignments of these nodes could be tolerated.
In the present scheme this difficulty is removed by the definition of a mesh of virtual nodes and
virtual elements that exactly oppose the actual mesh of either body. Thus each contact node, or
element, of body A will find a virtual node, or element, of body B touching it. Likewise the mesh of
body B will find the virtual mesh of body A touching it. These virtual nodes provide the comparison for
the definition of p, the penetration error, E,, the equilibrium error, and K, the opposing nodal
stiffness. The use of the properties of these virtual nodes and elements gives the analyst freedom in the
design of meshes and lifts the restriction that actual nodes must stay close together in pairs. The
properties of a virtual mesh can be calculated at need.
The needed properties are found in a simple way. For the sake of discussion, assume that a node of
body A has penetrated into the bulk of body B. For a node embedded in another body the virtual node
is located on the surface of the opposing body at a point where an inward normal would point to the
embedded node. Penetration depth is measured as the separation of the real (embedded) node of A
and the virtual (surface) node of B as shown in Fig. 3.

2.5. Equilibrium corrections

For this we must find the stiffness and forces of each pair of a node and its opposing virtual node. The
forces and stiffness of the actual node are available directly.
To find the stiffness of the virtual B node, first find the coordinates of the virtual B node in the actual
B element. Then interpolate between the nodal stiffnesses of the actual B nodes. Since the nodal
stiffnesses are available as Castigliano style total stiffness, they represent the resistance of a body to
motion in the degrees of freedom. Interpolation serves well to find an estimate of the resistance to
motion at a location (virtual node of body B) between these known points (actual nodes in body B).
To find the forces acting on the virtual B node, first the actual forces on the actual B nodes are
redistributed as tractions. This gives an approximate description of the contact tractions acting on body
B. The virtual B element is chosen to occupy the same region and uses the same shape functions as the
actual A element, but it is on the surface of body B. The B tractions are integrated over the virtual B
element, weighted by the shape functions in the ordinary manner of finding nodal forces from tractions.
Each virtual elemental force vector is accumulated so that the complete virtual nodal forces are
available when all of the virtual elements have been assembled.
The forces and stiffnesses on both the node and the virtual node are re-expressed in a coordinate set
with the first coordinate in the normal direction. Equilibrium checks are made in the normal and
tangential directions. If a node is found to be in tension in the normal direction it is released. Otherwise
corrections are applied to distribute the force imbalance as in Eq. (2.7).
 W.R.C. Underhill et al. I Comput. Methods Appl. Mech. Engrg. 143 (1997) 229-247 241

BOOY B

PEICIRATION
DEPTH

Fig. 3. Location of a virtual node.

2.6. Compatibility corrections

This is accomplished in two stages. In the first stage, the position of the contact surface is estimated
and corrective displacements proposed. The penetration depth of a node is the distance between the
real node and the virtual node. In the case of a node that has penetrated, there is a region where
overlapping volumes must be resolved. In the case of a node that is outside of the opposing body, but
has not been released (still has compressive contact forces acting on it), this reveals that the current
estimate of contact forces is too high. These forces must be reduced to allow the node to approach the
contact surface. In either case, the proposed direction of approach to the contact surface is always along
the normal. The distance to approach such that equilibrium should be maintained, to first order, is
given in Eq. (2.4).
In the second stage, the corrective force increment is found from these desired displacements. The
forces given in Eq. (2.4) could be used, but these ignore interactions between nodes. It is more
effective to use the displacements and convert them to forces by multiplication by the stiffness matrix.
That way all the interactions between nodes are preserved in the force increments that correct
compatibility errors. Since the stiffness matrix is already factored and available from the linear
prediction stage, there is very little cost for this way of converting the compatibility correcting
displacements to forces.

3. Results

Four examples are shown to illustrate the success of the contact algorithm. The first example is a
Hertz problem. This is to demonstrate that the results are realistic and accurate. The calculated results
are compared to the theoretical solution. The second example is a punch problem. A shallow V shaped
wedge is pushed onto a flat plate. The elastic properties of the wedge and the plate are varied to
explore the effects of the relative stiffness of the bodies in contact. The third example is two rectangular
blocks of different sizes. This shows the behaviour in the presence of sharp corners. The results are
242 W.R.C. Underhill et ul. i C‘ompur. Method.5 Appl. Mech. Engrg. 143 (lYY7) 229-247

robust even in the presence of such challenging geometry. The last example is a weak gasket pushed
against a stiff gasket and then retracted. The displacement history shows that the method is reversible
where that would be expected for physical reasons.

3.1. Hertz example

This is a plane strain example. An infinitely long right circular cylinder of isotropic linearly elastic
material is pushed against a rigid flat frictionless surface. The cylinder has unit radius; the elastic
material has a modulus of 10h and a Poisson’s ratio of 0.2. The ‘rigid’ surface is modelled as a line of
elements, completely fixed on the lower surface, with an elastic modulus of 3 x 10’. That is, the
material is 3 x lo3 times as stiff as that of the cylinder. Together with the support conditions, the lower
body is effectively rigid. The mesh is shown in the original configuration in Fig. 4. Since the problem is
symmetric about the central plane, only half of it is modelled. At the end of motion a force of 1530
units has been applied to the half model. The final configuration of elements near the contact region is
shown in Fig. 5. The normal contact stress is plotted against distance from the central plane in Fig. 6.
The solid quarter circle shows the theoretical results. The horizontal lines show the range of the
calculated results. The vertical scale is normal contact stress divided by the maximum normal contact
stress. The uncertainties in the data on this scale are negligible. The horizontal axis is horizontal
position divided by the radius of contact. As the calculated results can only say that the radius of
contact is somewhere in element 207 a range is shown. The left end of each datum (shown with a hollow
square) uses the extreme outer edge of element 207 (node 208) as the upper limit of the radius of
contact. The right end of each datum (shown with a filled square) uses the extreme inner edge of
element 207 (node 207) as the lower limit of the radius of contact.
The data follow the theoretical results quite well. This shows that the method is adequate and
accurate for classical problems. This suggests that it should perform well in analyses of practical
interest.

3.2. Wedge example

In this example a wide shallow wedge is driven slightly into a flat plate. The flat plate is of unit
thickness and is rigidly fixed on its lower surface. The wedge varies from unit thickness at the outer
edges to a thickness of 1.3 at the center. Initially the point of the wedge is just touching the plate. The
central portion of the wedge has an imposed downward motion of 0.3. The problem is symmetric with
respect to the central plane, so only the right half is modelled. In five runs the geometry and fixed

Fig. 4. Hertz example: initial configuration

Fig. 5. Hertz example: final configuration of contact region.

 W.R.C. Underhill et al. I Comput. Methods Appi. Mech. Engrg. 143 (1997) 229-247 243

NORMAL CONTACT STRESS

MAXlhUvl NORNALCONTACT STRESS

0.0 0.5 110

RADIAL DISTANCE

CONTACT RADIUS

Fig. 6. Hertz example: normal stresses in contact region.

displacements stay the same while the elastic modulus of the materials is varied. This will show whether
there is any practical limit on the relative stiffness of the two materials. The five cases examined are:
(1) The wedge is 1000 times stiffer than the plate.
(2) The wedge is twice as stiff as the plate.
(3) The wedge and the plate are of equal stiffness.
(4) The plate is twice as stiff as the wedge.
(5) The plate is 1000 times as stiff as the wedge.
The initial configuration and mesh is shown in Fig. 7. For the sake of visual comparison Fig. 8 shows the
cases 1, 3 and 5; this allows comparison of the extreme cases and the equal stiffness case, Fig. 9 shows

Fig. 7. Wedge example: initial configuration.

Fig. 8. Wedge example: cases 1, 3 and 5: extreme cases.

CASE 3

CASE 2

CASE 1

Fig. 9. Wedge example: cases 1, 2 and 3: wedge stiffer than plate.

Fig. 10. Wedge example: cases 3, 4 and 5: plate stiffer than wedge.
244 W.K.C. Underhill et ul. i C‘omput. Methods Appl. Mech. Engrg. 143 (1997) 229-247

cases 1, 2 and 3; this compares cases where the wedge is as stiff or stiffer than the plate, Fig. 10 shows
cases 3, 4 and 5; these are the cases where the plate is at least as stiff as the wedge. The cases shown in
this example show sensible results. From them it may be inferred that the algorithm behaves well in the
presence of a wide range of stiffness ratios.

3.3. Block example

In the third example a larger block is forced against a smaller, fixed block. The materials in the
blocks are identical. The challenge to a contact calculation that this example presents is the sharp
corners of the smaller block. How well are the contact conditions met near such an irregular boundary?
The initial configuration is shown in Fig. 11, the linear prediction in Fig. 12 and the final in Fig. 13. The
model used is very crude. Even so it can be seen that reasonable results are produced in the presence of
sudden angular changes in the boundary.

3.4. Gasket example

This example demonstrates that the method behaves reversibly in the absence of other non-
linearities. The bodies in contact are a soft gasket of semicircular section and a similar stiff gasket
against which the soft one is pushed. The back side of each gasket has all of its displacements
prescribed. Both gaskets have unit radius and a thickness of 0.1. There is a vertical offset of 0.25
between the centers and a horizontal separation of 0.2. They are forced together by prescribed
horizontal motion of their backs such that each moves 0.25. This would cause considerable overlap if
contact were ignored. The initial and final positions that would be obtained by ignoring contact, are
shown in Fig. 14. The calculation, including contact, is advanced by small steps. At step 5 overlap
occurs for the first time. After iterating contact is resolved as shown in Fig. 15. At the maximum
deflection the resolved positions are shown in Fig. 16. Then the gaskets start to retreat from each other.
The overlay of the initial and final positions is very close as shown in Fig. 17. This demonstrates the
reversibility of the method.

Fig. 11. Block example: initial positions of blocks.

Fig. 12. Block example: linear prediction ignoring contact.

Fig. 13. Block example: final configuration of blocks

 W.R.C. Underhill et al. I Comput. Methods Appl. Mech. Engrg. 143 (1997) 229-247 245

Fig. 14. Gasket example: initial position and closest approach that would be obtained if contact were ignored.
Fig. 15. Gasket example: initial contact resolved.

li==i /
INITIAL POSITION

FINAL POSITION

Fig. 16. Gasket example: resolved positions at maximum deflection.

Fig. 17. Gasket example: initial and final configurations.

4. Discussion

This is a non-linear problem of some severity. There are certain aspects of non-linear finite element
calculations that should be mentioned: step size, mesh refinement, tolerances, convergence and
relaxation.
The first let us look at step size. For non-linear calculations to converge, the step size must usually be
small. For accurate results, the step size is often strictly limited; otherwise the calculated history of
loads and displacements may diverge from the intended problem to an unacceptable extent. Further,
there are cost penalties as well as accuracy disadvantages to the use of large steps. It is common for the
use of larger steps to require disproportionately more iterations. Given that the method is stable, there
is a certain problem-dependent limit, within which any step will converge quickly. It is not clear a priori
where this limit is. There is a necessary conceptual limit in these contact calculations in particular. One
of the requisites is that the penetration or emergence of nodes must be shallow. This is necessary both
to keep the step size small and to keep the direction of the normal to the surface well defined. The
limitation that has been used in the calculations presented here is that no node may penetrate all the
way through the layer of elements that form the outer rind of the body. If a node is found embedded in
an element that is in the bulk of the body’s volume, then an error is declared and calculations are
stopped. This puts a limit on the size of step that can be undertaken in the linear prediction. If the
elements on the surface are thin, then the steps must be kept very small. This is of practical importance.
 246 W.H.C‘. Underhill ct rd. i Cbmnput. Methods Appl. Mrch. Engrg. 143 (1997) “Y-247

For such an analysis, it is possible that the requirement that nodes not penetrate past the surface layer
of elements may place a more restrictive upper bound on the size of a step than the limitations imposed
by the need to define the outward normal or the desire for rapid convergence.
Mesh refinement is a matter for thought in designing a finite element model. How well any algorithm
can follow the details of contact in a small region is limited by the fineness of the mesh. Also, the
relative refinement of the mesh on either side of the contact surface is not unimportant to this method.
If on one side the mesh were very fine and on the other side the mesh were very crude, then there
would be implications for the numerical integrations of the virtual elements. According to the smaller
elements, there would be very slow variation in the opposing contact traction as compared to the
internal coordinates and low-order integration would suffice. For the larger elements, the opposing
contact tractions may seem to vary rapidly. For the large elements, it may be necessary to use an
integration scheme of quite high order. Good results are observed to occur when the elements on either
side of contact are similar in size. A size ratio of 3:l is easily tolerated. For larger ratios it is suggested
that care be taken that the surface integrals are not of too low an order. This is not an important
calculational burden, but it does affect performance.
The method uses two tolerances. If a node has penetrated, or emerged, only a very small distance,
then this is considered a negligible error. No compatibility correction is found for it. Instead, its
displacements are considered free variables in the backsubstitution pass that calculates the forces
needed to move the body to the contact surface. The other tolerance is on the equilibrium errors at
pairs of real and virtual nodes. If this error is sufficiently small, then further corrections are not
attempted. The behaviour of the method is insensitive to the equilibrium tolerance so long as it is small
compared to the nodal loads. There is an interaction between element size and penetration tolerance.
This tolerance must always be small compared to the thickness of an element or else the linear
prediction stage will tend to cause nodes to penetrate past the surface layer of elements. If very detailed
results are needed, then this tolerance can be set so low as to enforce near perfect compatibility.
The criterion for convergence is another matter that affects non-linear calculations strongly. AS used
in these calculations, a step is considered to have converged when either of two criteria is met. The first
possibility is that the corrective forces for an iteration are less than some specified fraction of the largest
for the step. The other possibility is that the displacements that resulted from applying the corrective
forces are small when compared to the largest incremental displacement for the step. If the force
criterion causes iterations to cease, then that increment of force is not applied. Experience so far
suggests that 1% gives good results. It is possible to get results with a less stringent convergence limit,
but it is not guaranteed that the contact conditions will be as well met. This is particularly true where
the bodies involved have greatly different elastic moduli and the stiffer body has sharp irregularities on
a scale similar to the size of the elements of the body with the more flexible material. In that case, it is
possible that the sharp points of the stiff material will not be completely ejected from the flexible
material. For a very approximate analysis this may not matter. For an investigation that requires
accuracy in the fine details, it is better to reduce the size of the elements in the flexible material and to
set the limits of convergence to be very fine.
Another consideration is the use of relaxation factors. Experiments show that for many analyses it is
effective to under-relax the first few iterations. It is sometimes effective to over-relax later iterations.
The most useful application of relaxation factors arises in the release of nodes that are found to have
tensile contact forces. It is possible when step sizes are large for the compatibility corrections of the first
iteration to apply tensile forces to nodes near the edge of the contact region. For physical reasons these
must be released. However, it is possible, if all of the tension is released at once, that there may be an
unwonted disturbance to the calculations. If, instead, the release of this excess traction is done over two
or more iterations, then the magnitude of subsequent corrections is much less and the overall
convergence is greatly accelerated. Trials to date suggest that releasing half of the force over each of
two iterations results in disturbances which are easily tolerated by the calculation.
It has been shown that the method is accurate, that it behaves well under a wide variety of conditions
of relative stiffness, that geometric irregularities do not disturb it unduly and that it is reversible. This
algorithm is sufficiently robust and rapidly convergent to be very suitable for practical applications.
 W.R.C. Underhill et al. I Comput. Methods Appl. Mech. Engrg. 143 (1997) 229-247 247

References

[l] S.P. Timoshenko and J.N. Goodier, Theory of Elasticity, 3rd edition (McGraw-Hill, Toronto, 1970) 403-422.
[2] K.L. Johnson, Tribology group nominated lecture: One hundred years of Hertz contact, Proc. Inst. Mech. Engrs. 196 (1982)
363-378.
[3] S.K. Pascoe and J.E. Mottershead, Two new finite element contact algorithms, Comput. Struct. 32 (1989) 137-144.
[4] D. Osmont, A finite element code for the computation of the dynamic response of structures involving contact effects,
Comput. Struct. 20 (1985) 555-561.
[5] E. Zolti, A finite element procedure to time dependent contact analysis, Comput. Struct. 17 (1983) 555-561.
[6] Wang-Quan Cheng, Fu-Wei Zhu and Ji-Wei Luo, Computational finite element analysis and optimal design for multibody
contact system, Comput. Methods Appl. Mech. Engrg. 71 (1988) 31-39.
[7] G. Yagawa and H. Hirayama, A finite element method for contact problems related to fracture mechanics, Int. J. Numer.
Methods Engrg. 20 (1984) 2175-2195.
[8] A. Francavilla and O.C. Zienkiewicz, A note on numerical computation of elastic contact problems, Int. J. Numer. Methods
Engrg. 9 (1975) 913-924.
[9] T.D. Sachdeva and C.V. Ramakrishnan, A finite element solution for the two-dimensional elastic contact problems with
friction, Int. J. Numer. Methods Engrg. 17 (1981) 1257-1271.
[lo] B. Torstenfelt, Contact problems with friction in general purpose finite element computer programs, Comput. Struct. 16
(1983) 487-493.
[ll] R.B. Haber, A mixed Eulerian-Lagrangian displacement model for large-deformation analysis in solid mechanics, Comput.
Methods Appl. Mech. Engrg. 43 (1984) 277-292.
[12] R.B. Haber and B.H. Harandria, Computational strategies for nonlinear and fracture mechanics problems, Comput. Struct.
20 (1985) 193-201.
[13] R.J. Gu, Moving finite element analysis for the elastic beams in contact problems, Comput. Struct. 24 (1986) 571-579.
[14] S.C. Shyu, T.Y. Chang and F. Salle, Friction-contact analysis using a mixed finite element method, Comput. Struct. 32
(1989) 223-242.
[15] T.Y. Chang, A.F. Saleeb and S.C. Shyu, Finite element solutions two-dimensional contact problems based on a consistent
mixed formulation, Comput. Struct. 27 (1987) 455-466.
[I61 J.C. Simo, P. Wriggers and R.L. Taylor, A perturbed Lagrangian formulation for the finite element solution of contact
problems, Comput. Methods Appl. Mech. Engrg. 50 (1985) 163-180.
[17] J. Padovan, S. Tovichakchaikul and I. Zeid, Finite element analysis of steadily moving contact fields, Comput. Struct. 18
(1984) 191-200.
[I81 W. Zhong and S. Sun, A finite element method for elasto-plastic structures and contact problems by parametric quadratic
programming, Int. J. Numer. Methods Engrg. 26 (1988) 2723-2738.
[I91 M.E. Plesha, R. Ballarini and A. Parulekar, Constitutive model and finite element procedure for dilatant contact problems,
J. Engrg. Mech. 115 (1989) 2649-2668.