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

HAL
open science

# An introduction to the unilateral contact problems in mechanics

Raul Feijoo, Helio Barbosa

► To cite this version:

Raul Feijoo, Helio Barbosa. An introduction to the unilateral contact problems in mechanics. Latin American VI Congreso on Methods for Computacionales Ingenieria, Oct 1985, Santa Fé, Argentina. hal-01295252

HAL Id: hal-01295252

https://hal.science/hal-01295252

Submitted on 1 Apr 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_2.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
AN INTRODUCTION TO THE UNILATERAL CONTACT
PROBLEMS IN MECHANICS

Raúl A. Feijão

Hélio J.C. Barbosa

Laboratório Nacional de Computação Científica
Rua Lauro Müller 455
22290 Rio de Janeiro, RJ
Brazil

RESUMEN

Se presenta en este trabajo el problema del contacto unilateral den-
tro de la teoría lineal de la elastoestática. El problema es formulado
como una inecuación variacional o como el problema de la minimización de
un funcional definido en un convexo. La obtención de soluciones aproxima
das es realizada a través del método de los elementos finitos y técnicas
de programación matemática. Por último son mostrados algunos ejemplos nu-
méricos.

ABSTRACT

In this paper the unilateral contact problem within the elastostatic
linear theory is presented. The problem is formulated as a variational
inequality or as a constrained minimization problem. The accomplishment
of approximate solutions is done using the finite element method and
mathematical programming techniques. At last some numerical examples are
shown.
# INTRODUCTION

In the several areas of the engineering sciences are frequently found structures which are in contact with their supports (or other components) but without being perfectly bounded to it. In order to emphasize the possibility that the structure will lose contact with the support and/or the possibility of slippage are not excluded one uses the expression *unilateral contact*.

As classical examples of the several engineering problems related to the analysis and design of structures with unilateral supports, one may enumerate piping in petrochemical plants, piping in nuclear power plants, joints between mechanical components, etc.

As one will notice in this paper, the unilateral contact leads to a non-linear problem, independently of the characteristics of the materials behaviour with which the structural elements in contact are build up. As a result most work related to numerical solutions for contact problems in elasticity is done according to one of the two following approaches:

i) Incremental techniques which almost always require some sort of iterative procedure and/or the introduction of special artificial interface elements [1,2].

ii) Direct formulation based upon variational principles that lead to an optimization problem which is solved by mathematical programming techniques [3-15].

Efforts along (i) seem to be motivated by the desire to reduce development costs by introducing special elements and procedures in existing finite element computer codes for linear and non-linear analysis. However the resulting algorithms often lack a convergence proof.

In this paper the second approach will be followed since, as one shall see later on, it enables a more accurate formulation from both the mechanical and mathematical point of view, where the numerical algorithms arise in a much more natural way and where mathematical results concerning existence and uniqueness of the solution as well as convergence of numerical algorithms are available [7-19].

On this presentation, which should be considered as a simple introduction to the unilateral contact problem, one will analyse this problem in the context of the classical elasticity theory (displacements and infinitesimal deformations) and where only equilibrium problems in terms of displacements will be presented. The formulation in terms of stresses or mixed formulations [3,4,5] and dynamic problems [7,12] will not be considered in this introduction.

In order to establish the variational formulation which one has referred to, one uses the principle of virtual work where the unilateral characteristic of the kinematical restriction will be considered.

Therefore the equilibrium will be characterized by a variational inequality instead of a variational equation (as proposed by the other formulations (i)) which, for the type of material that has been adopted, is equivalent to the minimization of a functional defined in a convex
set.

As naturally suggested by this variational formulation the numerical algorithm to obtain approximate solutions will consist in the redefinition of this problem in a subspace of finite dimension. This will lead to a mathematical programming problem.

Historically this problem has been formulated by Signorini in 1933 [6], but it was recently from 1970 on, that the problem was deeply studied from the mechanical and mathematical points of view. In the works [7-15] the reader will also find a vast bibliography.

# PART 1. FRICTIONLESS UNILATERAL CONTACT PROBLEMS IN THE ELASTOSTATIC STRUCTURAL ANALYSIS

## THE EQUILIBRIUM PROBLEM

Consider a body which in its undeformed state occupies an open, bounded, connected subset of the three-dimensional Euclidean point-space, Figure 1. The boundary $\partial\Omega$ of $\Omega$ is assumed to be regular and it consists of three open, disjoint parts $\partial\Omega_u$, $\partial\Omega_f$ and $\partial\Omega_c$, i.e.:

$$ \partial\Omega = \partial\Omega_u \cup \partial\Omega_f \cup \partial\Omega_c $$

$$ \partial\Omega_u \cap \partial\Omega_f = \partial\Omega_u \cap \partial\Omega_c = \partial\Omega_f \cap \partial\Omega_c = \emptyset \quad (1) $$

$\emptyset$: empty set

<figure><img src="image_3.png" /><figcaption>Figure 1</figcaption></figure>

On $\partial\Omega_u$ the displacements are prescribed and will be considered nule in order to simplify the presentation:

$$ u = 0 \quad \text{on} \quad \partial\Omega_u $$

On $\partial\Omega_f$ the surface forces are also prescribed and characterized by the vector value function $a$. On $\partial\Omega_c$ the body is assumed to be supported without friction by a rigid unilateral support S.

Since the contact is unilateral, the actual displacement $u$ of the
body when submitted to the action of the surface loads a and body forces b, will be such that:

$$ u \cdot n = u_n \le 0 \qquad \psi_x e^{\partial \Omega_c} \qquad (2) $$

where n is the unit outward normal vector to the candidate contact boundary $\partial\Omega_c$. On the other hand at those points of $\partial\Omega_c$ which remain in contact with the support the reaction will be such that:

$$ \lambda_n = \lambda_n u_n \qquad \psi_x e^{\partial \Omega_c} \qquad (3) $$

since it is supposed that the friction is nule.

Due to the unilateral nature of the contact, the reaction points away from the support. Hence:

$$ \lambda_n \le 0 \quad \text{if} \quad u_n = 0 \qquad (4) $$

$$ \lambda_n = 0 \quad \text{if} \quad u_n < 0 $$

The expressions (2) and (4) can be rewritten in the following way:

$$ \lambda_n \le 0, \quad u_n \le 0, \quad \lambda_n u_n = 0 \qquad \psi_x e^{\partial \Omega_c} \qquad (5) $$

known as the *complementarity condition*.

The elastostatic equilibrium problem therefore will be:

P1) Find $u$ sufficiently regular so that satisfy:

- Equilibrium equations

$$ \mathrm{div}(\mathcal{D}E(u)) + b = 0 \quad \text{in } \Omega $$

- Boundary conditions

$$ \mathcal{D}E(u)n = a \quad \text{on} \quad \partial\Omega_f $$

$$ u = 0 \quad \text{on} \quad \partial\Omega_u $$

$$ \lambda_n = \mathcal{D}E(u)n \cdot n \le 0, \quad u_n = u \cdot n \le 0, \quad \lambda_n u_n = 0 \quad \text{on} \quad \partial\Omega_c $$

where:

$$ \mathcal{E}(\cdot) = \frac{1}{2} (\nabla(\cdot) + \nabla(\cdot))^T : \text{ is the strain operator} $$

$\nabla$ : is the gradient operator

and where $D$ is the fourth-order elasticity tensor field which satisfies the usual properties of symmetry and ellipticity:

$$ D(x)A \cdot B = A \cdot D(x) B $$
where $\Omega$ is a constant.

Let $D(x)A \cdot A \geq \alpha A \cdot A$ be a second order tensor. We define $\Psi x \in \Omega$, $\Psi A, B \in \text{Sym}$, $\text{Sym}$ is the space of symmetric 2nd order tensors.

Comparing (P1) with the classical elastostatic problem of equilibrium (without unilateral restrictions) one can observe the difficulty which rises from the fact that $\partial\Omega_c$ is in itself unknown. In other words the solution of (P1) should also give us the contact region.

## VARIATIONAL FORMULATION OF THE EQUILIBRIUM PROBLEM

Instead of considering the problem (P1) it will be shown its equivalence to the variational formulation which is the principle of virtual work.

Therefore, one will consider the set:

$$ K = \{v; v \text{ sufficiently regular in } \Omega, v=0 \text{ on } \partial\Omega_u, v_n = v \cdot n \le 0 \text{ on } \partial\Omega_c \} $$

where one vaguely assumes that regularity is sufficient to render meaningful, at least in some generalized sense, the operations introduced below.

From the definition of $K$ it follows that:

$$ v \equiv 0 \in K $$

$$ \text{if } v \in K \rightarrow \lambda v \in K \quad \forall \lambda \ge 0 $$

$$ \text{if } v_1, v_2 \in K \rightarrow v = \theta v_1 + (1-\theta)v_2 \in K \quad \forall \theta \in [0,1] $$

hence $K$ is a convex cone and the solution $u$ of the problem (P1) also belongs to the set $K$:

$$ u \in K $$

(6)

From the mechanical point of view $K$ is the set of all kinematical admissible displacement fields.

Taking $v \in K$ arbitrarily one can multiply (P1)$_1$ by $v$ and integrating over $\Omega$ one has:

$$ \int_{\Omega} \mathrm{div}(\mathcal{D}E(u)) \cdot v d\Omega + \int_{\Omega} b \cdot v \, d\Omega = 0 \quad \forall v \in K $$

From the divergence theorem and from the fact that $u$ is the solution of (P1), the last expression gives:

$$ \int_{\Omega} \mathcal{D}E(u) \cdot E(v) d\Omega = \int_{\Omega} b \cdot v \, d\Omega + \int_{\partial\Omega_c} a \cdot v \, d\partial\Omega + \int_{\partial\Omega_c} \lambda \cdot v \, d\partial\Omega \quad \forall v \in K \quad (7) $$

