# Boundary minimum principles for the unilateral contact problems

P P Panagiotoopoulos, P P Lazaridis

► To cite this version:

P P Panagiotoopoulos, P P Lazaridis. Boundary minimum principles for the unilateral contact problems. International Journal of Solids and Structures, 1987, 23 (11), pp.1465-1484. 10.1016/0020-7683(87)90064-3. hal-01406291

HAL Id: hal-01406291

https://hal.science/hal-01406291

Submitted on 1 Dec 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# BOUNDARY MINIMUM PRINCIPLES FOR THE UNILATERAL CONTACT PROBLEMS

P. D. PANAGIOTOPOULOS

Department of Civil Engineering, Aristotle University, Thessaloniki, Greece; Institute for Technical Mechanics, R.W.T.H. Aachen, Federal Republic of Germany

and

P. P. LAZARIDIS

Department of Civil Engineering, Aristotle University, Thessaloniki, Greece

**Abstract**—The aim of the present paper is the derivation of boundary variational principles for the unilateral contact problem. Using the inequality constrained principles of minimum potential and complementary energy and the equivalent variational inequality formulations we derive first saddle point formulations for the problems using appropriate Langrangian functions. An elimination technique allows the formulation of two minimum “principles” on the boundary with respect to the unknown normal displacements and reactions of the contact region, respectively. It is also shown that these two minimum problems are equivalent to multivalued boundary integral equations involving symmetric operators. The theory is illustrated by numerical examples solved both by the FEM and the BEM.

## 1. INTRODUCTION

In recent years a large number of structural problems involving unilateral constraints have been studied [1, 2]. The convexity of the arising superpotentials [3] leads for a large class of problems to variational inequalities expressing the principle of virtual or complementary virtual work in its inequality form. These variational inequalities are equivalent to the minimum of potential or complementary energy which after discretization permits the numerical treatment by using an appropriate non-linear programming algorithm. One important class of unilateral problems is the unilateral contact problems arising when an elastic body is in ambiguous contact with a rigid or deformable support. The term “ambiguous” means that we do not know a priori which parts of the body are in contact with the support and which not and this fact renders the problem unilateral. We pay special attention to the case of a rigid support since the unilateral contact problem with a deformable support can be reduced to the rigid support problem. This is the problem posed by Signorini in 1933 and studied by Fichera in 1963 [4, 5].

For the numerical treatment we have to solve an inequality constrained quadratic programming (QP) problem either with respect to the displacements (minimum of potential energy) or with respect to the stresses (minimum of complementary energy). Several techniques have been applied solving directly or indirectly a minimization problem (see e.g. Refs. [6, 7] and the references cited in Ref. [2]). However, all these techniques have as a main disadvantage that due to the use of QP algorithms, either the size of the structure to be solved is small (substructuring is necessary), or major changes of the existing general FE programs are necessary and considerable computer time is needed. In order to avoid all these cases an active constraint strategy method was proposed in Refs. [8, 9]. This method is based on a “translation” in the language of structural analysis, of the QP theorems of Theil and Van de Panne and uses completely the existing classical linear structural analysis computer programs. Thus unilateral problems with a large number of unknowns could be safely solved. However, the main disadvantage of this method is, generally speaking, the low level of automatization concerning the determination of each linear substructure, i.e. the determination of the active constraints.
In order to avoid this difficulty we develop here a more delicate method, based on Lagrangian formulations, which leads to new inequality constrained minimization problems with respect to the unknown reactions (or displacements) of the contact zone. Thus we obtain discrete QP problems having symmetric full matrices but with a drastically reduced number of unknowns in comparison to the initially formulated minimum problem of the potential or the complementary energy. It is worth noting that this method reveals many of the properties of the classical boundary integral equation methods and from this standpoint the proposed method can be seen as a first attempt towards the effective calculation of inequality problems by boundary integral "equations". The developed method in the present paper holds for a large class of inequality (or unilateral) problems (i.e. for problems leading to variational inequalities) besides the unilateral contact problem: we mean here all the unilateral problems for which the minimum of the potential or the complementary energy is constrained by linear inequality constraints and does not include non-differentiable superpotentials (cf. in the friction problem the complementary energy formulation). The fact that non-linear constraints can be linearized as well to yield the same situation as before extends considerably the validity of the proposed method.

## 2. PRIMAL, DUAL AND MIXED FORMULATIONS OF THE PROBLEM

We consider a three-dimensional linear elastic body. The procedure followed is general and holds also for shells, plates, beams, etc. generally speaking for all structures permitting a formulation of the equilibrium problem in terms of a Lagrangian formulation.

Let $\Omega$ be a subset of the three-dimensional Euclidean space $\mathbb{R}^3$ with a boundary $\Gamma$. $\Omega$ is occupied by a linear elastic body and is referred to an orthogonal Cartesian coordinate system $0x_1x_2x_3$. $\Gamma$ is decomposed into three mutually disjoint parts $\Gamma_U$, $\Gamma_F$ and $\Gamma_S$. We demand that on $\Gamma_U$ (resp. $\Gamma_F$) the displacements (resp. the tractions) are given and that on $\Gamma_S$ the Signorini-Fichera boundary conditions hold. We denote also by $n = \{n_i\}$ the outward unit normal vector to $\Gamma$ by $S = \{S_i\} = \{\sigma_{ij}n_j\}$ the traction vector on the boundary, where $\sigma = \{\sigma_{ij}\}$ is the stress vector. Further $S_N$ (resp. $S_T$) is the normal (resp. the tangential) component of $S$ with respect to $\Gamma$; $u_N$ and $u_T$ are the corresponding components of the displacement $u$. Let $\varepsilon = \{\varepsilon_{ij}\}$ be the strain tensor (assumption of small strains) and $C = \{C_{ijk}\} i,j,h,k = 1,2,3$ Hooke's tensor of elasticity obeying the well-known symmetry

$$C_{ijk} = C_{jik} = C_{khi} \qquad (1)$$

and ellipticity conditions

$$C_{ijk} \cdot \varepsilon_{ij} \cdot \varepsilon_{hk} \ge c \cdot \varepsilon_{ij} \cdot \varepsilon_{hk}, \quad \forall \varepsilon = \{\varepsilon_{ij}\} \in \mathbb{R}^3, \ c \ \text{const.} > 0. \qquad (2)$$

On $\Gamma_U$ we assume that

$$u_i = U_i, \quad U_i = U_i(x) \qquad (3)$$

and on $\Gamma_F$

$$S_i = F_i, \quad F_i = F_i(x). \qquad (4)$$

The Signorini-Fichera boundary condition reads: if $u_N < 0$ then $S_N = 0$ (no contact), if $u_N = 0$ then $S_N \le 0$ (contact), or equivalently (see Fig. 1(a))

<figure><img src="image_2.png" /><figcaption>Fig. 1. Reaction-displacement diagrams for the unilateral contact laws.</figcaption></figure>
$$S_N \le 0, \quad u_N \le 0, \quad u_N S_N = 0 \quad \text{on } \Gamma_S. \qquad (5)$$

Relations (5) must be completed with the condition in the tangential direction

$$S_{T_i} = C_{T_i}, \quad C_{T_i} = C_{T_i}(x) \quad \text{on } \Gamma_S \qquad (6)$$

where $C_T$ is a prescribed tangential force distribution. Besides law (5) we can consider the more general case depicted in Fig. 1(a) (dotted lines) which corresponds to the unilateral contact with a rigid support at distance $h$ from the body.

It is interesting to note that the more realistic unilateral contact laws with a linearly deformable support of Fig. 1(b) (Winkler spring) lead to the same variational expressions as the ones of Fig. 1(a). Indeed, it is sufficient to enlarge the body $\Omega$ by fictitious linear-elastic springs of zero length along $\Gamma_S$ with spring constant $k$ and assuming that the laws of Fig. 1(a) hold. The foregoing holds also for any physically meaningful combination of the diagrams of Fig. 1(a) with the linear $S_N - u_N$ law.

Further we have

$$\sigma_{ij,j} + f_i = 0 \quad \text{in } \Omega \qquad (7)$$

$$\varepsilon_{ij} = \varepsilon_{ij}(u) = \frac{1}{2}(u_{i,j} + u_{j,i}) \quad \text{in } \Omega \qquad (8)$$

$$\sigma_{ij} = C_{ijkk} \cdot \varepsilon_{hk} \quad \text{in } \Omega \qquad (9)$$

where $f = \{f_i\}$ represents the volume force vector and the comma denotes partial derivation. We denote further by $\tilde{V}$ the vector space of the displacement $v_i$ and let $V$ be the set of kinematically admissible displacement fields