The difficulty remains since $\partial\Omega_c$ is not known and therefore one can
not evaluate the integral on $\partial\Omega_c$ of the expression (7). Observing that:

$$ \lambda_n \le 0 \text{ and } \nu_n \le 0 $$

one will have:

$$ \lambda \cdot v = \lambda_n \nu_n \ge 0 \quad \forall x \in \partial\Omega_c $$

and (7) can be substituted by the variational inequality:

$$ \int_{\Omega} \mathcal{D}E(u) \cdot \mathcal{E}(v) d\Omega \ge \int_{\Omega} b \cdot v \, d\Omega + \int_{\partial\Omega_f} a \cdot v \, d\partial\Omega \quad \forall v \in \mathcal{K} \quad (8) $$

In particular if $v=u$ the (8) is transformed in equality due to the existing complementarity condition (5) between $\lambda$ and $u$:

$$ \int_{\Omega} \mathcal{D}E(u) \cdot \mathcal{E}(u) d\Omega = \int_{\Omega} b \cdot u \, d\Omega + \int_{\partial\Omega_f} a \cdot u \, d\partial\Omega $$

and from this expression and (8) one arrives at:

$$ \text{P2)} \quad \int_{\Omega} \mathcal{D}E(u) \cdot \mathcal{E}(v-u) d\Omega \ge \int_{\Omega} b \cdot (v-u) d\Omega + \int_{\partial\Omega_f} a \cdot (v-u) d\partial\Omega \quad \forall v \in \mathcal{K} \quad (9) $$

The above result shows that the solution $u$ of (P1) satisfies (P2). One now will show that if $u$ satisfies (P2), therefore $u$ is the solution for (P1). In order to prove this one can apply the divergence theorem to (P2):

$$ \begin{aligned} & -\int_{\Omega} [\operatorname{div}(\mathcal{D}E(u))+b] \cdot (v-u) d\Omega + \int_{\partial\Omega_f} \mathcal{D}E(u)n \cdot a \cdot (v-u) d\partial\Omega + \\ & \qquad + \int_{\partial\Omega_c} \mathcal{D}E(u)n \cdot (v-u) d\partial\Omega \ge 0 \quad \forall v \in \mathcal{K} \end{aligned} \quad (10) $$

Taking:

$$ v = u + w, \quad w \in C_o(\Omega) = \{h; h \text{ suff., reg., } h=0 \text{ on } \partial\Omega\} $$

then $v \in \mathcal{K}$ and (10) is reduced to:

$$ -\int_{\Omega} [\operatorname{div}(\mathcal{D}E(u))+b] \cdot w d\Omega \ge 0 \quad \forall w \in C_o(\Omega) $$

As $C_o(\Omega)$ is a vector space, the previous inequality will occur for both $w$ and $-w$, therefore:

$$ \int_{\Omega} [\operatorname{div}(\mathcal{D}E(u))+b] \cdot w d\Omega = 0 \quad \forall w \in C_o(\Omega) $$
hence:

$$ \mathrm{div}(\mathcal{D}\mathcal{E}(u)) + b = 0 \quad \text{in } \Omega \quad (11) $$

An analogous reasoning to the previous one made but now adopting fields $w \in C_1(\Omega) = \{w; w \text{ suff. reg.}, w=0 \text{ on } \partial\Omega_u \text{ and } \partial\Omega_c\}$ will lead to:

$$ \mathcal{D}\mathcal{E}(u)n = a \quad \text{on } \partial\Omega_f \quad (12) $$

From the previous results (11) and (12), the expression (10) is reduced to:

$$ \int_{\partial \Omega_c} \mathcal{D}\mathcal{E}(u)n \cdot (v-u) d\partial \Omega \geq 0 \quad \forall v \in \mathbb{K} \quad (13) $$

Let $(\cdot)_t$ and $(\cdot)_n$ be the respective tangential and normal components of any vector field on $\partial\Omega_c$, taking:

$$ v = u + w \quad ; \quad w \in C_2(\Omega) = \{w; w \text{ suff. reg.}, w_n = 0 \text{ on } \partial\Omega_c\} $$

the previous expression leads to:

$$ \int_{\partial \Omega_c} \mathcal{D}\mathcal{E}(u)n_t v_t d\partial \Omega \geq 0 \quad \forall v_t \in C_2(\Omega) $$

Again, as $v_t$ are elements of the space $C_2(\Omega)$, one obtains

$$ \mathcal{D}\mathcal{E}(u)n_t = 0 \quad \text{on } \partial\Omega_c $$

which tells us that the reaction $\lambda$ associated to the unilateral kinematical restrictions on $\partial\Omega_c$ must be normal to the boundary, i.e.:

$$ \lambda = \lambda_n = \mathcal{D}\mathcal{E}(u)n \quad \text{on } \partial\Omega_c \quad (14) $$

and (13) is reduced to:

$$ \int_{\partial \Omega_c} \lambda_n (v_n - u_n) d\partial \Omega \geq 0 \quad \forall v_n \in \mathbb{K} \quad (15) $$

Remembering that $v \in \mathbb{K}$ implies that $v_n \le 0$ and $\lambda v \in \mathbb{K}$, $\forall \lambda \ge 0$, the above expression leads to:

$$ \int_{\partial \Omega_c} \lambda_n (\lambda v_n - u_n) d\partial \Omega \geq 0 \quad \forall v_n \in \mathbb{K}, \forall \lambda \geq 0 \quad (16) $$

Considering $\lambda$ arbitrarily large, the expression $\lambda v_n - u_n$ has values arbitrarily negatives, therefore (16) leads to:

$$ \lambda_n \le 0 \quad \text{on } \partial\Omega_c \quad (17) $$

Finally, considering at (15) $v=0$ and $y=2u$ one will have:
$$ \int_{\partial \Omega_c} \lambda_n u_n d\partial \Omega = 0 $$

Because $\lambda_n \le 0$ (eq. 17) and $u_n \le 0$ on $\partial \Omega_c$ (u e $\mathbb{K}$), the previous expression leads to:

$$ \lambda_n u_n = 0 \quad \text{on } \partial \Omega_c \qquad (18) $$

The results (11), (12), (14), (17) and (18) show that the solution $u$ of the problem (P2) is also the solution of the problem (P1).

Therefore one arrives at the variational characterization of the contact problem without friction:

P2) $u$ is the solution for the elastostatic equilibrium problem with frictionless unilateral contact kinematical restrictions if and only if $u$ is the solution of the following variational inequality problem: find $u \in \mathbb{K}$ such that:

$$ \int_{\Omega} \mathcal{D}E(u) \cdot E(v-u) d\Omega \geq \int_{\Omega} b \cdot (v-u) d\Omega + \int_{\partial \Omega_f} a \cdot (v-u) d\partial \Omega \quad \forall v \in \mathbb{K} $$

The reader can observe that (P2) is nothing else but the Principle of Virtual Work extended to the case of unilateral contact kinematical restrictions. It is also important to observe that (P2) is a non-linear problem due to the characterization of $\mathbb{K}$ and therefore (P2) is out of the classical elasticity domain.

For the elastic material we adopted one can define the potential energy function $\Phi$:

$$ \Phi(u) = \frac{1}{2} \int_{\Omega} \mathcal{D}E(u) \cdot E(u) d\Omega $$

which is convex:

$$ \Phi(v) - \Phi(u) \geq \int_{\Omega} \mathcal{D}E(u) \cdot E(v-u) d\Omega $$

and where the equality holds if and only if $v-u \in N(E)$, i.e. iff $v-u$ is a rigid body motion. Substituting in (P2) one arrives at the equivalent constrained minimization problem:

$$ (P3) \quad F(u) = \min\{F(v); v \in \mathbb{K}\} $$

$$ F(v) = \Phi(v) - \int_{\Omega} b \cdot v \, d\Omega - \int_{\partial \Omega_f} a \cdot v \, d\partial \Omega $$

Functional $F$ expresses the potential energy of the body and is the sum of the elastic energy and the energy of the external loads. Proposition (P3) may be regarded as a generalization of the Minimum Potential Energy Principle of classical elasticity.
## EXISTENCE AND UNICITY OF THE SOLUTION

In this section some results which are well known [7,8,10,11,12, 15,16,28] and establish the existence and unicity of the solution of (P3) will be presented.

Following J.T. Oden et al. [27] one will consider:

V a real Hilbert space

F: V→R a functional defined on V

K a non empty closed subset of V

Then, there exists a unique solution u ∈ K of the problem:

$$F(u) = \min\{F(v); v \in K\}$$

whenever the following four conditions hold:

1. K is convex

2. F is strictly convex; i.e. for θ ∈ (0,1) and u≠v

$$F(\theta u + (1-\theta)v) < \theta F(u) + (1-\theta)F(v)$$

3. F is differentiable on K; i.e. for each u ∈ K there exists an operator DF(u): V→V' such that

$$\lim_{\alpha \to 0^+} \frac{\partial F(u+\alpha v)}{\partial \alpha} = \langle DF(u), v \rangle \quad \forall v \in V$$

where V' is the dual space of V, and <·, ·> denotes duality pairing on V'×V. $\langle DF(u), v \rangle$ is the "first variation of F at u in the direction v".

4. F is coercive, i.e. for v ∈ K

$$F(v) \to +\infty \quad \text{as} \quad \|v\| \to +\infty$$

where $\| \cdot \|$ is the norm on V

and the solution u can also be characterized as the solution of the variational inequality:

$$\langle DF(u), v-u \rangle \ge 0 \quad \forall v \in K$$

For the unilateral contact problem the functional F is given by:

$$F(u) = \frac{1}{2} a(v,v) - l(v)$$

where:

$$a(v,v) = \int_{\Omega} DE(v) \cdot E(v) d\Omega$$
$$ \ell(v) = \int_{\Omega} b \cdot v \, d\Omega + \int_{\partial\Omega_f} a \cdot v \, d\partial\Omega $$

Then, as a consequence of the Korn's inequality [12] and from the properties of the elasticity tensor the bilinear symmetric form $a(\cdot, \cdot)$: $\nabla x \to \mathbb{R}$ is continuous and $V$-elliptic, i.e. there exists positive constants $M$ and $m$ such that:

$$ a(u,v) \leq M \|u\|_{V} \|v\|_{V} \quad \forall u, v \in V $$

$$ a(u,u) \geq m \|u\|_{V}^{2} \quad \forall u \in V $$

where:

$$ V = \{v; v \in (H^1(\Omega))^3, v=0 \text{ on } \partial\Omega_u\} $$

Also, $\ell$ is a continuous linear functional on $V$:

$$ \ell(v) \leq \| \ell \|_{V} \| v \|_{V} $$

Taking into account all these properties, one verifies that $F$ is strictly convex, differentiable and coercive. Moreover, since:

$$ K = \{v; v \in V; v_n \leq 0 \text{ on } \partial\Omega_c\} $$

is a nonempty closed convex subset of $V$ one follows that there exists a unique solution $u$ of the problem (P3).

## EXAMPLES OF FRICTIONLESS UNILATERAL CONTACT PROBLEMS

### The Beam Problem

Consider the plane bending of a beam with discrete frictionless unilateral supports. An elastic material is assumed as well as infinitesimal strains (and displacements) and the hypothesis that plane sections remain plane and normal to the axis of the beam after deformation.

Reference will be made to the beam schematically shown in Figure 2, which is subjected to the loading system $\ell$ comprised of a distributed loading $q$, concentrated forces $F_i$ and concentrated moments $m_i$. At point A a frictionless unilateral support initially not in contact with the beam is assumed.

<figure><img src="image_4.png" /><figcaption>Figure 2</figcaption></figure>
The first step in the formulation is to define the space V of all possible displacements and its subset K of all admissible displacements.

w is said to be a possible transversal displacement, i.e. w ∈ V, if w is sufficient regular so that one can calculate the strains associated to w. From mechanical point of view w must be continuous (otherwise the beam breaks) with continuous first derivative (rotation of the cross-section) w'-dw/dx (plastic hinges are not allowed since the material is assumed to be elastic) and the second derivative (curvature) must be sectionally continuous, the discontinuities occurring at those points where concentrated moments are applied or EI is discontinuous (E is the Young's modulus and I is the moment of inertia of the cross-section). Formally $V=H^2(B)$, where $H^2(B)$ is the Hilbert space of functions w such that $w,v'$ and $w''$ are square integrables in the interval B of the real line.

If $w \in V$ is such that satisfies the kinematical constraints then $w \in K$. In the case of the example shown in Figure 2 w must be zero at both ends of the interval B and $w \ge 0$ at the point A. Then:

$$K = \{w; w \in H_0^2(B) \text{ and } w(A) \ge 0\} \quad (19)$$

where:

$$H_0^2(B) = \{w; w \in H^2(B) \text{ and } w(0)=w(L)=0\}$$

If $\bar{w}$ denotes an arbitrary admissible displacement, $\bar{w} \in K$, such that $w(A)=0$, then:

$$K = \bar{w} + \text{Var}$$

where:

$$\text{Var} = \{\bar{v}; \bar{v} \in H_0^2(B), \bar{v}(A) \ge 0\}$$

is a convex cone with the vertex at the origin and K is then a linear variety of a convex cone.

The second step in the formulation consists in defining the strain operator which is given according to the Bernoulli theory by:

$$E(\cdot) = - \frac{d^2(\cdot)}{dx^2}$$

A kinematically admissible displacement w is then said to be rigid if $w \in K$ and $E(w)=0$. The set (subspace) of all $w \in V$ such that $E(w)=0$ is denoted by $N(E)$ and is characterized by the functions (equivalence class)

$$w = b + cx, \quad b, c \in \mathbb{R}$$

In the case of the example shown in Figure 2 $K \cap N(E) = \{0\}$ where 0 is the identically zero function. However one should note that for other
boundary conditions that intersection may contain displacements w other than θ.

Finally, the constitutive equation for elastic beam is given by:

$$M = EI(-\frac{d^2w}{dx^2})$$

Now the equilibrium problem of plane bending of an elastic beam with frictionless unilateral supports can be stated as follows: find u e K such that:

$$a(u,v-u) - \mathcal{L}(v-u) \le 0 \quad \forall v \in K \quad (20)$$

where:

$$a(u,v) = \int_{0}^{L} EI \frac{d^2u}{dx^2} \frac{d^2v}{dx^2} dx$$

$$\mathcal{L}(v) = \int_{0}^{L} qv \, dx + \sum_{i} P_{i}v(x_{i}) + \sum_{i} m_{i}v'(x_{i})$$

or equivalent, find u e K such that

$$F(u) = \min_{veK} F(v) \quad (21)$$

where:

$$F(v) = \frac{1}{2} a(v,v) - \mathcal{L}(v)$$

and K given by (19) in the case of Figure 2. For more general discrete unilateral supports constraining the i-th component of the displacement field u at the j-th support Pj, the convex set K can be written as:

$$K = \{v; v \in V, a_{ij} \le v_i(P_j) \le b_{ij}, P_j = 1, 2, \dots, m\} \quad (22)$$

where m is the total number of discrete supports.

## The Punch Problem

If a deformable body is indented by a rigid body of specified shape one has the so called rigid indentation or punch problem. For plane problems where the punch is loaded by forces directed along the x-axis (Figure 3) with resultants P and M, the configuration of the system can be defined by the triple (α,θ,u) where α and θ are respectively the depth of indentation and the angle of rotation of the punch and u is the displacement field of the deformable body.

Since one is dealing with infinitesimal deformations, the linearized non-interpenetration kinematical restriction can be written as:
$$ \{u(x_1, \phi(x_1)) - \begin{bmatrix} -\theta\psi(x_1) \\ \alpha + x_1\theta \end{bmatrix}\} \cdot n \le (\psi(x_1) - \phi(x_1))n_2 \quad \forall (x_1, \phi(x_1)) \in \partial\Omega_c \quad (23) $$

where $\psi$ describes the shape of the punch, $n(n_1, n_2)$ is the unit outward normal vector to the candidate contact boundary $\partial\Omega_c$, described by $\phi$, of the deformable body. In this case, the set $K$ is given by:

$$ K = \{(\alpha, \theta, u) \in \mathbb{R} \times \mathbb{R} \times \mathbb{V}; \text{ such that (23) must be satisfied and also the kinematical restriction on } \partial\Omega_u \} \quad (24) $$

and again, the equilibrium problem of the frictionless punch problem can be stated as follows:

$$ \text{Find}(\alpha, \theta, u) \in K \text{ such that } F(\alpha, \theta, u) = \min_{(\alpha^*, \theta^*, v^*) \in K} F(\alpha^*, \theta^*, v^*) \\ (\alpha^*, \theta^*, v^*) \in K \quad (25) $$

where:

$$ F(\alpha^*, \theta^*, v^*) = \int_{\Omega} \frac{1}{2} \left( DE(v^*) \cdot E(v^*) d\Omega - Pa^* - M\theta^* \right) \quad (26) $$

<figure><img src="image_5.png" /><figcaption>Figure 3</figcaption></figure>

## Contact Without Friction between Elastic Bodies

If the contact is between two bodies $B^1$ and $B^2$ then it will be assumed that a common unit normal $n_{12}$ can be defined along the candidate contact boundaries $\partial\Omega_c^1$ and $\partial\Omega_c^2$ and in such case one has:

$$ K = \{(v^1, v^2) \in \mathbb{V} \times \mathbb{V}; v^i = 0 \text{ on } \partial\Omega_u^i, (v^1 - v^2) \cdot n_{12} \le 0 \text{ on } \partial\Omega_c^i, i=1,2\} \quad (27) $$
where s is the initial gap (on the direction $n_{12}$) between the two bodies.

## APPROXIMATE SOLUTIONS

To obtain approximate solutions for the minimization problem P3), finite dimensional approximation $F_h$ and $K_h$ are constructed and the finite element method is chosen for the spacial discretization due to its generality and widespread use in computer programs. One is then led to the following quadratic programming problem defined in $\mathbb{R}^n$:

$$ \min_{u_h, \mathbf{K}_h} \{ F_h = \frac{1}{2} u_h^\top \mathbf{K}_h u_h - u_h^\top \mathbf{f}_h \} \quad (28) $$

where n is the number of degrees of freedom, $K_h$ is the standard stiffness matrix, $\mathbf{f}_h$ is the vector of equivalent nodal loads, $u_h$ is the unknown nodal displacement vector and h is the parameter associated to the mesh that will be dropped from now on for ease of notation.

One way to construct an approximation for $\mathbf{K}$, is to approximate the field v by the interpolation functions of the finite element method and then to enforce the non-interpenetration condition at the nodal points belonging to $\partial\Omega_c$. This was the technique adopted for the numerical examples presented here. As a result, for most problems, the constraint set $\mathbf{K}$ will be described by a set of m linear inequalities:

$$ Au \leq c \quad (29) $$

where A is a $m \times n$ matrix. However, the important particular case of discrete unilateral supports (22) - which are often used in piping systems - may lead to constraints that can be written as:

$$ a \leq u \leq b \quad (30) $$

In the following, various alternatives for the solution of the minimization problem are discussed beginning with the case of the constraint set (30).

### a) The constraint set $a \leq u \leq b$

This kind of constraint allows for the direct application of a very simple iterative algorithm: Gauss-Seidel with relaxation and projection (GSRP), see Glowinski et al [8] for details, which can be described as follows:

#### i) Initialization

- Choose $u^0$ admissible, i.e. $a \leq u^0 \leq b$

- Pick $w \in (0,2)$

#### ii) Iteration

For $k=0,1,2,...$ execute:

For $i=1,2,...,n$ execute:
$$ \cdot u_i^* = \{ f_i - \sum_{j=1}^{i-1} K_{ij} u_j^{k+1} - \sum_{j=i+1}^{n} K_{ij} u_j^k \} / K_{ii} $$