$$V = \{v | v = \{v_i\}, \ v_i \in \tilde{V}, \ v_i = U_i \quad i = 1, 2, 3 \text{ on } \Gamma_V\} \qquad (10)$$

without taking into account the constraints on $\Gamma_S$. Due to the Signorini-Fichera boundary conditions the kinematically admissible set of displacements takes the form

$$K = \{v | v = \{v_i\}, \ v_i \in V, \ V_N \le 0 \quad \text{on } \Gamma_S\}. \qquad (11)$$

We denote by $(f, v)$ the work done by the force $f = \{f_i\}$ for the displacement $v = \{v_i\}$ on $\Omega$ and by $[f, v]_\Gamma$ the corresponding work on $\tilde{\Gamma} \subset \Gamma$ (i.e. $\int_\Omega f_i v_i d\Omega$), etc. Further let

$$a(u, v) = (C\varepsilon(u), \varepsilon(v)) = \int_{\Omega} C_{ijkk} \cdot \varepsilon_{ij}(u) \cdot \varepsilon_{hk}(v) \, d\Omega \qquad (12)$$

be the bilinear form of elasticity and let $\Pi$ be the potential energy

$$\Pi(v) = \frac{1}{2}a(v, v) - (f, v) - [C_T, v_T]_{\Gamma_S} - [F, v]_{\Gamma_F}. \qquad (13)$$

It is well known[4, 5] that the problem

$$\Pi(u) = \min \{\Pi(v) | v \in K\} \qquad (14)$$