$$ \cdot u_i^{k+1} = P_i[(1-\omega)u_i^k + \omega u_i^*] $$

until

$$ \| u^{k+1} - u^k \| / \| u^k \| \le \epsilon $$

where $P_i[\cdot]$ is the projection operator for the interval $[a_i, b_i]$, that is:

$$ P_i[\alpha] = \min[b_i, \max(a_i, \alpha)], \quad \alpha \in \mathbb{R} $$

and $\epsilon$ is a suitable tolerance.

### b) The constraint set $Au \le c$

In this case the primal problem:

$$ \min \left\{ \frac{1}{2} u^T K u - u^T f \right\} \qquad (31) $$

$Au \le c$

is equivalent to the saddle-point problem:

$$ \min_{u} \max_{\lambda \ge 0} \left\{ \frac{1}{2} u^T K u - u^T f + (Au-c) \cdot \lambda \right\} \qquad (32) $$

where the Lagrange multiplier $\lambda \in \mathbb{R}^m$ has been introduced in order to release the constraint $Au \le c$. As the minimization over $u$ is unconstrained, it is attained by:

$$ u = K^{-1} (f - A^T \lambda) $$

when $K$ is positive-definite (no rigid motions allowed). Substituting in (32) one is led to the dual problem:

$$ \min_{\lambda \ge 0} \left\{ \frac{1}{2} \lambda^T P \lambda - \lambda^T d \right\} \qquad (33) $$

where $P$ is a $m \times m$ symmetric matrix and $d$ is a $m$-vector given by:

$$ P = AK^{-1}A^T, \quad d = AK^{-1}f - c \qquad (34) $$

The resulting quadratic programming problem has a simpler constraint set and is usually much smaller than (31) as $m$ is usually much smaller than $n$. If $\lambda^*$ is the solution of (33) the solution $u^*$ of (31) is given by:

$$ u^* = K^{-1} (f - A^T \lambda^*) \qquad (35) $$

For the solution of a quadratic programming problem (33) two
possibilities are considered here. The first is the use of the Gauss-Seidel algorithm with relaxation and projection, and the second is the use of Lemke's algorithm to solve the linear complementary problem associated to (33) (See Bazaraa and Shetty [17], Glowinski et al [8] and Cottle [18] for details about the algorithm and pivoting methods on which the Lemke's algorithm is based).

In fact, if one takes the standard quadratic programming problem:

$$ \begin{gathered} \min \frac{1}{2} u \cdot Qu - u \cdot b \\ \text{Au} \le c, \ u \ge 0 \end{gathered} \tag{36} $$

where $Q$ is a symmetric positive-semidefinite $n \times n$ matrix, $A$ is an $m \times n$ matrix of rank $m$, $c \in \mathbb{R}^m$, $u$ and $b \in \mathbb{R}^n$, and denoting the Lagrangian multiplier vectors of the constraints $\text{Au} \le c$ and $u \ge 0$ by $\lambda \in \mathbb{R}^m$ and $\rho \in \mathbb{R}^n$, respectively, and denoting the vector of slack variables by $y \in \mathbb{R}^m$ then, the Kuhn-Tucker conditions for (36) could be written as:

$$ \begin{gathered} w - Mz = q \\ w \cdot z = 0 \\ w \ge 0 \\ z \ge 0 \end{gathered} \tag{37} $$

where:

$$ M = \begin{bmatrix} 0 & -A \\ A^T & Q \end{bmatrix}, \quad q = \begin{bmatrix} c \\ -b \end{bmatrix}, \quad w = \begin{bmatrix} y \\ \rho \end{bmatrix}, \quad z = \begin{bmatrix} \lambda \\ u \end{bmatrix} $$

which is a linear complementarity problem solvable in a finite number of steps by Lemke's algorithm.

In particular, as the constraint set $\text{Au} \le c$ is absent in (33) the matrix $M$ for Lemke's algorithm would be $P$ itself and $q=-b$, $w=\rho$ and $z=u$ in this case.

Instead of solving the dual problem (33) one could think of solving the saddle-point problem given by (32) where the solution $(u^*, \lambda^*)$ must satisfy:

$$ \begin{gathered} Ku^* - f + A^T \lambda^* = 0 \\ (Au^*-c) \cdot \lambda^* = 0 \\ \lambda^* \ge 0 \\ \text{Au}^* - c \le 0 \end{gathered} $$

Uzawa's algorithm, which is quite general, can be applied here and in this case can be described as follows (see Kikuchi and Oden [9]):
i) Initialization

- Choose $\lambda^0 \ge 0$

- Find $u^0$: $K u^0 = f - A^T \lambda^0$

ii) Iteration

a. Set $\lambda^{n+1} = \max\{0, \lambda^n + \gamma(Au^n - c)\}$

b. Find $u^{n+1}$: $K^{-1}u^{n+1} = f - A^T \lambda^{n+1}$

c. Repeat a) and b) until $||\lambda^{n+1} - \lambda^n|| / ||\lambda^n|| \le \epsilon$ where $\epsilon$ is a suitable tolerance

The parameter $\gamma$ must be positive and sufficiently small and was set, for the numerical examples presented here, to 0.005 times the minimum coefficient of the diagonal of $K$.

In order to reduce the size of the problem a sub-structuring technique could be used. For this, constrained degrees of freedom denoted by $u_e$ and unconstrained ones $u_i$ are segregated and the functional $F$ is rewritten as:

$$F(u_i, u_e) = \frac{1}{2} [u_i^T u_e^T] \begin{bmatrix} K_{ii} & K_{ie} \\ K_{ie}^T & K_{ee}^T \end{bmatrix} \begin{bmatrix} u_i \\ u_e \end{bmatrix} - [u_i^T u_e^T] \begin{bmatrix} f_i \\ f_e \end{bmatrix} \quad (38)$$

As the minimization over $u_i$ is unconstrained it is attained by:

$$u_i = K_{ii}^{-1} (f_i - K_{ie} u_e) \quad (39)$$

Substituting into $F(u_i, u_e)$ one obtains the reduced primal problem:

$$\min \left\{ \frac{1}{2} u_e^T f_i - u_e^T f_i - u_e^T f_e \right\} \quad (40)$$

Subjected to the m constraints associated to $u_e$ and where:

$$K_{ee}^* = C^T K_{KC} , \quad f_e^* = C^T f_i , \quad C = \begin{bmatrix} -K_{ii}^{-1} & K_{ie} \\ I_{ee} & 0 \end{bmatrix}$$

and $I_{ee}$ is the m×m identity matrix. One could now choose one of the alternative schemes already described in order to solve the reduced primal problem.

At this point it is interesting to note that for the rigid punch problem considered here one could write:
$$ u = \begin{bmatrix} u_i \\ u_e \\ \alpha \\ \theta \end{bmatrix}, \quad K = \begin{bmatrix} K_{ii} & K_{ie} & 0 \\ K_{T} & K_{ee} & 0 \\ K_{ie} & K_{ee} & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad f = \begin{bmatrix} 0 \\ 0 \\ P \\ M \end{bmatrix} $$

As the global stiffness matrix K is not positive-definite the dual problem cannot be written as in (33). It is more convenient then to construct a reduced primal problem by condensation of all degrees of freedom not related to the contact surface and solve the reduced problem by Lemke's algorithm. The condensation process is always possible provided that the deformable body is properly restrained.

From the point of view of computer implementation, the solution of the unilateral contact problem by the Gauss-Seidel algorithm with relaxation and projection (restricted to constraint sets of the type $a \approx b$) seems to be the simplest one. However, substantial computer savings may be achieved if a condensation process is performed before the iterative phase begins. More efficient are the pivoting methods, which were developed first in the theory of Linear and Quadratic Programming and then extended to the linear complementarity problem (Lemke's algorithm). One of the most interesting properties of the Lemke's algorithm is that if gives the exact solution of the discrete problem in a finite number of steps.

The fact that the class of unilateral contact problems reported here can be associated with constrained minimization problems provides the possibility of employing the classical descent algorithms used in the minimization of functionals: steepest descent, gradient, conjugate gradient, etc. (See e.g. [8,16,17]). Another class of solution algorithms is that of penalty algorithm (see [8,9,16]).

## NUMERICAL EXAMPLES

In this section some numerical examples (reported by the authors at [20-25]) are analysed in order to show the feasibility of the preceding variational formulations and algorithms.

The first example consists of an elastic beam, schematically shown in Figure 4, which has been modelled by 8 beam element. The beam is built-in at node 1 and has unilateral supports at nodes 3,5,7 and 9. Four load cases have been considered, all of them consisting of the same vertical concentrated load P applied at nodes 2,4,6 and 8, respectively.

For each load case the deformed configuration is shown in Figure 5 and the corresponding support reactions are listed in Table 1.

<table><caption>Table 1. Support reactions for example 1</caption><thead><tr><th>Reactions</th><th>M<sub>1</sub></th><th>R<sub>1</sub></th><th>R<sub>3</sub></th><th>R<sub>5</sub></th><th>R<sub>7</sub></th><th>R<sub>9</sub></th></tr></thead><tbody><tr><td>Load Case 1</td><td>0.966</td><td>0.755</td><td>0.228</td><td>0.018</td><td>0.0</td><td>-0.001</td></tr><tr><td>Load Case 2</td><td>0.141</td><td>0.012</td><td>0.453</td><td>0.565</td><td>0.0</td><td>-0.029</td></tr><tr><td>Load Case 3</td><td>-0.138</td><td>-0.052</td><td>0.0</td><td>0.591</td><td>0.490</td><td>-0.030</td></tr><tr><td>Load Case 4</td><td>-0.115</td><td>-0.029</td><td>0.0</td><td>0.0</td><td>0.587</td><td>0.442</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 4</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5</figcaption></figure>
All these results agree with the exact solution and correspond to the solution of the primal problem by the GSRP algorithm and the termination criterion adopted was:

$$ \max_i \frac{|u_i^k - u_{i-1}^{k+1}|}{\max_i |u_i^k|} < 0.0001 $$

The relationship between the number of iterations and the overrelaxation parameter w is also shown in Figure 5.

The second example is that of the tridimensional piping shown in Figure 6 which was subjected to a temperature increase of 800°F. The structure is built-in at nodes 1 and 8 and has unilateral supports (with gaps) at nodes 4 and 6. The support reactions are listed in Table 2 and correspond to w=1.85 and a tolerance ε=0.00001. The number of iterations performed was 197.

<table><caption>Table 2. Support reactions for example 2</caption><thead><tr><th>Node</th><th>R<sub>x</sub></th><th>R<sub>y</sub></th><th>R<sub>z</sub></th><th>M<sub>x</sub></th><th>M<sub>y</sub></th><th>M<sub>z</sub></th></tr></thead><tbody><tr><td>1</td><td>4036</td><td>4371</td><td>334</td><td>37648</td><td>-34400</td><td>-117105</td></tr><tr><td>4</td><td>0</td><td>-5087</td><td>2218</td><td>0</td><td>0</td><td>0</td></tr><tr><td>6</td><td>-2329</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>8</td><td>-1704</td><td>715</td><td>-2552</td><td>82848</td><td>131289</td><td>27075</td></tr></tbody></table>

The third example is that of a long circular cylinder resting on a rigid and frictionless horizontal support and subjected to a vertical compressive distributed load as indicated in Figure 7(a). The cylinder is analysed under uniform pressure q on the top and a state of plane-strain is assumed. The cylinder is also assumed to be homogeneous, isotropic, linearly elastic with Young's modulus E=1000 and Poisson's ratio v=0.3 and has a radius R=8. The discrete model adopted is that shown in Figure 7(a) where 136 four node isoparametric finite elements are used resulting in 304 degrees of freedom. Four load cases were considered corresponding to the distributed loadings of q=3.75, q=6.25, q=12.5 and q=30.

In Figure 7(b) normalized contact stresses are shown together with the results given by the classical Hertz solution [26]. Contact stresses were calculated by averaging element nodal stresses obtained by solving the dual problem by Lemke's algorithm.

Similar results were obtained solving the dual problem by the GSRP algorithm with ε=0.001. The number of iterations required is shown in Table 3.

Solution of the primal problem by the GSRP algorithm with w=1.8 and ε=0.0001 required 275,266,307 and 321 iterations respectively for load cases 1 to 4. Setting ε=0.001 the number of iterations required for the solution of the reduced primal problem is shown in Table 4.
<figure><img src="image_8.png" /><figcaption>Figure 6</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7</figcaption></figure>
<table><caption>Table 3. Number of iterations required for the solution of the dual problem by the GSRP algorithm</caption><thead><tr><th colspan="5">LOAD CASE</th></tr><tr><th>w</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>0.6</td><td>17</td><td>19</td><td>20</td><td>24</td></tr><tr><td>0.8</td><td>12</td><td>16</td><td>19</td><td>26</td></tr><tr><td>1.</td><td>12</td><td>15</td><td>20</td><td>30</td></tr><tr><td>1.2</td><td>12</td><td>16</td><td>23</td><td>39</td></tr><tr><td>1.4</td><td>15</td><td>22</td><td>32</td><td>56</td></tr></tbody></table>

<table><caption>Table 4. Number of iterations required for the solution of the reduced primal problem by the GSRP algorithm</caption><thead><tr><th colspan="5">LOAD CASE</th></tr><tr><th>w</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>1.2</td><td>37</td><td>33</td><td>28</td><td>20</td></tr><tr><td>1.4</td><td>25</td><td>23</td><td>19</td><td>13</td></tr><tr><td>1.6</td><td>16</td><td>16</td><td>16</td><td>15</td></tr><tr><td>1.8</td><td>29</td><td>30</td><td>31</td><td>31</td></tr></tbody></table>

Typical total CPU time in seconds (the algorithms were implemented on FORTRAN IV, compiled without optimization and run under Michigan Terminal System-MTS-in an IBM 370/158) for the 4 load cases considered in this example were around:

34 sec. : for the solution of the dual problem by Lemke's and GSRP algorithms.

47 sec. : for the solution of the reduced primal problem by the GSRP algorithm.

573 sec. : for the solution of the primal problem by the GSRP algorithm.

Uzawa's algorithm was also applied to this example and setting $ε=0.0001$ converged with 129, 119, 109 and 109 iterations for load cases 1 to 4. The time required was about 5 times the CPU time used by Lemke's algorithm.

The fourth example is that of a circular plate of radius $R=60$ and constant thickness $h=4$ resting on a frictionless rigid horizontal foundation. The plate is subjected to its own weight and to an upward vertical load $P=100$ uniformly distributed on a disk of radius $r=8$ concentric with the plate.

The discrete model adopted corresponds to a uniform mesh of 120
axisymmetric 4-node isoparametric finite elements being 30 along the radius and 4 along the thickness (Figure 8) that results in 304 degrees of freedom.

<figure><img src="image_10.png" /><figcaption>Figure 8</figcaption></figure>

The material of the plate has a specific weight $\rho=0.8085$ and elastic constants $E=2.1 \times 10^6$ and $\nu=0.3$.

If one adopts the theory of plate bending [26] the value of the radius of the circumference inscribing the raised part takes the value $a=44$ and the deflection at the center is $u_c=1.337 \times 10^{-4}$.

For the axisymmetric solid discrete model adopted here the values found were $a=42$ (the first 2! nodes from the center are raised) and $u_c=1.347 \times 10^{-4}$, given by the solution of the dual problem by the Lemke's algorithm.

Essentially the same results are obtained after 795 iterations with $\epsilon=0.0001$ and $\omega=1$, using the GSRP algorithm.

The solution of the primal problem by the GSRP algorithm with $\epsilon=0.0001$ and $\omega=1.95$ required 1203 iterations yielding $u_c=1.276 \times 10^{-4}$ and poor results for the reactions.

Setting $\epsilon=0.0001$ and $\omega=1.9$ the solution of the reduced primal problem by the GSRP algorithm leads, after 191 iterations, to results close to those obtained by Lemke's algorithm.

Difficulties arise with the use of Uzawa's algorithm in this example. To begin with the automatic choice for the parameter which worked well in the third example did not work at all here and several trials were made before the termination test with $\epsilon=0.0001$ could be satisfied.

The values of $u_c$ corresponding to $\gamma=12000$. (1123 iterations) and $\gamma=8000$ (1747 iterations) were $u_c=1.358 \times 10^{-4}$ and $u_c=1.357 \times 10^{-4}$ respectively. However the Lagrange multipliers obtained did not agree very well with those given by Lemke's algorithm.

CPU time in seconds for all the alternatives were respectively 23, 34, 260, 34, 211 ($\gamma=12000$) and 313 ($\gamma=8000$).
The last example consists in the indentation of a rectangular block by a rigid solid with a cylindrical contact surface of radius R=8 as shown schematically in Figure 9(a). The material properties are: Young's modulus E=1000, Poisson's ratio v=0.3. A state of plane-strain is assumed and five "load" cases were considered, corresponding to prescribed values for the depth of indentation α=0.1,0.2,0.3,0.4 and α=0.5 while θ was prescribed as zero.

Due to the symmetry of the problem only one-half of the block was discretized by means of 202 four-node isoparametric finite elements resulting in 440 degrees of freedom. Figure 9(b) displays the adopted mesh before and after deformation (for α=0.8) and Figure 9(c) shows the relationship between α and the total applied force P.

Finally Figure 9(d) shows the normalized contact pressure obtained from Lagrange multipliers compared to the Hertz solution in solid lines. These results correspond to the condensation of all degrees of freedom not related to contact and solution of the reduced primal problem by Lemke's algorithm.

## PART II. UNILATERAL CONTACT PROBLEM WITH FRICTION

### THE VARIATIONAL FORMULATION

Consider again the same body taken at Part 1, but now one supposes there is an initial gaps between the body and the rigid foundation s on ∂Ω_c. Also, on ∂Ω_c friction boundary conditions are assumed to hold and in the notation of Part 1 they read:

if $|σ_t| < η|σ_n|$ then $u_t=0$ on $∂Ω_c$

if $|σ_t| = η|σ_n|$ then $∃λ≥0$ s.t. $u_t = -λσ_t$ on $∂Ω_c$

The classical formulation of the Signorini problem with friction is: find the displacement field u which satisfies the equilibrium equations and boundary conditions:

$$ \operatorname{div}(DE(u)) + b = 0 \quad \text{in } \Omega \qquad (42) $$

$$ DE(u)n = a \quad \text{on } \partial\Omega_c \qquad (43) $$

$$ u = 0 \quad \text{on } \partial\Omega_u \qquad (44) $$

$$ u \cdot n - s \le 0 $$

$$ u \cdot n - s < 0 \quad DE(u)n = 0 $$

$$ u \cdot n - s = 0 + \sigma_n(u) = DE(u)n \cdot n \le 0 \quad \text{and} \quad \text{on } \partial\Omega_c \qquad (45) $$

$$ |σ_t(u)| < η|σ_n(u)| \quad \Rightarrow \quad u_t = 0 $$

$$ |σ_t(u)| = η|σ_n(u)| \quad \Rightarrow \quad ∃λ≥0 \text{ s.t. } u_t = -λσ_t $$

In order to give the variational formulation, one defines the set of admissible displacements:

$$ K = \{ v \in V; v=0 \text{ on } \partial\Omega_u, v \cdot n - s \le 0 \text{ on } \partial\Omega_c \} $$
<figure><img src="image_11.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 9</figcaption></figure>
K is a nonempty closed convex subset of $V=(H^1(\Omega))^3$.

From eqs. (42)-(44) one obtains by means of the divergence theorem the relation:

$$ 
\begin{aligned}
a(u,v-u) = \ell(v-u) &+ \int_{\partial\Omega_c} \sigma_t(u)(v_t-u_t)d\partial\Omega + \\ 
&+ \int_{\partial\Omega_c} \sigma_n(u)(v_n-u_n)d\partial\Omega
\end{aligned} 
\quad (46) 
$$

for every $v$ such that $v=0$ on $\partial\Omega_u$.

From (45)$_{1,2,3}$ the following inequality results:

$$ \sigma_t(u)(v_t-u_t) + \eta|\sigma_n(u)|(|v_t|-|u_t|) \ge 0 \quad \text{on } \partial\Omega_c \quad (47) $$

which combined to (46) implies the variational inequality:

$$ 
\begin{aligned}
a(u,v-u) \ge \ell(v-u) &- \int_{\partial\Omega_c} \eta|\sigma_n(u)|(|v_t|-|u_t|)d\partial\Omega + \\ 
&+ \int_{\partial\Omega_c} \sigma_n(u)(v_n-u_n)d\partial\Omega
\end{aligned} 
\quad (48) 
$$

for any $v$ such that $v=0$ on $\partial\Omega_u$.

From (45)$_{1,2,3}$ and $v \in K$ one has the following variational inequality problem:

$$ a(u,v-u) + j(u,v) - j(u,u) \ge \ell(v-u) \quad \forall v \in K \quad (49) $$

where:

$$ j(u,v) = \int_{\Omega} \eta |\sigma_n(u)| |v_t| d\Omega \quad (50) $$

It can be shown [7] that the classical problem (42)-(45) is formally equivalent to the problem (49). Since the normal component of the stress density vector on $\partial\Omega_c$ is defined only as linear form, $|\sigma_n(u)|$ has no mathematical meaning.

The issue of existence and uniqueness of solutions for (49) is still open. For a particular situation Necas et al [29] showed the existence of solutions to (49) provided that $\eta$ is sufficiently small. Duvaut [30] introduced the idea of non-local friction law and established an existence result for any friction and also uniqueness for the case of small friction. Oden and Pires [31, 32, 33] proposed a class of nonlocal laws as well as numerical algorithms for obtaining approximate solutions for contact problem. See also [12] and the results reported by M. Cocu [34].

However, for the special case where the normal stress $\sigma_n$ is
prescribed along $\partial\Omega_c$, ($\sigma_n = \bar{F}_n$), Duvaut and Lions [7] established the existence and uniqueness of the solution. In this case the contact surface $\partial\Omega_c$ is known in advance and $u_n$ is not prescribed on $\partial\Omega_c$. The boundary conditions on $\partial\Omega_c$ reduce to:

$$|\sigma_t| < g + u_t = 0$$

$$|\sigma_t| = g + u_t = -\lambda\sigma_t \quad \text{for some } \lambda \ge 0$$

where $g=\eta|\bar{F}_n|$ is given and represents the maximum tangential stress that can be developed due to friction along $\partial\Omega_c$.

Introducing the functionals

$$j_g(v) = \int_{\partial\Omega_c} g|v_t| d\partial\Omega \qquad f_n(v) = \int_{\partial\Omega_c} \bar{F}_n v_n d\partial\Omega \quad (51)$$

where $g$ is a given positive function on $\partial\Omega_c$ and $\bar{F}_n$ is a given normal stress distribution on $\partial\Omega_c$ and defining the subspace $V$

$$V = \{v \in V: v|_{\Gamma_u} = 0\}$$

the PVW for the special problem of friction with prescribed normal stress can be stated as:

Find $u \in V$ such that

$$-\int_{\Omega} \mathcal{D}E(u) \cdot E(v-u)d\Omega - j_g(v) + j_g(u) + \ell(v-u) + \\ + f_n(v-u)d\Gamma \le 0 \quad \forall v \in V \quad (52)$$

It can be shown [7] that to solve this inequality is equivalent to solving the following minimization problem

$$\inf_{v \in K} [F(v) + j_g(v) - f_n(v)] \quad (53)$$

For Signorini's problem with friction, inequality (49), the following iterative procedure can be envisaged:

i) Solve Signorini's problem without friction, problem (P2) or (P3).