characterizes the position of equilibrium. Problem (14) has one and only one solution (for $v_i \in H^1(\Omega)$—the Sobolev space

$$C_{ijkk} \in L^\infty \quad \text{and} \quad f_i, C_{T_i}, F_i \in L^2, U_i \in H^{1/2}$$

which equivalently, satisfies the variational inequality
$$u \in K, \quad a(u, v-u) - l(v-u) \ge 0 \quad \forall v \in K \quad (15)$$

with $l(v) = (f, v) + [F, v]_{\Gamma_r} + [C_T, v_T]_{\Gamma_s}$. Relations (14) or (15) are the primal formulations of the BVP.

For the mixed formulation we introduce the convex subset (which is closed in the previously mentioned functional framework)

$$L = \{ \mu_N | \mu_N \le 0 \text{ on } \Gamma_S \}. \quad (16)$$

Let $u_0$ be a kinematically admissible displacement field, i.e. such that $u_{0i} = U_i$ on $\Gamma_u$, and let

$$\bar{u} = u - u_0, \quad \bar{v} = v - u_0 \quad (17)$$

where

$$\bar{u}, \bar{v} \in V_0 = \{v | v = \{v_i\}, \ v_i \in \bar{V}, \ v_i = 0 \text{ on } \Gamma_v\}. \quad (18)$$

Then relations (15) take the form: find $u = \bar{u} + u_0 \in K$ such that

$$a(\bar{u}, \bar{v} - \bar{u}) - l(\bar{v} - \bar{u}) + a(u_0, \bar{v} - \bar{u}) \ge 0 \quad \forall \bar{v} = \bar{v} + u_0 \in K \quad (19)$$

and relation (14) becomes

$$\tilde{\Pi}(\bar{u}) = \min \{ \tilde{\Pi}(\bar{v}) | \bar{v} \in \tilde{K} \} \quad (20)$$

where

$$\tilde{\Pi}(\bar{v}) = \frac{1}{2} a(\bar{v}, \bar{v}) - l(\bar{v}) + a(u_0, \bar{v}) = \Pi(v) - \Pi(u_0) \quad (21)$$

and

$$\tilde{K} = \{ \bar{v} | \bar{v}_N + u_{0N} \le 0 \text{ on } \Gamma_S \}. \quad (22)$$

Note that more generally $u_0$ may represent any other initial strain field (temperature distributions, given dislocations, etc.). We denote by $\tilde{I}$ the functional

$$\tilde{I}(\bar{v}) = l(\bar{v}) - a(u_0, \bar{v}). \quad (23)$$

Here $\mu_N$ is the Lagrange multiplier for the problem. Through this Lagrange multiplier we introduce the boundary condition (5) on $\Gamma_S$. Indeed as it is easy to verify

$$\inf_{\mathcal{R}} \tilde{\Pi}(\bar{v}) = \inf_{v_0} (\tilde{\Pi}(\bar{v}) + I_R(\bar{v}))$$

where

$$I_R(\bar{v}) = \begin{cases} 0 & \text{if } \bar{v}_N + u_{0N} \le 0 \text{ on } \Gamma_S \\ \infty & \text{otherwise.} \end{cases} \quad (24)$$

But

$$I_R(\bar{v}) = \sup_{\mu_N \le 0} [[-\mu_N, \bar{v}_N + u_{0N}]_{\Gamma_S}] \quad (25)$$

and thus
$$ \inf_{\mathcal{R}} \bar{\Pi}(\bar{v}) = \inf_{\bar{v} \in V_0, \mu_N \in L} \sup \{\bar{\Pi}(\bar{v}) - [\mu_N, \bar{v}_N + u_{0N}]_{\Gamma_S}\}. \quad (26) $$

The Lagrangian of the problem is a real-valued function $\mathcal{L}$ on $V_0 \times L$ defined by the relation

$$ \mathcal{L}(\bar{v}, \mu_N) = \frac{1}{2}a(\bar{v}, \bar{v}) - [\mu_N, \bar{v}_N + u_{0N}]_{\Gamma_S} - [C_T, \bar{v}_T]_{\Gamma_S} - [F, \bar{v}]_{\Gamma_F} - (f, \bar{v}) + a(u_0, \bar{v}). \quad (27) $$

Thus the mixed variational formulation of the Signorini-Fichera problem now reads (cf. Refs [10-12]): find the saddle point $\{\bar{w}, \lambda_N\} \in V_0 \times L$ of $\mathcal{L}$ on $V_0 \times L$, i.e.

$$ \mathcal{L}(\bar{w}, \mu_N) \le \mathcal{L}(\bar{w}, \lambda_N) \le \mathcal{L}(\bar{v}, \lambda_N) \quad \forall \bar{v} \in V_0, \mu_N \in L. \quad (28) $$

It can be proved in the previously mentioned functional framework by the methods given on p. 57 of Ref. [13] that problem (28) has a unique solution $\{\bar{w}, \lambda_N\} \in V_0 \times L$ such that $\bar{w} = \bar{u} \in \bar{K}$ and $\lambda_N = S_N(\bar{u})$ on $\Gamma_S$.

The dual formulation with respect to the stresses results easily from problem (28) if we consider the statically admissible set of stresses

$$ \Lambda = \{\sigma | \sigma = \{\sigma_{ij}\}, \sigma_{ij} = \sigma_{ji}, \sigma_{ij,j} + f_i = 0 \text{ in } \Omega, S_N \in L, S_{T_i} = C_{T_i} \text{ on } \Gamma_S, S_i = F_i \text{ on } \Gamma_F\}. \quad (29) $$

The stress field at the position of equilibrium $\sigma$ is characterized by the minimum of the complementary energy[14]

$$ \Pi^c(\sigma) = \frac{1}{2} A(\sigma, \sigma) - [U, \sigma]_{\Gamma_U} \quad \text{over } \Lambda, \text{ i.e.} \quad (30) $$

$$ \Pi^c(\sigma) = \min \{\Pi^c(\tau) | \tau \in \Lambda\}. \quad (31) $$

Here $\frac{1}{2}A(\sigma, \sigma) = (c\sigma, \sigma)$ where $c = \{c_{ijk}\}$ is the inverse tensor to $C$.

### 3. FORMULATION WITH RESPECT TO THE TRACTIONS OF THE CONTACT AREA

If $\{\bar{w}, \lambda_N\} \in V_0 \times L$ is the solution of problem (28) then we can write problem (28) as it is obvious equivalently in the form

$$ \mathcal{L}(\bar{w}, \lambda_N) = \inf_{v_0} \sup_L \mathcal{L}(\bar{v}, \mu_N) = \sup_L \inf_{v_0} \mathcal{L}(\bar{v}, \mu_N). \quad (32) $$

We set

$$ \inf_{v_0} \mathcal{L}(\bar{v}, \mu_N) = \bar{\Pi}_1(\mu_N) \quad (33) $$

assuming for the moment that $\mu_N \in L$ is given. Then problem (33) is equivalent to the following bilateral problem: find $\tilde{u} = \tilde{u}(\mu_N) \in V_0$ such that

$$ a(\tilde{u}, \tilde{v}) - [\mu_N, \tilde{v}_N]_{\Gamma_S} - (f, \tilde{v}) - [F, \tilde{v}]_{\Gamma_F} - [C_T, \tilde{v}_T]_{\Gamma_S} + a(u_0, \tilde{v}) = 0 \quad \forall \tilde{v} \in V_0. \quad (34) $$

Obviously problem (34) is the expression of the principle of virtual work for a fictive structure resulting from the initial unilateral one by eliminating the unilateral constraints on $\Gamma_S$ and by adding the corresponding reactions $\mu_N$. The position of equilibrium of this structure is characterized by the minimization problem (33) for the potential energy of this fictive structure. The solution $\tilde{u}$ of problem (34) can be considered, due to the linearity of problem (34), as the sum of $\tilde{u}_1 \in V_0$ and $\tilde{u}_2 \in V_0$ which are solutions of the two following bilateral problems:
<figure><img src="image_3.png" /><figcaption>Fig. 2. The problem decomposition (force method).</figcaption></figure>

$$a(\tilde{u}_1, \tilde{v}) - (f, \tilde{v}) - [F, \tilde{v}]_{\Gamma_F} - [C_T, \tilde{v}_T]_{\Gamma_S} + a(u_0, \tilde{v}) = 0 \quad \forall \tilde{v} \in V_0 \quad (35)$$

$$a(\tilde{u}_2, \tilde{v}) - [\mu_N, \tilde{v}_N]_{\Gamma_S} = 0 \quad \forall \tilde{v} \in V_0 \quad (36)$$

respectively. Obviously both problems describe the equilibrium configuration of two bilateral structures resulting from the initial one by ignoring the unilateral support and assuming the appropriate boundary parts with zero loading.

In the case of problem (35) the structure is loaded by the forces $f$ in $\Omega$ and $C_T$ on $\Gamma_S$ tangentially and $F$ on $\Gamma_F$, whereas on $\Gamma_S$ the normal loading is zero. Moreover, the initial displacement field $u_0$ is considered. In the case of problem (36) the structure is subjected to normal forces $\mu_N$ on $\Gamma_S$ and we assume zero forces in $\Omega$ on $\Gamma_F$ and on $\Gamma_S$ in the tangential direction. Accordingly the solutions $\tilde{u}_1$ and $\tilde{u}_2$ are uniquely determined, as it is well known from the classical (bilateral) elasticity (cf. e.g. Ref. [12]). For these two classical bilateral structures we can write the solution in terms of Green's operator $G$. This operator is the same for the two structures due to the "same" boundary conditions holding in both cases in the form of Fig. 2. Thus we may write in both cases the solution of the problem as follows:

$$\tilde{u}_1 = G(\tilde{I}), \quad \tilde{u}_2 = G(\mu_N), \quad \tilde{u} = \tilde{u}_1 + \tilde{u}_2, \quad \tilde{e} = \{f, F, C_T, u_0\}. \quad (37)$$

Note that the yet unknown force distribution $\mu_N$ must be admissible in the sense of relation (16), i.e.

$$\mu_N \in L. \quad (38)$$

Thus from problems (33), (35) and (36) we obtain by setting $\tilde{v} = \tilde{u}_1$ in problem (35) and $\tilde{v} = \tilde{u}_2$ in problem (36) that

$$\begin{aligned} \tilde{f}_1(\mu_N) &= \mathcal{L}(\tilde{u}, \mu_N) = \frac{1}{2}a(\tilde{u}, \tilde{u}) - [\mu_N, \tilde{u}_N + u_{0N}]_{\Gamma_S} - [C_T, \tilde{u}_T]_{\Gamma_S} - (f, \tilde{u}) - [F, \tilde{u}]_{\Gamma_F} + a(u_0, \tilde{u}) \\ &= -[\mu_N, \tilde{u}_N]_{\Gamma_S} - \frac{1}{2}[\mu_N, \tilde{u}_N]_{\Gamma_S} - \frac{1}{2}(f, \tilde{u}_1) - \frac{1}{2}[F, \tilde{u}_1]_{\Gamma_F} \\ &\quad - \frac{1}{2}[C_T, \tilde{u}_T]_{\Gamma_S} - [\mu_N, u_{0N}]_{\Gamma_S} + \frac{1}{2}a(u_0, \tilde{u}_1) \\ &= -[\mu_N, [G(\tilde{I})]_N]_{\Gamma_S} - \frac{1}{2}[\mu_N, [G(\mu_N)]_N]_{\Gamma_S} + \frac{1}{2}a(u_0, G(\tilde{I})) \\ &\quad - \frac{1}{2}(f, G(\tilde{I})) - \frac{1}{2}[F, G(\tilde{I})]_{\Gamma_F} - \frac{1}{2}[C_T, G(\tilde{I})]_{\Gamma_S} - [\mu_N, u_{0N}]_{\Gamma_S}. \end{aligned} \quad (39)$$

Further we denote by $\beta$ the bilinear form

$$\beta(\mu_N, v_N) = [\mu_N, [G(v_N)]_N]_{\Gamma_S} \quad (40)$$

and by $\gamma$ the linear form
$$ \gamma(\mu_N) = -[\mu_N, [G(\mathcal{I})]_N]_{\Gamma_S} - [\mu_N, u_{0N}]_{\Gamma_S}. \quad (41) $$

Thus

$$ \begin{aligned} \tilde{\Pi}_1(\mu_N) = & -\frac{1}{2}\beta(\mu_N, \mu_N) + \gamma(\mu_N) - \frac{1}{2}(f, G(\mathcal{I})) - \frac{1}{2}[F, G(\mathcal{I})]_{\Gamma_F} \\ & - \frac{1}{2}[C_T, G(\mathcal{I})]_{\Gamma_S} + \frac{1}{2}a(u_0, G(\mathcal{I})). \end{aligned} \quad (42) $$

From eqns (32), (33) and (42) we can formulate the following minimization problem with respect to the unknown boundary tractions $\mu_N$:

$$ \min \{ \Pi_1(\mu_N) = \frac{1}{2}\beta(\mu_N, \mu_N) - \gamma(\mu_N) \mid \mu_N \in L \}. \quad (43) $$

The term $-\frac{1}{2}(f, G(\mathcal{I})) - \frac{1}{2}[F, G(\mathcal{I})]_{\Gamma_F} - \frac{1}{2}[C_T, G(\mathcal{I})]_{\Gamma_S} + \frac{1}{2}a(u_0, G(\mathcal{I}))$ does not depend on $\mu_N$ and for this reason is omitted. We denote obviously again by $\lambda_N$ the solution of this problem. From eqns (33), (27) and (34) we obtain another expression for $\tilde{\Pi}_1$

$$ \begin{aligned} \tilde{\Pi}_1(\mu_N) &= \inf_{v_0} \mathcal{L}(\bar{v}, \mu_N) = \mathcal{L}(\tilde{u}, \mu_N) \\ &= \frac{1}{2}a(\tilde{u}, \tilde{u}) - [\mu_N, \tilde{u}_N + u_{0N}]_{\Gamma_S} - [F, \tilde{u}]_{\Gamma_F} - (f, \tilde{u}) - [C_T, \tilde{u}_T]_{\Gamma_S} + a(u_0, \tilde{u}) \\ &= \frac{1}{2}a(\tilde{u}, \tilde{u}) - a(\tilde{u}, \tilde{u}) - [\mu_N, u_{0N}]_{\Gamma_S} = -\frac{1}{2}a(\tilde{u}, \tilde{u}) - [\mu_N, u_{0N}]_{\Gamma_S}. \end{aligned} \quad (44) $$

Thus for $\mu_N = \lambda_N$ we have obviously $\tilde{u} = \bar{w}$ and thus

$$ \tilde{\Pi}_1(\lambda_N) = -[\frac{1}{2}a(\bar{w}, \bar{w}) + [\lambda_N, u_{0N}]_{\Gamma_S}]. \quad (45) $$

Further we shall show that

$$ \lambda_N = S_N(\tilde{u}) \quad \text{on} \quad \Gamma_S \quad (46) $$

where $\tilde{u} \in \tilde{K}$ is the solution of the primal problem (20). Indeed let us choose $\lambda_N \in L$ such as to satisfy eqn (46) where $\tilde{u}$ is a solution of eqn (20) and $\bar{w} = \tilde{u} \in \tilde{K}$. Then $(\tilde{u}, \lambda_N)$ is the unique solution of the saddle point problem (28).

The proof will be completed if we prove, first, that the chosen $\lambda_N$ is a minimizer of $\tilde{\Pi}_1$ over $L$ and secondly that the solution of eqn (43) is unique. Indeed from eqns (28) and (33) we obtain that

$$ \mathcal{L}(\bar{w}, \lambda_N) = \inf_{\bar{v} \in V_0} \mathcal{L}(\bar{v}, \lambda_N) = \tilde{\Pi}_1(\lambda_N). \quad (47) $$

But

$$ \mathcal{L}(\bar{w}, \lambda_N) = \sup_L \inf_{v_0} \mathcal{L}(\bar{v}, \mu_N) = \sup_L \tilde{\Pi}_1(\mu_N) \quad (48) $$

and thus

$$ \tilde{\Pi}_1(\lambda) = \inf_L \tilde{\Pi}_1(\mu). \quad (49) $$

It remains to show the uniqueness of the solution of problem (43). Indeed in the functional framework introduced before, it can be shown also that problem (43) has a unique solution. The proof uses functional analytic theories and is omitted here (see in this context Ref. [15]).
<figure><img src="image_4.png" /><figcaption>Fig. 3. The Betti theorem.</figcaption></figure>

It is interesting to note the fact that the quadratic form $\beta(\mu_N, \mu_N)$ is symmetric. Indeed from eqn (40) using the linearity of $G$ and the reciprocal theorem of Betti (p. 391 in Ref. [16]) we can easily verify that (see also Fig. 3)

$$ \beta(\mu_N, v_N) = \beta(v_N, \mu_N). \qquad (50) $$

We close this section by giving two equivalent formulations of problem (43). The first is a variational inequality: find $\lambda_N \in L$ such as to satisfy

$$ \beta(\lambda_N, \mu_N - \lambda_N) - \gamma(\mu_N - \lambda_N) \ge 0 \quad \forall \mu_N \in L \qquad (51) $$

and the second a multivalued integral equation on the boundary part $\Gamma_S$ of the structure which reads

$$ \gamma - \frac{1}{2} \operatorname{grad} \beta(\lambda_N, \lambda_N) \in \partial I_L(\lambda_N) \qquad (52) $$

where $I_L(\lambda_N)$ is the indicator of the admissible set $L$, i.e.

$$ I_L(\lambda_N) = \begin{cases} 0 & \text{if } \lambda_N \le 0 \quad (\text{i.e. } \lambda_N \in L) \\ \infty & \text{otherwise } (\text{i.e. } \lambda_N \notin L). \end{cases} \qquad (53) $$

Relation (51) is equivalent to problem (43) by using the well-known results in the theory of variational inequalities (see Refs [1, 2, 13] and p. 40 in Ref. [4]). Relation (52) is equivalent to inequality (51) by the definition of the subdifferential $\partial$ (see e.g. Ref. [13]). Finally we pay some attention to the other boundary conditions depicted in Fig. 1. They give rise only with slight modifications to the same expressions as the classical Signorini-Fichera boundary conditions studied, which correspond to Fig. 1(a)(1). Thus for Fig. 1(a)(2) $K$ in relation (11) is defined by the inequality $u_N \ge h$. For boundary conditions of Fig. 1(b) the only difference is that we have to consider the energy of the added fictive springs with constant $k$. This amounts to say that Green's operator in eqns (37) must be taken for the enlarged structure with the attached fictitious linear springs.

### 4. FORMULATION WITH RESPECT TO THE DISPLACEMENTS OF THE CONTACT AREA

Let us consider again the saddle-point formulation (32). In the previous section, in order to derive a variational formulation on the boundary with respect to the unknown unilateral contact tractions, we have considered $\sup_{L} \inf_{v_0} \mathcal{L}(\bar{v}, \mu_N)$ with $\mu_N$ as given.

Here we will work on $\inf_{v_0} \sup_L \mathcal{L}(\bar{v}, \mu_N)$ or on an equivalent saddle-point formulation by considering $\bar{v}$ as given. We note that from eqns (17)
$$
\begin{align*}
\inf_{v_0} \sup_L \mathcal{L}(\bar{v}, \mu_N) &= \sup_L \inf_{v_0} \mathcal{L}(\bar{v}, \mu_N) = \sup_L \inf_{v_0} [\frac{1}{2}a(\bar{v}, \bar{v}) - [\mu_N, \bar{v}_N + u_{0N}]_{\Gamma_S} \\
&\quad - [C_T, \bar{v}_T]_{\Gamma_S} - [F, \bar{v}]_{\Gamma_F} - (f, \bar{v}) + a(u_0, \bar{v})] \\
&= \sup_{\mu_N \in L} \inf_{v \in V} [\frac{1}{2}a(v - u_0, v - u_0) - [\mu_N, v_N]_{\Gamma_S} - [C_T, v_T - u_{0T}]_{\Gamma_S} \\
&\quad - [F, v - u_0]_{\Gamma_F} - (f, v - u_0) + a(u_0, v - u_0)] \\
&= \sup_{\mu_N \in L} \inf_{v \in V} [\frac{1}{2}a(v, v) + \frac{1}{2}a(u_0, u_0) - a(v, u_0) - [\mu_N, v_N]_{\Gamma_S} \\
&\quad - [C_T, v_T - u_{0T}]_{\Gamma_S} - [F, v - u_0]_{\Gamma_F} - (f, v - u_0) + a(u_0, v) - \frac{1}{2}a(u_0, u_0)] \\
&= \inf_{v \in V} \sup_{\mu_N \in L} [\frac{1}{2}a(v, v) - [\mu_N, v_N]_{\Gamma_S} - [C_T, v_T]_{\Gamma_S} - [F, v]_{\Gamma_F} \\
&\quad - (f, v) - \Pi(u_0)] = \inf_{v \in V} \sup_{\mu_N \in L} \mathcal{L}(v, \mu_N) - \Pi(u_0). \tag{54}
\end{align*}
$$

Here we have introduced the modified Lagrangian

$$
\mathcal{L}(v, \mu_N) = \frac{1}{2}a(v, v) - [\mu_N, v_N]_{\Gamma_S} - [C_T, v_T]_{\Gamma_S} - [F, v]_{\Gamma_F} - (f, v). \quad (55)
$$

We recall that

$$
\sup_{\mu_N \le 0} -[\mu_N, v_N]_{\Gamma_S} = \begin{cases} 0 & \text{if } v_N \le 0 \\ \infty & \text{if } v_N > 0 \end{cases} = I_k(v) \quad (56)
$$

and write

$$
\inf_{v \in V} \sup_{\mu_N \in L} \mathcal{L}(v, \mu_N) = \inf_{v \in V} \left[ \frac{1}{2} a(v, v) + I_k(v) - l(v) \right] = \inf_{v \in K} \Pi(v). \quad (57)
$$

But

$$
\frac{1}{2}a(v,v) = \frac{1}{2}(C \cdot \varepsilon(v) \cdot \varepsilon(v)) = \sup_{\tau \in \Sigma} ((\tau, \varepsilon(v)) - \frac{1}{2}A(\tau, \tau)) \quad (58)
$$

where $\Sigma$ denotes the set of all symmetric stress tensors. Indeed by simple derivation with
respect to $\tau$ we find that the supremum is attained for

$$
\epsilon_{ij}(v) = c_{ijhk} \cdot \sigma_{hk} \quad \text{or} \quad \sigma_{ij} = C_{ijhk} \cdot \epsilon_{hk}(v) \tag{59}
$$

which yields eqns (58). Thus

$$
\inf_{v} \sup_{L} \mathcal{L}(v, \mu_N) = \inf_{v \in K} \Pi(v) = \inf_{v \in K} \sup_{\tau \in \Sigma} \left\{ (\tau, \varepsilon(v)) - (f, v) - [F, v]_{\Gamma_F} - [C_T, v_T]_{\Gamma_S} - \frac{1}{2} A(\tau, \tau) \right\}. \quad (60)
$$

We note also that for every $v \in V$ the Green–Gauss theorem implies that

$$
(\tau, \varepsilon(v)) = -(\tau_{ij,j}, v_i) + [T_N, v_N]_{\Gamma_S} + [T_T, v_T]_{\Gamma_S} + [T, v]_{\Gamma_F} + [T, U]_{\Gamma_U} \quad (61)
$$

where $T, T_N$ and $T_T$ are the boundary tractions corresponding to $\tau$. Note that eqn (61) is the principle of virtual work for the free structure. From eqns (60) and (61) we obtain that
$$
\begin{align}
& \inf_{v \in K} \sup_{\tau \in \Sigma} \left[ -(\tau_{ij,j} + f_i, v_i) + [T-F, v]_{\Gamma_f} + [T, U]_{\Gamma_u} + [T_N, v_N]_{\Gamma_S} + [T_T - C_T, v_T]_{\Gamma_S} - \frac{1}{2} A(\tau, \tau) \right] \nonumber \\
& \qquad = \inf_{v \in K} \sup_{\tau \in \Sigma_i} \left[ [T, U]_{\Gamma_v} + [T_N, v_N]_{\Gamma_S} - \frac{1}{2} A(\tau, \tau) \right] \tag{62}
\end{align}
$$

where

$$
\Sigma_1 = \{\tau | \tau = \{\tau_{ij}\}, \quad \tau_{ij} = \tau_{ji}, \quad \tau_{ij,j} + f_i = 0 \text{ in } \Omega, \quad T_i = f_i \text{ on } \Gamma_f\}
$$

$$
T_{T_i} = C_{T_i} \text{ on } \Gamma_S. \quad (63)
$$

Indeed

$$
\begin{align}
& \sup_{\tau \in \Sigma} \left[ -(\tau_{ij,j} + f_i, v_i) + [T-F, v]_{\Gamma_f} + [T_T - C_{T,v}]_{\Gamma_S} \right] \nonumber \\
& \qquad = \begin{cases} 0 & \text{if } \tau_{ij,j} + f_i = 0, \quad T_i = F_i, \quad T_{T_i} = C_{T_i} \\ \infty & \text{otherwise.} \end{cases} \tag{64}
\end{align}
$$

Thus

$$
\inf_{v \in K} \sup_L \mathcal{P}(v, \mu_N) = \inf_{v \in K} \Pi(v) = \inf_{v \in K} \sup_{\tau \in \Sigma_1} L(v_N, \tau) \quad (65)
$$

where the new Lagrangian $L$ is given by

$$
L(v_N, \tau) = [T, U]_{\Gamma_v} + [T_N, v_N]_{\Gamma_S} - \frac{1}{2} A(\tau, \tau). \quad (66)
$$

Note also that due to the duality[2, 13] between the primal and the dual problem (14) and
(30) the relation

$$
\begin{align}
\inf_{v \in K} \Pi(v) &= \sup_{\tau \in \Lambda} -\Pi^c(\tau) = \sup_{\tau \in \Sigma_1} \left(-\frac{1}{2} A(\tau, \tau) + [T, U]_{\Gamma_v} - I_L(\tau)\right) \nonumber \\
&= \sup_{\tau \in \Sigma_1} \inf_{v \in K} \left(-\frac{1}{2} A(\tau, \tau) + [T, U]_{\Gamma_v} + [T_N, v_N]_{\Gamma_S}\right) = \sup_{\tau \in \Sigma_1} \inf_{v \in K} (v_N, \tau) \tag{67}
\end{align}
$$

holds. Here $\Lambda$ is defined in eqns (29)

$$
I_L(T) = \begin{cases} 0 & \text{if } T_N \le 0 \quad (\text{i.e. if } T_N \in L) \\ \infty & \text{otherwise} \quad (\text{i.c. if } T_N \notin L) \end{cases} \tag{68}
$$

and

$$
\inf_{v_N \le 0} [T_N, v_N]_{\Gamma_S} = -I_L(T). \tag{69}
$$

Thus

$$
\inf_{v \in K} \sup_{\tau \in \Sigma_1} L(v_N, \tau) = \sup_{\tau \in \Sigma_1} \inf_{v \in K} L(v_N, \tau) = L(u_N, \sigma). \quad (70)
$$

The last equality in (70) can be easily verified by writing the foregoing relations for the
actual solution {u, σ, S}, etc. (i.e. eqn (65), etc. without inf sup, inf, etc.). A useful remark
is that besides L and another Lagrangian has also been obtained, namely the expression
on the right-hand side of eqn (60). This Lagrangian can be obtained directly if one applies
the saddle-point theory of Ref. [13].
From the above relations and the fact that the saddle-point problem (28) has a unique solution we can easily verify that the corresponding saddle-point problem with respect to $L$ (cf. eqn (60)) admits also a unique solution $\sigma \in \Sigma_1$ and $u_N = u_N(\sigma) \le 0$ (cf. also Ref. [13]).
Let us denote further

$$ \sup_{\tau \in \Sigma_1} L(v_N, \tau) = \tilde{\Pi}_2(v_N) \quad \text{or} \quad \inf_{\tau \in \Sigma_1} -L(v_N, \tau) = -\tilde{\Pi}_2(v_N) \quad (71) $$

assuming that $v_N \le 0$ is prescribed. Then eqns (71) correspond to the variational inequality:
find $\tilde{\sigma} = \tilde{\sigma}(v) \in \Sigma_1$ such that

$$ A(\tilde{\sigma}, \tau - \tilde{\sigma}) - [U, T - \tilde{S}]_{\Gamma_U} - [v_N, T_N - S_N]_{\Gamma_S} \ge 0 \quad \forall \tau \in \Sigma_1. \quad (72) $$

Inequality (72) expresses the principle of complementary virtual work for the structure (the unilateral displacements on $\Gamma_S$ are assumed as given and equal to $v_N$).

Let us introduce a stress field $\sigma_0$ which is statically admissible in the sense of $\Sigma_1$, i.e. it satisfies the equation of equilibrium and the static boundary conditions on $\Gamma_F$ and on $\Gamma_S$ in the tangential direction. Then we make the substitutions

$$ \tilde{\sigma} = \tilde{\sigma} - \sigma_0, \quad \tilde{\tau} = \tau - \sigma_0. \quad (73) $$

We verify easily that $\tilde{\sigma}, \tilde{\tau} \in \Sigma_0$ where

$$ \Sigma_0 = \{\tilde{\tau} | \tilde{\tau} = \{\tilde{\tau}_{ij}\}, \tilde{\tau}_{ij} = \tilde{\tau}_{ji}, \tilde{\tau}_{ij,j} = 0 \text{ in } \Omega, \tilde{T}_i = 0 \text{ on } \Gamma_F, \tilde{T}_{\tau_i} = 0 \text{ on } \Gamma_S\}. \quad (74) $$

Then inequality (72) takes the form: find $\tilde{\sigma} = \tilde{\sigma}(v) \in \Sigma_0$ such as to satisfy

$$ A(\tilde{\sigma} + \sigma_0, \tilde{\tau} - \tilde{\sigma}) - [U, \tilde{T} - \tilde{S}]_{\Gamma_U} - [v_N, \tilde{T}_N - \tilde{S}_N] \ge 0 \quad \forall \tilde{\tau} \in \Sigma_0. \quad (75) $$

Here $\tilde{S}, \tilde{T}$ are the boundary tractions corresponding to $\tilde{\sigma}, \tilde{\tau}$, respectively. Due to the fact that $\Sigma_0$ is a linear space we may easily show (set $\tilde{\tau} - \tilde{\sigma} = \pm \phi \in \Sigma_0$) that inequality (75) is equivalent to the variational equality (i.e. to a classical bilateral problem)

$$ A(\tilde{\sigma}, \tilde{\tau}) - [U, \tilde{T}]_{\Gamma_U} - [v_N, \tilde{T}_N]_{\Gamma_S} + A(\sigma_0, \tilde{\tau}) = 0 \quad \forall \tilde{\tau} \in \Sigma_0. \quad (76) $$

We note that

$$ A(\sigma_0, \tilde{\tau}) = \int_{\Omega} \varepsilon_{0ij} \tilde{\tau}_{ij} \, d\Omega = - \int_{\Omega} \tilde{u}_{0i} \tilde{\tau}_{ij,j} \, d\Omega + [\tilde{u}_{0N}, \tilde{T}_N]_{\Gamma_S} + [\tilde{u}_{0T}, \tilde{T}_T]_{\Gamma_S} \\ + [\tilde{u}_0, \tilde{T}]_{\Gamma_F} + [\tilde{u}_0, \tilde{T}]_{\Gamma_U} = [\tilde{u}_{0N}, \tilde{T}_N]_{\Gamma_S} + [\tilde{u}_0, \tilde{T}]_{\Gamma_U} \quad \forall \tilde{\tau} \in \Sigma_0 \quad (77) $$

where $\varepsilon_0 = c \cdot \sigma_0$ and $\tilde{u}_0$ is a displacement field corresponding to $\sigma_0$. Note that we are free to assume that $\sigma_0$ is the unique solution of a bilateral problem having on $\Gamma_U$ and on $\Gamma_S$ zero displacements. In this case eqn (77) takes the form

$$ A(\sigma_0, \tilde{\tau}) = 0 \quad \forall \tilde{\tau} \in \Sigma_0. \quad (78) $$

If $\sigma_0$ is such that on $\Gamma_U$ eqn (71) is satisfied and on $\Gamma_S$, $\tilde{u}_{0N} = 0$, then

$$ A(\sigma_0, \tilde{\tau}) = [U, \tilde{T}]_{\Gamma_U} \quad \forall \tilde{\tau} \in \Sigma_0. \quad (79) $$

Further we follow the general case of eqn (77).

Then $\tilde{\sigma}$ in eqn (76) can be written generally as the sum $\tilde{\sigma}_1 + \tilde{\sigma}_2$ where $\tilde{\sigma}_1$ and $\tilde{\sigma}_2$ are solutions of
<figure><img src="image_5.png" /><figcaption>Fig. 4. The problem decomposition (displacement method).</figcaption></figure>

$$A(\bar{\sigma}_1, \bar{\tau}) - [v_N, \bar{T}_N]_{\Gamma_S} = 0 \quad \forall \bar{\tau} \in \Sigma_0 \quad (80)$$

$$A(\bar{\sigma}_2, \bar{\tau}) - [U - \tilde{u}_0, \bar{T}]_{\Gamma_v} + [\tilde{u}_{0N}, \bar{T}_N]_{\Gamma_S} = 0 \quad \forall \bar{\tau} \in \Sigma_0 \quad (81)$$

respectively. Both eqns (80) and (81) are respectively expressions of the “principle” of complementary virtual work for fictive bilateral structures resulting from the previous one in the following way: for eqn (80) (resp. eqn (81)) we assume a structure $\Omega$ under the action of “given” displacements $v_N$ (resp. $-\tilde{u}_{N0}$) on $\Gamma_S$, zero forces in $\Omega$, on $\Gamma_F$ and tangentially on $\Gamma_S$, and zero (resp. $U - \tilde{u}_0$) displacements on $\Gamma_U$. Due to the fact that both variational equalities correspond to bilateral structures, $\bar{\sigma}_1$ and $\bar{\sigma}_2$ are uniquely determined[12]. Obviously if eqn (79) holds, eqn (76) takes the simplified form of eqn (80), i.e. $\bar{\sigma} = \bar{\sigma}_1$. From eqns (80) and (81) we obtain that

$$\bar{\sigma}_1 = H(v_N), \quad \bar{\sigma}_2 = H(U, \tilde{u}_0), \quad \bar{\sigma}_1 + \bar{\sigma}_2 = \bar{\sigma}, \quad \bar{\sigma} \in \Sigma_0 \quad (82)$$

where $H$ denotes the Green operator corresponding to the boundary value problems (80) and (91).

For both bilateral structures corresponding to eqns (80) and (81), $H$ is the same because for both structures the “same” boundary conditions hold (Fig. 4).

We now combine eqns (71) and (73), then from eqns (77) (for $\bar{\tau} = \bar{\sigma}$), (80) (for $\bar{\tau} = \bar{\sigma}_1$), and (81) (for $\bar{\tau} = \bar{\sigma}$ and $\bar{\tau} = \bar{\sigma}_1$), we obtain that

$$
\begin{align}
\tilde{\Pi}_2(v_N) &= L(v_N, \bar{\sigma}) = [\tilde{S}, U]_{\Gamma_v} + [\tilde{S}_N, v_N]_{\Gamma_S} - \frac{1}{2}A(\bar{\sigma}, \bar{\sigma}) \nonumber \\
&= [\tilde{S} + S_0, U]_{\Gamma_v} + [\tilde{S}_N + S_{0N}, v_N]_{\Gamma_S} - \frac{1}{2}A(\bar{\sigma}, \bar{\sigma}) - \frac{1}{2}A(\sigma_0, \sigma_0) - A(\bar{\sigma}, \sigma_0) \nonumber \\
&= [\tilde{S}, U]_{\Gamma_v} + [\tilde{S}_N, v_N]_{\Gamma_S} - \frac{1}{2}A(\bar{\sigma}, \bar{\sigma}) - \frac{1}{2}A(\sigma_0, \sigma_0) - [u_0, \tilde{S}]_{\Gamma_v} - [u_{0N}, \tilde{S}_N]_{\Gamma_S} \nonumber \\
&\quad + [S_0, U]_{\Gamma_v} + [S_{0N}, v_N]_{\Gamma_S} \nonumber \\
&= [\tilde{S}_2, U]_{\Gamma_v} + [\tilde{S}_{2N}, v_N]_{\Gamma_S} + \frac{1}{2}[\tilde{S}_{1N}, v_N]_{\Gamma_S} - \frac{1}{2}[u_{0N}, \tilde{S}_{2N}]_{\Gamma_S} - \frac{1}{2}[u_0, \tilde{S}_2]_{\Gamma_v} \nonumber \\
&\quad - \frac{1}{2}[U, \tilde{S}_2] + [S_0, U]_{\Gamma_v} + [S_{0N}, v_N]_{\Gamma_S} - \frac{1}{2}A(\sigma_0, \sigma_0) \nonumber \\
&= [H(U, u_0), U]_{\Gamma_v} + [[H(U, u_0)]_N, v_N]_{\Gamma_S} + \frac{1}{2}[[H(v_N)]_N, v_N]_{\Gamma_S} \nonumber \\
&\quad - \frac{1}{2}[u_{0N}, [H(U, u_0)]_N]_{\Gamma_S} - \frac{1}{2}[u_0, H(U, u_0)]_{\Gamma_v} - \frac{1}{2}[U, H(U, u_0)] \nonumber \\
&\quad + [S_0, U]_{\Gamma_v} + [S_{0N}, v_N]_{\Gamma_S} - \frac{1}{2}A(\sigma_0, \sigma_0). \tag{83}
\end{align}
$$

Let us introduce the bilinear form
$$ \delta(v_N, w_N) = [[H(v_N)]_N, w_N]_{\Gamma_S} \quad (84) $$

and the linear form (H is a linear operator)

$$ \zeta(w_N) = -[[H(U, u_0)]_N + S_{0N}, w_N]_{\Gamma_S} \quad (85) $$

We denote also by $R(U, u_0, \sigma_0)$ the remaining constant terms in eqn (83). Thus

$$ \tilde{\Pi}_2(v_N) = \sup_{\tau \in \Sigma_1} L(v_N, \tau) = L(v_N, \tilde{\sigma}) = \frac{1}{2}\delta(v_N, v_N) - \zeta(v_N) + R(U, u_0, \sigma_0) \quad (86) $$

and from eqn (70) we see that we can formulate the following minimization problem with respect to the unknown boundary displacement $v_N$:

$$ \min \{ \Pi_2(v_N) = \frac{1}{2}\delta(v_N, v_N) - \zeta(v_N) | v_N \le 0 \} \quad (87) $$

Using the functional framework introduced in Section 2 we can show that the minimum problem (87) admits a solution which is uniquely determined. The proof surpasses the purposes of the present paper and is omitted here.

Further we shall show that if $u_N$ is the solution of problem (87) then $u_N$ corresponds to the stress field $\sigma$, i.e.

$$ u_N = u_N(\sigma) \quad (88) $$

where $\sigma \in \Lambda$ is the solution of the minimum complementary energy problem (22): let us choose $u_N$ such as to satisfy eqn (88), where $\sigma$ is the solution of eqn (31). Then $(u_N, \sigma) \in K \times \Sigma_1$ is the unique solution of the saddle-point problem (70) and we can verify that $u_N(\sigma)$ is a solution of the minimum problem (87). Then the uniqueness of the solution of problem (87) implies eqn (88). Again we can verify by using Betti's theorem that $\delta$ is symmetric, i.e. that (Fig. 4)

$$ \delta(v_N, w_N) = \delta(w_N, v_N). \quad (89) $$

We close this section by giving the analogous formulation to relations (51) and (52). The first is a variational inequality: find $u_N \in K$ such as to satisfy

$$ \delta(u_N, v_N - u_N) - \zeta(v_N - u_N) \ge 0 \quad \forall v_N \in K. \quad (90) $$

The second is a multivalued integral equation on the boundary part $\Gamma_S$ of the structure which reads

$$ \zeta - \frac{1}{2} \operatorname{grad} \delta(u_N, u_N) \in \partial I_k(u_N). \quad (91) $$

Each one of the last two formulations is equivalent to the minimum problem (87).

Note that if we choose as $\sigma_0$ a stress field satisfying eqn (78) or (79) then the expressions in eqn (81) are simplified. We get again a quadratic minimization problem of the form of (87). We arrive at the same result if the unilateral contact boundary conditions with a deformable support hold. In this case solutions (80) are obtained with respect to the enlarged structure by adding the appropriate linear springs.

## 5. DISCRETIZATION OF THE PROBLEM—REMARKS AND NUMERICAL APPLICATIONS

As we have shown in the previous sections two minimum problems hold on the part of the boundary of the body subjected to the unilateral contact conditions, problem (43) with respect to the unknown reactions and problem (87) with respect to the unknown displacements. They have been derived from appropriate Lagrangian formulations of the
problem. Obviously the same method can be applied to discretized structures and minimum principles analogous to problems (43) and (87) can be obtained. It should be noted here that in the framework of the discretized theory such minimum problems with respect to the unknown boundary displacements have already been formulated through elimination of the internal degrees of freedom by means of the superelement technique[17]. An analogous variational principle has been obtained in Ref. [7] for a discretized structure using a Cholesky decomposition technique. For a continuous structure the minimum problems (43) and (87) have to be discretized (see in this respect also Ref. [15]). To illustrate this procedure let us discretize problem (43) for a plane polygonal elastic body: we assume that $\Omega \subset \mathbb{R}^2$ is a polygonal domain and let for simplicity $U = 0$ on $\Gamma_U$. Let $\{D_h\}$, $h \to 0$ positive, be a family of finite element discretizations of $\Omega$. To any finite element discretization $d_h \in \{D_h\}$ we attribute a finite dimensional subspace $\tilde{V}_h$ where

$$ \tilde{V}_h = \{v_h | v_{hi} \in C(\Omega), v_{hi} \text{ over } T \in P(T) \ \forall T \in d_h, v_h = U \text{ on } \Gamma_U\}. \quad (92) $$

Here $C(\Omega)$ denotes the space of continuous functions on $\Omega$ and $P(T)$ denotes a space of polynomials over the finite element $T$. The degree of the polynomial depends on the finite element scheme chosen. We introduce further a family $\{D_H\}$, $H \to 0$ positive, of partitions on $\Gamma_S$. The nodes of $\{D_H\}$ generally do not need to coincide with the nodes of $\{D_h\}$ on $\Gamma_S$. If they coincide we write symbolically $h \equiv H$. Let for instance $P(T)$ be the space of linear polynomials and let us assume a family of triangular finite elements, regular with respect to their angles (cf. p. 138 in Ref. [18]). Let $\{D_H\}$ be a regular family of partitions on $\Gamma_S$ in the sense that, if $a_i^{(H)}$, $a_2^{(H)}$, ..., $a_m^{(H)}$ are the nodes of the partition $D_H$, then there exists $c > 0$ such that $\min H_i \max H_i > C$ where $H_i$ is the length of $a_i^{(H)}$ $a_{i+1}^{(H)}$. If $P_0(a_i^{(H)}, a_{i+1}^{(H)})$ is the set of all constant functions defined on $a_i^{(H)} \cdot a_{i+1}^{(H)}$ then we introduce the space

$$ E_H = \{ \mu_{HN} | \mu_{HN} \text{ over } a_i^{(H)} a_{i+1}^{(H)} \in P_0(a_i^{(H)}, a_{i+1}^{(H)}), \ i = 1, 2, \dots, m \} \quad (93) $$

and the set

$$ L_H = \{ \mu_{HN} | \mu_{HN} \in E_H, \ \mu_{HN} \le 0 \text{ on } \Gamma_S \}. \quad (94) $$

It is obvious, that now we are in the position to formulate directly, by means of Galerkin's method the finite dimensional problem corresponding to problem (43): find $\lambda_{HN} \in L_H$ such as to minimize $\Pi_1(\mu_{HN})$ over $L_H$. However, the numerical treatment of this problem is difficult because the explicit form of Green's operator $G$ is known only in special cases. Therefore, we introduce the stiffness matrix $K_h$, which results for the discretization of $\Omega$, or equivalently, by formulating the bilinear form $a(u, v)$ on $\tilde{V}_h$. The inverse $K_h^{-1}$ plays the role of $G$ in the approximation scheme and thus the finite dimensional problem reads: find $\lambda_{HN} \in L_H$ such that

$$ \Pi_{1h}(\lambda_{HN}) = \min \{ \Pi_{1h}(\mu_{HN}) | \mu_{HN} \in L_H \} \quad (95) $$

where

$$ \Pi_{1h}(\mu_{HN}) = \frac{1}{2}\mu_{HN}^T \cdot [K_h^{-1}(\mu_{HN})]_N + \mu_{HN}^T [K_h^{-1}(\bar{I}_h)]_N \quad (96) $$

and $\bar{I}_h \in \tilde{V}_h$ is the discretized form of $\bar{I}$, i.e. (equality of works)

$$ (\bar{I}_h, v_h) = (\bar{I}, v_h) \quad \forall v_h \in \tilde{V}_h. \quad (97) $$

We can easily verify that, if $(\mu_{HN}, v_{hN}) = 0 \forall v_h \in \tilde{V}_h$ implies $\mu_{HN} = 0$. $\Pi_{1h}$ is a strictly positive definite quadratic function of $\mu_{HN}$ and therefore by the well-known theorem of quadratic optimization the discretized minimum problem (95) admits a unique solution.

Analogously we may proceed for the case of the minimum problem (87). Until now we have applied a direct mathematical discretization to the minimum problem (43) and (87)
which seems to be quite cumbersome. However, a careful observation of the mechanical meaning of problems (43) and (87) and of the method for their derivation facilitates the procedure considerably and constitutes the great advantage of the proposed method.

Let us consider first the minimum problem (43). In order to calculate the discrete form of $\Pi_1$ we proceed as is obvious from problems (35), (36), (40) and (41) as follows. We consider first the structure $\Omega_h$ obtained from the initial one by assuming only the kinematical constraints on $\Gamma_v$. Then we solve this underconstrained structure for unit load $\mu_{N_1} = 1$ on the first node of $\Gamma_s$ and we obtain the corresponding normal displacements of all the $m$-nodes of $\Gamma_s$. These displacements constitute the first column of a matrix $B$. We repeat this procedure for the second node, etc. and thus we form the whole symmetric matrix $B$. The normal displacements of the nodes of $\Gamma_s$ for the same structure under the given external actions constitute a vector $g$. Then the solution of the discrete quadratic programming (QP) problem

$$ \min \left\{ \frac{1}{2} \mu^T B \mu - g^T \mu \mid \mu \le 0 \right\} \qquad (98) $$

where $\mu = \{\mu_{N_1}, \dots, \mu_{N_m}\}$, supplies the unknown normal reactions on $\Gamma_s$. Once the optimal solution is obtained the displacement and stress fields of the whole structure are calculated by back substitution.

We proceed analogously for problem (87), as it becomes obvious by considering eqns (80), (81), (84) and (85). We consider first a structure $\Omega_0$ obtained from the initial one by assuming $U_i = 0$ on $\Gamma_u$, $T_T = 0$ on $\Gamma_s$ and $T_i = 0$ on $\Gamma_f$. Then we solve this structure by imposing a unit normal displacement on the first node of $\Gamma_s$ and zero normal displacements on the other nodes of $\Gamma_s$. The solution of this kinematically overconstrained structure $\Omega'_0$ gives the corresponding normal reactions of all the $m$-nodes of $\Gamma_s$. They constitute the first column of a matrix $D$. We repeat this procedure for the second node, etc. and thus we form the whole symmetric matrix $D$. The normal reactions of the nodes of $\Gamma_s$ for a structure $\Omega''_0$ having $U_i - \tilde{u}_{0i}$ displacements on $\Gamma_u$, $T_T = 0$ on $\Gamma_s$, $T_i = 0$ on $\Gamma_f$, and $\tilde{u}_{0N}$ normal displacements on $\Gamma_s$ together with the reactions $S_{0N}$, which include the influence of the loading $C_T$, $F$, $f$ of the nodes of $\Gamma_s$, constitute a vector $z$. Then the solution of the discrete QP problem

$$ \min \left\{ \frac{1}{2} v^T D v - z^T v \mid v \le 0 \right\} \qquad (99) $$

where $v = \{v_{N_1}, \dots, v_{N_m}\}$ gives the unknown normal displacements on $\Gamma_s$. If, more specifically, $\sigma_0$ is such that on $\Gamma_u$ eqn (1) is fulfilled and on $\Gamma_s \tilde{u}_{0N} = 0$, then eqn (85) takes the simplest form

$$ \zeta(w_N) = -[S_{0N}, w_N]_{\Gamma_S} \qquad (100) $$

since $\tilde{\sigma}_2 = 0$ and thus $z$ is obtained simply from a structure loaded by the external actions (i.e. $f, C_T, F_i$), having the given displacements on $\Gamma_u$ and satisfying $\tilde{u}_{0N} = 0$ on $\Gamma_s$. It is to be noted that in the foregoing procedure we can have also as external "actions" initial strains. Note the boundary variational principles (43) and (87) exhibit a kind of inherent duality in the sense of the duality between the "force" method and "displacement" method in the classical theory of elasticity. The obtained QP problems are symmetric positive definite and have full matrices but a small number of unknowns (i.e. the normal reactions or displacements on $\Gamma_s$) compared to the QP formulation of the problem without elimination of the internal DOF (see, e.g. Refs [6, 8, 9, 14, 19]).

It is also very important to note that the matrices $D$, $B$ and the vectors $z$ and $g$ are obtained from the solution of bilateral structures under unit imposed displacements or forces by applying either the FEM or the BEM as, for instance, it is advised for infinite media. Also any analytical approach available (i.e. in plate problem diagrams for "Einflussfelder") can be used. Then the QP problems obtained, can be solved by any available QP algorithms. The whole procedure, i.e. the formulation of matrices $D$, $B$, $z$ and $g$ and
the QP algorithm, can be fully automated and no external intervention is needed as it is, Refs [8, 9].

The proposed method can be considered also as an integral “equation” method for the treatment of the unilateral contact problem. Indeed eqns (52) and (81) are multivalued boundary integral equations. An extension of the direct BEM, as it is known for bilateral problems, to the present unilateral contact problem would lead to non-symmetric linear complementarity problems (LCP). Indeed the relations connecting $S_i$ and $u_i$ on $\Gamma$ involve non-symmetric matrices (cf. e.g. p. 98 in Ref. [20]) and these relations combined with relations (5) do not allow formulation of a symmetric LCP. Analogous results are obtained with the indirect BEM (cf. also in this context p. 160 in Ref. [2]). The numerical solution of non-symmetric LCPs presents considerable difficulties due to the lack of well functioning algorithms. On the other hand the multivalued integral equation approach as presented here, leads to symmetric positive definite LCPs: indeed the minimum problem (98) (resp. (99)) is equivalent to the LCP (cf. p. 352 in Ref. [2])

$$ B\mu - g \le 0, \quad \mu \le 0, \quad \mu^T(B\mu - g) = 0 \qquad (101) $$

resp.

$$ Dv - z \le 0, \quad v \le 0, \quad v^T(Dv - z) = 0. \qquad (102) $$

Therefore, besides the well developed QP algorithms also the corresponding LCP algorithms can be applied for the numerical treatment. It is worth noting also that boundary minimum problems can be formulated for all classes of monotone unilateral constraints on the boundary, such as, e.g. friction problems, plastic hinge problems, and in the interior of the body, such as e.g. in plasticity, locking materials, etc., with the difference that in these cases (98) and (99) involve also additional convex non-differentiable energy functionals. This theory will be presented elsewhere.

As a first example we have treated the space frame of Fig. 5. The springs denote the unilateral constraints. Besides the external forces beams 51-55 of the frame are subjected to a uniform temperature distribution $t_s = 30°C$.

The minimum problem (98) was solved by the iterative algorithm of Hildreth and d'Esopo. Matrix B is a full symmetric 25 × 25 matrix the elements of which are obtained by applying unit loads in the directions of the unilateral constraints as already described. The corresponding bilateral structures were calculated by the STRESS program. As a result we obtain that only constraints 2, 5 and 18 (in the x-direction) are active. Note that the same structure solved by a direct QP approach[6, 14, 19] would lead to a QP problem involving a [6 × 36]² symmetric banded matrix with 6 × 36 unknowns and 25 inequality conditions, which for this example needs eight times the time needed with the present method on a PRIME 2250 computer. In the later case the use of a better QP algorithm is necessary. This could reduce computation time for the direct QP approach (but also in the proposed method as well). Analogous is the advantage of this method compared with the linearization method proposed in Refs [8, 9] which is based on the Theil and Van de Panne[21] optimality criterion; this is obvious especially for the present example, for which due to the relatively large number of unilateral reactions, user's intuition does not conduct the algorithm quickly to the solution.

As a second example we treat the plane elasticity problem of Fig. 6. Due to the fact that all the supports are unilateral we can use the minimization problem (99) and we determine the unknown unilateral displacements.

The resulting overconstrained bilateral structures are solved by a BE computer program using constant displacement boundary elements. It is to be noted that for this example the linearization procedure of Refs [8, 9] cannot be used due to the fact that all the constraints are unilateral. For the same reasons the minimization problem (100) cannot be applied as well. Finally the example of Fig. 7 was solved by both variational principles (98) and (99). The intermediate bilateral structures were solved by applying the BE technique as before. Both methods give an active constraint, constraint 15, and therefore the remaining results
<figure><img src="image_6.png" /><figcaption>Fig. 5. First example: (a) data of the problem, displacement and bending moment diagrams of the lateral frame without (in (b) and (c)) and with (in (d) and (e)) contact springs.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

b.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Second example: (a) data of the problem; (b) final deformation scheme.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Third example: (a) data of the problem; (b) stress diagrams.</figcaption></figure>
are the same. It is obvious that the proposed method has great advantages for large structures with a relatively small number of unilateral constraints. Moreover, the resulting QP problems with a small number of unknowns which must be solved do not compel us to use a modern and powerful QP solver; we content ourselves with a simple and relatively old iterative algorithm such as the Hildreth and d'Esopo algorithm. Finally, it should be noted that the whole method due to its structure based on Lagrangian formulations (cf. also Ref. [22]) is very well suited for parallel computation (calculations of the matrices of (98) or (99); use of an iterative QP algorithm), which does not happen in this extent for the direct QP approach and the linearization method.

## REFERENCES

1. G. Duvaut and J. L. Lions, *Les inéquations en Mécanique et en Physique*. Dunod, Paris (1972).

2. P. D. Panagiotopoulos, Inequality problems in mechanics and applications. In *Convex and Nonconvex Energy Functions*. Birkhäuser, Boston (1985).

3. J. J. Moreau, La notion de sur-potentiel et les liaisons unilatérales en élastostatique. C.R. Acad. Sci. Paris 267A, 954 (1968).

4. G. Fichera, The Signorini elastostatic problem with ambiguous boundary conditions. Proc. Int. Conf. Application of the Theory of Functions in Continuum Mechanics, Vol. I. Tbilisi (1963).

5. G. Fichera, Problemi elastostatici con vincoli unilaterali: il problema di Signorini con ambigue condizioni al contorno. Mem. Accad. Naz. Lincei VIII(7), 91 (1964).

6. M. Frémond, Etude des structures viscoélastiques stratifiées soumises à des charges harmoniques et de solides élastiques resposant sur ces structures. Thèse de doctorat d'Etat, Univ. Paris VI (1971).

7. C. Bisbos, A Cholesky condensation method for unilateral contact problems. S.M. Archives 11, 1-23 (1986).

8. P. D. Panagiotopoulos and D. Talaslidis, A linear analysis approach to the solution of certain classes of variational inequality problems in structural analysis. Int. J. Solids Structures 16, 991 (1980).

9. D. Talaslidis and P. D. Panagiotopoulos, Linear finite element approach to unilateral contact problem in structural dynamics. Int. J. Numer. Meth. Engng 18, 1505 (1982).

10. J. Haslinger, Mixed formulation of elliptic variational inequalities and its approximation. Aplik. Mat. 26, 462 (1981).

11. J. Haslinger and I. Hlavaček, Approximation of the Signorini problem with friction by a mixed finite element method. J. Math. Analysis Applic. 86, 99 (1982).

12. I. Hlavaček, J. Haslinger, J. Nečas and J. Lovíšek, *Riešenie variačnych nerovnosti v mechanike*. ALFA, Bratislava (1982).

13. I. Ekeland and R. Teman, *Convex Analysis and Variational Problems*. North Holland, Amsterdam; American Elsevier, New York (1976).

14. P. D. Panagiotopoulos, A nonlinear programming approach to the unilateral contact and friction-boundary value problem in the theory of elasticity. Ing.-Arch. 44, 421 (1975).

15. J. Haslinger and P. D. Panagiotopoulos, The reciprocal variational approach to the Signorini problem with friction. Approximation results. Proc. R. Soc. Edinb. 98, 250 (1984).

16. I. S. Sokolnikoff, *Mathematical Theory of Elasticity*. McGraw-Hill, New York (1956).

17. Nguyen Dang Hung and Géry de Saxcé, Frictionless contact of elastic bodies by finite element method and mathematical programming technique. Comput. Struct. 11, 55 (1980).

18. G. Fix and G. J. Strang, *An Analysis of the Finite Element Method*. Prentice Hall, Englewood Cliffs, New Jersey (1973).

19. G. Dupuis and A. Probst, Etude d'une structure élastique soumise à des conditions unilatérales. J. Méc. 6, 1 (1967).

20. P. K. Banerjee and R. Butterfield, *Boundary Element Methods in Engineering Science*. McGraw-Hill, New York (1981).

21. H. Künzi and W. Krelle, *Nichtlineare Programmierung*. Springer, Berlin.

22. O. C. Zienkiewicz, J. P. Vilotte and S. Toyoshima, Iterative method for constrained and mixed approximation. An inexpensive improvement of FEM performance. Comp. Meth. Appl. Mech. Engng 51, 3 (1985).