ii) With the normal stress in $\partial\Omega_c$ found in (i) solve the special friction problem with prescribed normal stress, inequality (52).

iii) Tangential stress found in (ii) are then used as additional loads in Signorini's problem without friction (i) and the steps (i), (ii) and (iii) are repeated until convergence is (hopefully) achieved.

Introducing the functional:
$$f_t(v) = \int_{\partial\Omega_c} \bar{F}_t \cdot v_t d\partial\Omega$$

where $\bar{F}_t$ is a given distribution of tangential forces along $\partial\Omega_c$, the procedure described above can be written as:

1) Given $\bar{F}_t^{k-1}$ find $u^k$ solution of the minimization problem

$$\inf_{veK} [F(v)-f_t(v)]$$

2) Calculate $\bar{F}_n^k = \sigma_n(u^k)$ and $g=\eta|\sigma_n(u^k)|$

3) Find $u^*$ solution of the minimization problem

$$\inf_{veV} [F(v)+j_g(v)-f_n(v)]$$

4) Calculate $\bar{F}_t^k = \sigma_t(u^k)$ and repeat all steps for $k=2,3,...$ until convergence is achieved.

The procedure just described, whose convergence has not been formally proved yet, involves two minimization problems. In the first one the main difficulty is due to the constraint set $K$ while in the second one the difficulty arises due to the non-differentiability of $j_g(v)$.

Panagiotopoulos [35] follows the scheme described above solving both minimization problems by non-linear programming techniques; Campos, Oden and Kikuchi [36] adopt a penalization technique in the first problem and a regularization technique in the second one. Some other possibilities are presented by Raous [37] and Haslinger and Panagiotopoulos [3]. The basic idea used here is duality [7]. The first minimization problem is substituted by the equivalent saddle-point problem:

$$\inf_{veV} \sup_{\lambda_n \ge 0} L_1(v, \lambda_n) \qquad (54)$$

where the constraint set $K$ is absent. The Lagrangean $L_1(v, \lambda_n)$ is given by:

$$L_1(v, \lambda_n) = F(v) - f_t(v) + \int_{\partial\Omega_c} \lambda_n (v \cdot n - s) d\partial\Omega$$

In the second minimization problem the non-differentiable functional $j_g(v)$ is replaced by:

$$\sup_{\lambda_t \in \Lambda} \int_{\partial\Omega_c} \lambda_t \cdot v_t d\partial\Omega$$

where:
$$ \Lambda = \{\lambda=(\lambda_1, \lambda_2, \lambda_3); \sum_{i=1}^{3} \lambda_i^2(x) \le g^2(x), x \in \partial\Omega_c \} $$

and one is led to the equivalent saddle-point problem:

$$ \inf_{v \in V} \sup_{\lambda_t \in \Lambda} L_2(v, \lambda_t) \qquad (55) $$

with:

$$ L_3(v, \lambda_t) = F(v) - f_n(v) + \int_{\partial \Omega_c} \lambda_t \cdot v \, d\partial \Omega $$

The Lagrange multipliers $\lambda_n$ and $\lambda_t$ can be interpreted, by duality, respectively as the normal stress on $\partial\Omega_c$ and the tangential stress, due to friction, on $\partial\Omega_c$.

## APPROXIMATE SOLUTIONS

To obtain approximate solutions for the problems formulated in the preceding section the finite element method is used to construct finite-dimensional approximation spaces. For plane problems the following interpolation scheme can be adopted:

$$ v = \phi q, \quad \lambda_n = \psi p, \quad \lambda_t = \psi t $$

where $\phi$ is the matrix of interpolation functions for the displacements field $v$ in terms of the nodal unknowns $q$ and $\psi$ is a row-vector with interpolation functions for the Lagrange multipliers $\lambda_n$ and $\lambda_t$ in terms of the parameters $p$ and $t$. The global interpolants are constructed from local bilinear interpolants associated to a four-node quadrilateral isoparametric finite element. The interpolation of $\lambda_t$ and $\lambda_n$ is done by means of piecewise constant functions along the sides of the elements on $\partial\Omega_c$. In this way problems (54) and (55) are approximated by:

$$ \min_{q \ p \ge 0} \max_{p \ge 0} \frac{1}{2} q \cdot Kq - q \cdot F - q \cdot F_t + q \cdot Mp + S \cdot p \qquad (56) $$

and

$$ \min_{q} \max_{|\psi t| \le g} \frac{1}{2} q \cdot Kq - q \cdot F - q \cdot F_n + q \cdot At \qquad (57) $$

where $K$ is the standard stiffness matrix, $F$, $F_t$ and $F_n$ are vectors of nodal loads which are equivalent, respectively, to the applied load system $\ell$, tangential loads due to friction and normal reactions in the contact surface. The matrices $M$ and $A$ and the vector $S$ are given by:

$$ M = \int_{\partial \Omega_c} \phi^T \{ \sum_{n_y}^{n} x_i \} \psi d\partial \Omega, \quad A = \int_{\partial \Omega_c} \phi^T \{ \sum_{t_y}^{t} x_i \} \psi d\partial \Omega \qquad (58) $$
and

$$S = \int_{\partial \Omega_c} \psi^T s d\partial \Omega \qquad (59)$$

where $(n_K, n_y)$ and $(t_K, t_y)$ are the components of unitary vectors respectively outward normal and tangent to the boundary $\partial\Omega_c$. It is clear that $q \in \mathbb{R}^n$, $p \in \mathbb{R}^m$ and $t \in \mathbb{R}^m$ where $n$ is the number of degrees of freedom of the discrete model defined by the finite element mesh adopted and $m$ is the number of elements along the boundary $\partial\Omega_c$.

As the minimization over $q$ is unconstrained and $K$ is assumed positive-definite this variable can be eliminated using the stationarity conditions:

$$Kq - (F+F_t) + Mp = 0 \qquad (60)$$

and

$$Kq - (F+F_n) + At = 0 \qquad (61)$$

Combining (56) with (60) and (57) with (61) one obtains:

$$\min_{p>0} \frac{1}{2} p \cdot Pp - p \cdot d_1 \qquad (62)$$

$$-\frac{\min}{g \le t \le g} \frac{1}{2} t \cdot Tt - t \cdot d_2 \qquad (63)$$

where:

$$P = M^T K^{-1} M, \quad d_1 = M^T K^{-1} (F + F_t) - S$$

$$T = A^T K^{-1} A, \quad d_2 = A^T K^{-1} (F + F_n)$$

The vector $\bar{g} \in \mathbb{R}^m$ has its i-th entry equal to $\eta$ times the absolute value of the (prescribed) normal stress along the side of the i-th element in $\partial\Omega_c$.

Due to the type of constraints that arise in problems (62) and (63) a very simple numerical algorithm can be used: Gauss-Seidel with relaxation and projection (GSRP). Finally, from (60) nodal unknowns $q$ are obtained and element stresses can be computed.

**Remark 1.** Although a piecewise constant interpolation for the Lagrange multipliers $\lambda_n$ and $\lambda_t$ has been used to obtain the finite-dimensional approximations (56) and (57) it is important to note that matrices $M$ and $A$ in (56) and (57) result from the approximation of the integrals:

$$\int_{\partial \Omega_c} \lambda_n (v \cdot n - s) d\partial \Omega \quad \text{and} \quad \int_{\partial \Omega_c} \lambda_t \cdot u_t d\partial \Omega$$
and, as such, can take different forms according to the numerical scheme adopted in the approximation of these integrals. In fact, in addition to the piecewise constant interpolation scheme already mentioned another scheme was tried which resulted in an improved performance of the numerical algorithm. The idea was to take concentrated Lagrange multipliers in the nodal points along $\partial\Omega_c$. In this case, M, A and S can still be given by (58) and (59) provided we take the entries of row-vector $\psi$ as Dirac's delta "functions" associated to the nodal points along $\Gamma_c$.

**Remark 2.** The problem of contact with Coulomb friction between two deformable bodies and the problem of indentation of a deformable body by a rigid one can both be treated along the same basic lines. Indeed work is under way in this direction and the results will be reported soon.

## A NUMERICAL EXAMPLE

This section describes the results of some numerical experiments performed with the algorithm discussed in the preceding section. The problem considered is that of a rectangular block pressed against a rigid horizontal foundation on which Coulomb's law of friction is assumed to hold. The block is also submitted to a horizontal uniformly distributed load as shown in Figure 10 and a state of plane strain is assumed. The material of the block is homogeneous and isotropic with Young's modulus $E=13000$ and Poisson's coefficient $\nu=0.2$. Due to the symmetry of the problem, only half of the block was discretized by means of 194 four-node isoparametric finite elements leading to a discrete model with 439 degrees of freedom.

The first load case considered here corresponds to $F=15$, $f=5$ and a coefficient of friction $\eta=1.0$. Figure 11 shows the deformed mesh amplified by a factor of 100. Normal and tangential nodal displacements as well as normal and tangential nodal reactions along the contact surface are displayed in Table 5 where regions with different behaviour are easily identified: adhesion (from node 1 to node 9), sliding (from node 10 to node 30) and a region where contact was lost (from node 31 to node 33).

<figure><img src="image_16.png" /><figcaption>Figure 10</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 11</figcaption></figure>
<table><caption>Table 5</caption><thead><tr><th>NODE</th><th>u<sub>n</sub></th><th>u<sub>t</sub></th><th>F<sub>n</sub></th><th>F<sub>t</sub></th></tr></thead><tbody><tr><td>1</td><td>0.</td><td>0.</td><td>5.23</td><td>0.</td></tr><tr><td>2</td><td>0.</td><td>0.</td><td>10.47</td><td>0.64</td></tr><tr><td>3</td><td>0.</td><td>0.</td><td>10.47</td><td>1.28</td></tr><tr><td>4</td><td>0.</td><td>0.</td><td>10.46</td><td>1.96</td></tr><tr><td>5</td><td>0.</td><td>0.</td><td>10.44</td><td>2.69</td></tr><tr><td>6</td><td>0.</td><td>0.</td><td>10.41</td><td>3.52</td></tr><tr><td>7</td><td>0.</td><td>0.</td><td>10.38</td><td>4.48</td></tr><tr><td>8</td><td>0.</td><td>0.</td><td>10.31</td><td>5.75</td></tr><tr><td>9</td><td>0.</td><td>0.</td><td>9.80</td><td>8.57</td></tr><tr><td>10</td><td>0.</td><td>0.0003</td><td>8.93</td><td>8.93</td></tr><tr><td>11</td><td>0.</td><td>0.0009</td><td>8.28</td><td>8.28</td></tr><tr><td>12</td><td>0.</td><td>0.0016</td><td>7.86</td><td>7.86</td></tr><tr><td>13</td><td>0.</td><td>0.0025</td><td>7.52</td><td>7.52</td></tr><tr><td>14</td><td>0.</td><td>0.0034</td><td>7.21</td><td>7.21</td></tr><tr><td>15</td><td>0.</td><td>0.0043</td><td>6.93</td><td>6.93</td></tr><tr><td>16</td><td>0.</td><td>0.0054</td><td>6.66</td><td>6.66</td></tr><tr><td>17</td><td>0.</td><td>0.0065</td><td>6.39</td><td>6.39</td></tr><tr><td>18</td><td>0.</td><td>0.0076</td><td>6.13</td><td>6.13</td></tr><tr><td>19</td><td>0.</td><td>0.0088</td><td>5.86</td><td>5.86</td></tr><tr><td>20</td><td>0.</td><td>0.0101</td><td>5.59</td><td>5.59</td></tr><tr><td>21</td><td>0.</td><td>0.0113</td><td>5.30</td><td>5.30</td></tr><tr><td>22</td><td>0.</td><td>0.0127</td><td>5.00</td><td>5.00</td></tr><tr><td>23</td><td>0.</td><td>0140</td><td>4.68</td><td>4.68</td></tr><tr><td>24</td><td>0.</td><td>0.0154</td><td>4.34</td><td>4.34</td></tr><tr><td>25</td><td>0.</td><td>0.0169</td><td>3.96</td><td>3.96</td></tr><tr><td>26</td><td>0.</td><td>0.0183</td><td>3.53</td><td>3.53</td></tr><tr><td>27</td><td>0.</td><td>0.0198</td><td>3.04</td><td>3.04</td></tr><tr><td>28</td><td>0.</td><td>0.0213</td><td>2.46</td><td>2.46</td></tr><tr><td>29</td><td>0.</td><td>0.0229</td><td>1.72</td><td>1.72</td></tr><tr><td>30</td><td>0.</td><td>0.0244</td><td>0.66</td><td>0.66</td></tr><tr><td>31</td><td>0.0001</td><td>0.0258</td><td>0.</td><td>0.</td></tr><tr><td>32</td><td>0.0003</td><td>0.0272</td><td>0.</td><td>0.</td></tr><tr><td>33</td><td>0.0006</td><td>0.0286</td><td>0.</td><td>0.</td></tr></tbody></table>
<table><caption>Table 6</caption><thead><tr><th>NODE</th><th>u<sub>t</sub></th><th>F<sub>n</sub></th><th>F<sub>t</sub></th></tr></thead><tbody><tr><td>1</td><td>0.</td><td>10.40</td><td>0.</td></tr><tr><td>2</td><td>0.</td><td>20.80</td><td>0.17</td></tr><tr><td>3</td><td>0.</td><td>20.80</td><td>0.34</td></tr><tr><td>4</td><td>0.</td><td>20.79</td><td>0.52</td></tr><tr><td>5</td><td>0.</td><td>20.77</td><td>0.70</td></tr><tr><td>6</td><td>0.</td><td>20.76</td><td>0.88</td></tr><tr><td>7</td><td>0.</td><td>20.74</td><td>1.08</td></tr><tr><td>8</td><td>0.</td><td>20.71</td><td>1.28</td></tr><tr><td>9</td><td>0.</td><td>20.68</td><td>1.50</td></tr><tr><td>10</td><td>0.</td><td>20.65</td><td>1.73</td></tr><tr><td>11</td><td>0.</td><td>20.61</td><td>1.99</td></tr><tr><td>12</td><td>0.</td><td>20.57</td><td>2.30</td></tr><tr><td>13</td><td>0.</td><td>20.52</td><td>2.66</td></tr><tr><td>14</td><td>0.</td><td>20.42</td><td>3.23</td></tr><tr><td>15</td><td>0.00003</td><td>20.17</td><td>4.03</td></tr><tr><td>16</td><td>0.00016</td><td>19.86</td><td>3.97</td></tr><tr><td>17</td><td>0.00036</td><td>19.63</td><td>3.93</td></tr><tr><td>18</td><td>0.00060</td><td>19.43</td><td>3.89</td></tr><tr><td>19</td><td>0.00088</td><td>19.24</td><td>3.85</td></tr><tr><td>20</td><td>0.00119</td><td>19.04</td><td>3.81</td></tr><tr><td>21</td><td>0.00153</td><td>18.84</td><td>3.77</td></tr><tr><td>22</td><td>0.00189</td><td>18.63</td><td>3.73</td></tr><tr><td>23</td><td>0.00228</td><td>18.40</td><td>3.68</td></tr><tr><td>24</td><td>0.00270</td><td>18.15</td><td>3.63</td></tr><tr><td>25</td><td>0.00314</td><td>17.86</td><td>3.57</td></tr><tr><td>26</td><td>0.00360</td><td>17.53</td><td>3.51</td></tr><tr><td>27</td><td>0.00409</td><td>17.14</td><td>3.43</td></tr><tr><td>28</td><td>0.00461</td><td>16.65</td><td>3.33</td></tr><tr><td>29</td><td>0.00515</td><td>16.03</td><td>3.21</td></tr><tr><td>30</td><td>0.00572</td><td>15.17</td><td>3.03</td></tr><tr><td>31</td><td>0.00633</td><td>13.88</td><td>2.78</td></tr><tr><td>32</td><td>0.00699</td><td>11.50</td><td>2.30</td></tr><tr><td>33</td><td>0.00776</td><td>3.64</td><td>0.73</td></tr></tbody></table>
A second load case, corresponding to $F=10$ and $f=15$, was analysed considering $\eta=0.2$ and the results are summarized in Table 6 where it can be seen that nodes 1 to 14 are in adhesion and nodes 15 to 33 are in a sliding condition.

Other load cases were also analysed and the results obtained agree with those found by Raous using a different algorithm [37].

## FINAL REMARKS

Concluding this introduction and following Prof. G. Del Piero's remarks it is important to emphasize that besides the unilateral problem which was seen before, there are other type of problems which are also associated to unilateral restrictions.

Among these problems which approximate solutions will be given by mathematical programming techniques one may enumerate fracture problems, problems which arise on non resisting tension materials (concrete, rocks, ceramics, soils, bricks, etc.), limited strength in tension, etc.

Problems associated to plasticity should also be mentioned. Elastic-plastic behavior is another important example of unilateral internal restrictions. Here the reader will find a wide variety of applications of mathematical programming problems [40].

Finally, the problems associated to structural optimization should be emphasized. Here the reader can observe that the numerical algorithms applied to optimization will from now on be applied to structural analysis [41].

## BIBLIOGRAPHY

[ 1 ] FRANCAVILLA, A.; ZIENKIEWICZ, O.C.; "A note on the numerical computation of elastic contact problem", Int.J.Num.Meth.Engng., vol. 9, 913-924, 1975.

[ 2 ] STADTER, J.T.; WEISS, R.D.; "Analysis of contact through finite element gaps", Computers and Structures, vol. 10, 867-873, 1979.

[ 3 ] HASLINGER, J.; PANAGIOTOPOULOS, P.D.; "The reciprocal variational approach to the Signorini problem with friction. Approximation results", Proceedings of the Royal Society of Edinburgh, 98A, 365-383, 1984.

[ 4 ] HASLINGER, J.; HLAVACEK, J.; "Approximation of the Signorini problem with friction by a mixed finite element method", J.Math. Anal.Appl., 86, 99-122, 1982.

[ 5 ] de SAXCE, G.; NGUYEN DANG HUNG; "Dual analysis of frictionless problems by displacement and equilibrium finite elements", Eng. Struct., vol. 6, 26-32, 1984.

[ 6 ] SIGNORINI, A., "Sopra alcune questioni di elastostatica", Atti Soc. It. per il Progresso delle Scienze, 1933.

[ 7 ] DUVAUT, G. LIONS, J.L., "Les inéquations en mécanique et en physique", Dunod, Paris, 1972.
[8] GLOWINSKI, R.; LIONS, J.L.; TREMOLIÈRES, R.; "Analyse numérique des inéquations variationnelles", Dunod, Paris, 1976.

[9] KIKUCHI, N.; ODEN, J.T.; "Contact problems in elasticity", TICOM Report 79-8, The Univ. of Texas at Austin, 1979.

[10] BRÉZIS, P.H.; "Problèmes unilatéraux", J.Math. Pures et Appl., 51, 1-168, 1972.

[11] KIKUCHI, N.; ODEN, J.T.; "Contact problems in elastostatics", Finite Elements. Special Problems in Solid Mechanics, Ed. J.T.Oden and G.F. Carey, vol. 5, chapter 4, 158-212, 1984.

[12] PANAGIOTOPOULOS, P.D.; "Inequality problems in Mechanics and Applications", Birkhäuser, 1985.

[13] DEL PIERO, G.; "Unilateral problems in structural analysis", Módulo III del 2º Curso de Mecânica Teórica e Aplicada, Laboratório Nacional de Computação Científica, July 1985, Brazil.

[14] DEL PIERO, G.; MACERI, F.; (Editores) Unilateral Problems in Structural Analysis, CISM Courses and Lectures, Springer, 1985.

[15] FICHERA, G.; "Boundary value problems of elasticity with unilateral constraints", Handbuch der Physik, Band via/2, 391-424, 1972.

[16] CEA, J.; "Optimisation, théorie et algorithmes", Dunod, Paris, 1971.

[17] BAZARAA, M.S.; SHETTY, C.M.; "Nonlinear programming: theory and algorithms", John Wiley and Sons, 1979.

[18] COTTLE, R.W.; "Some recent developments in linear complementarity theory", in: R.W. Cottle, F. Giannessi, and J.L. Lions (Editors), Variational Inequalities and Complementarity Problems. Theory and Applications, Wiley, 1980.

[19] LUENBERGER, D.G.; "Optimization by vector space methods", John Wiley, 1969.

[20] FEIJÓO, R.A.; BARBOSA, H.J.C.; "Análisis estático de vigas elásticas con apoyos unilaterales", Report n° 005/83, Laboratório de Computação Científica, Rio de Janeiro, Brazil.

[21] FEIJÓO, R.A.; BARBOSA, H.J.C.; "Static analysis of piping systems with unilateral supports", VII Congresso Brasileiro de Engenharia Mecânica, vol. D, 269-279, Uberlândia, dec, 1983.

[22] BARBOSA, H.J.C.; FEIJÓO, R.A.; "Numerical algorithms for contact problems in linear elastostatics", Conference on Structural Analysis and Design of Nuclear Power Plants, vol. 1, 231-244, Porto Alegre, Brazil, Oct., 1984.

[23] BARBOSA, H.J.C.; FEIJÓO, R.A.; "Um algoritmo para o problema de identação rígida em elastostática", V Congresso Latino-Americano de Métodos Computacionais para Engenharia, Salvador, Brazil, oct/nov, 1984. Report n° 025/84, LNCC/CNPq.
[24] BARBOSA, H.J.C.; FEIJÓO, R.A.; "Numerical algorithms for frictionless contact problems in linear elastostatics", 2nd *International Conference on Variational Methods in Engineering*, Southampton, England, jul, 1985.

[25] BARBOSA, H.J.C.; FEIJÓO, R.A.; "Um algoritmo numérico para o problema de Signorini com atrito de Coulomb", VIII Congresso Brasileiro de Engenharia Mecânica, São José dos Campos, Brazil, dec, 1985 (to be published).

[26] TIMOSHENKO, S.; GOODIER, J.N.; "Teoría de la elasticidad", Ed. Urmo, Bilbao, 1968.

[27] ODEN, J.T.; KIKUCHI, N.; "Finite element methods for constrained problems in elasticity", Int.J.Num.Meth.Engng., vol. 18, 701-725, 1982.

[28] ROMANO, G.; "Principi e metodi variazionali nella meccanica delle strutture e dei solidi", Parte I: Principi", Facoltà di Ingegneria dell'Università di Napoli, Italy, oct, 1984.

[29] NECAS, J.; JARUSEK, J.; HASLINGER, J.; "On the solution of the variational inequality to the Signorini problem with small friction", Boll.Un.Mat.Ital., 17B, 796-811, 1980.

[30] DUVAUT, G.; "Equilibre d'un solide elastique avec contact unilateral et frottement de Coulomb", C.R. Acad.Sc. Paris, t290, serie A, 263-265, 1980.

[31] ODEN, J.T.; PIRES, E.B.; "Numerical analysis of certain contact problems with non-classical friction laws", Computers and Structures, 16, 471-478, 1983.

[32] ODEN, J.T.; PIRES, E.B.; "Non local and non linear friction laws and variational principles for contact problems in elasticity", J. Appl.Mech., 50(1), 67-76, 1983.

[33] ODEN, J.T.; PIRES, E.B.; "Nonlocal friction in contact problems in plane elasticity", J.T. Oden and G.F. Carey (Editors), Finite Elements. Special Problems in Solid Mechanics, vol. v, 213-252, Prentice Hall, 1984.

[34] COCU, M.; "Existence of solutions of Signorini problems with friction", Int.J. of Engng.Sci., vol. 22, n9 5, 567-575, 1984.

[35] PANAGITOPOULOS, P.D.; "A nonlinear programming approach to the unilateral contact and friction boundary value problem in the theory of elasticity", Ing.Arch. 44, 421-432, 1975.

[36] CAMPOS, L.T.; ODEN, J.T.; KIKUCHI, N.; "A numerical analysis of a class of contact problems with friction in elastostatics", Comp. Meht.Appl.Mech.Engng., 34, 821-845, 1982.

[37] RAOUS, M.; "Contacts unilateraux avec frottement en viscoélasticité", in *Unilateral Problems in Structural Analysis*, Ed. G. Del Piero, F. Maceri, CISM Courses and Lectures, Springer, 1985.
[38] RAOUS, M.; "On two variational inequalities arising from a periodic viscoelastic unilateral problem", in *Variational Inequalities an Complementarity Problems*, Ed. Gianessi-Cottle-Lions, J. Wiley, 1979.

[39] RAOUS, M.; "Fissuration sous constraints alternées en viscoélasticité et viscoplasticité", *These*, Univ. de Provence, Marseille, 1980.

[40] COHN, M.Z.; MAIER, G.; (Editors) "Engineering Plasticity by Mathematical Programming", *Proceedings of the NATO Advanced Study Institute*, Univ. of Waterloo, Waterloo, Ontario, Canada, Aug, 1977, pub. by Pergamon Press, 1979.

[41] HAUG, E.J.; CEA, J.; "Optimization of Distributed Parameter Structures", *Proceedings of the NATO Advanced Study Institute on Optimization of Distributed Parameter Structural Systems*, Iowa, may-june, 1980, pub. by Sijthoff-Noordhoff.