# Multigrid methods for unilateral contact problems with friction

Frédéric Lebon, Michel Raous, Iulian Rosu

► To cite this version:

Frédéric Lebon, Michel Raous, Iulian Rosu. Multigrid methods for unilateral contact problems with friction. IUTAM-Symposium on Computational Methods in Contact Mechanics, Springer, pp.1-16, 2007, 10.1007/978-1-4020-6405-0_1. hal-00476788

HAL Id: hal-00476788

https://hal.science/hal-00476788v1

Submitted on 23 Mar 2021

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 CC BY 4.0 - Attribution - International License
# Multigrid methods for unilateral contact problems with friction

Frédéric Lebon, Michel Raous, Iulian Rosu

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

## To cite this version:

Frédéric Lebon, Michel Raous, Iulian Rosu. Multigrid methods for unilateral contact problems with friction. IUTAM-Symposium on Computational Methods in Contact Mechanics, Springer, pp.1-16, 2007. [hal-00476788](https://hal.archives-ouvertes.fr/hal-00476788)

HAL Id: hal-00476788

https://hal.archives-ouvertes.fr/hal-00476788

Submitted on 23 Mar 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.
# Multigrid Methods for Unilateral Contact Problems with Friction

Frédéric Lebon¹,², Michel Raous¹ and Iulian Rosu¹

¹Laboratoire de Mécanique et d'Acoustique – CNRS, 13402 Marseille Cedex 20, France
E-mail: {lebon, raous, rosu}@lma.cnrs-mrs.fr

²UFR Mathématiques, Informatique, Mécanique, Université de Provence, France

**Abstract.** In this paper, multigrid methods are tested on unilateral problems with friction. An optimal strategy is presented and efficiency of the solver is discussed on several examples.

**Key words:** Unilateral contact, friction, multigrid.

## 1 Introduction

Multigrid methods have been widely used in fluids mechanics when large numbers of degrees of freedom are involved. Usually the geometries are sufficiently simple to enable the generation of multiple overlapped meshes in an easy way (essentially in the context of finite difference methods). In nonlinear structure mechanics, the computational costs increase because of the treatment of nonlinearities and finite elements methods are dominant because of the complexity of the geometries. The present work investigates the ability of multigrid methods to reduce the computational times and analyzes the specific problems of formulation and implementation related to the treatment of nonlinearities in the context of finite element methods. This work is conducted on contact problems involving unilateral contact and friction between an elastic body and a rigid obstacle. The nonlinearities are stiff because the contact behavior laws are nonsmooth (the nonpenetration is characterized by the nonregularized Signorini conditions) and nondifferentiable because of the use of the nonregularized Coulomb law.

## 2 The Frictional Contact Problem

### 2.1 The Contact Model

Details of the formulation are to be found in [20, 28, 28, 30]. Considering 2D problems, let $\Omega$ be an open bounded set of $\mathbb{R}^2$, which is the interior of an elastic obstacle, with
a sufficiently smooth boundary $\Gamma = \Gamma_F \cup \Gamma_D \cup \Gamma_C$, in contact with a rigid body. The boundary $\Gamma_D$ (respectively $\Gamma_F$) is the part of $\Gamma$ on which the displacements (respectively the forces) are prescribed. Part $\Gamma_C$ is that on which the contact is possible with the rigid obstacle. Let $\sigma$ be the stress tensor, $K$ the elasticity tensor, $\varepsilon$ the strain tensor and $u$ the displacement.

* Unilateral contact conditions

On the contact boundary $\Gamma_C$, the displacement $u$ and the unknown contact force $F$ are written in local coordinates (n,t) where $n$ denotes the outward normal unit vector to $\Gamma_C$:

$$ \sigma \cdot n = F \quad \text{on } \Gamma_C \qquad (1) $$

$$ u = u_N n + u_T \qquad (2) $$

$$ F = F_N n + F_T \qquad (3) $$

On $\Gamma_C$, the unilateral conditions are written as follows where $d$ denotes the initial gap:

$$ u_N \le d \qquad (4) $$

$$ F_N \le 0 \qquad (5) $$

$$ (u_N - d) \cdot F_N = 0 \qquad (6) $$

This is a Signorini problem. Neither regularization nor penalty are used.

* The Coulomb friction law

It is written as follows where $\mu$ denotes the friction coefficient:

$$ \|F_T\| \le \mu |F_N| \quad \text{with if } \|F_T\| < \mu |F_N| \quad \text{then } \dot{u}_T = 0 \qquad (7) $$

$$ \begin{array}{l} \text{if } \|F_T\| = \mu |F_N| \text{ then } \dot{u}_T \text{ is colinear} \\ \text{and opposite to } F_T \end{array} \qquad (8) $$

As presented in [6, 29, 30], by introducing an incremental formulation, it was established that the problem to be solved at each step, is subsequently very similar to the static one. In that case, the unknowns are the displacement increments and there are extra terms both in the functional which has to be minimized and in the convex characterizing the constraints. In this paper, we consider only the static problem using the following displacement formulation for the friction ((9) and (10)). Formulation and numerical algorithms can easily be extended to quasistatic evolutive problems [6, 29]. As for the unilateral contact, no regularization (compliance) is used.

$$ \|F_T\| \le \mu |F_N| \quad \text{with if } \|F_T\| < \mu |F_N| \quad \text{then } u_T = 0 \qquad (9) $$

$$ \begin{array}{l} \text{if } \|F_T\| = \mu |F_N| \text{ then } u_T \text{ is colinear} \\ \text{and opposite to } F_T \end{array} \qquad (10) $$
## 2.2 Variational Formulation

The variational formulation of the Coulomb problem can be written as an implicit variational inequality [8]. By introducing a fixed point problem $P_1$ on the sliding limit $g$ [27, 29], the problem can be solved as a sequence of classical variational inequalities including nevertheless a nondifferentiable term. Because of the symmetry of the elasticity mapping, these problems can be solved as minimization problems under constraints (problem $P_2$).

**Problem $P_1$:** Find the sliding limit function $G$, fixed point of the application $g \to \mu|F_N(u_g)|$ where $u_g$ is solution to the problem $P_2$ depending on $g$.

**Problem $P_2$:** Let $\phi_1$ be the load applied in $\Omega$, $\phi_2$ be the load applied to part $\Gamma_F$ of the boundary, and $g$ be given,

find $u_g \in \mathbb{K} = \{v, v = 0 \text{ on } \Gamma_D \text{ and } v_N \le d \text{ on } \Gamma_C\}$ such that:

$$J(u_g) \le J(v) \quad \forall v \in \mathbb{K} \quad (11)$$

with

$$J(v) = \frac{1}{2}a(v, v) - (f, v) + j(v) \quad (12)$$

$$\text{where: } a(u, v) = \int_{\Omega} \text{grad}_s u K \text{grad}_s v \, dx \quad (13)$$

$$j(v) = \int_{\Gamma_C} g \|v_T\| \, dl \quad (14)$$

$$(f, v) = \int_{\Omega} \phi_1 v \, dx + \int_{\Gamma_F} \phi_2 v \, dl \quad (15)$$

## 3 Preconditioner for Mixed Formulation

In this section, we adopt the mixed formulation proposed in [1]. The problem is split into two parts involving the pair $x = (u, \lambda)$, where $u$ denotes the displacement vector and $\lambda$ the contact nodal forces, say a differentiable part $G$ and a nondifferentiable one $F$:

$$G(x) + F(x) = 0 \quad (16)$$

This last problem is solved using the Generalized Newton Method (GNM):

$$x^{i+1} = x^1 - (K_i + J_i)^{-1}(G(x^i) + F(x^i)) \quad (17)$$
where $K_i = \partial G(x^i)$ and $J_i \in \partial F(x^i)$. The tangent matrix $A = K_i + J_i$ is nonsymmetric (with zeros on the diagonal) and the previous linear system can be solved by the Conjugate Gradient Squared Method (CGS) associated with the coarse/fine preconditioner proposed in [2]. The basic idea of the preconditioner is to assume that the tangent matrix $A$ has been obtained through the refinement of a coarser mesh. The matrix $A$ is split into four parts (index $c$ refers to coarse level and $f$ to the fine one):

$$A = \begin{bmatrix} A_{ff} & A_{fc} \\ A_{cf} & A_{cc} \end{bmatrix} \qquad (18)$$

A block factorization of matrix $A$ gives:

$$A = \begin{bmatrix} A_{ff} & 0 \\ A_{cf} & A_{cc} - A_{cf}A_{ff}^{-1}A_{fc} \end{bmatrix} \begin{bmatrix} I & A_{ff}^{-1}A_{fc} \\ 0 & I \end{bmatrix} \qquad (19)$$

We postulate the following preconditioner:

$$C = \begin{bmatrix} A_{ff} & 0 \\ A_{cf} & A_{cc}^* \end{bmatrix} \begin{bmatrix} I & A_{ff}^{-1}A_{fc} \\ 0 & I \end{bmatrix} \qquad (20)$$

where $A_{cc}^*$ is the tangent matrix on the coarse grid. The coarse/fine preconditioner requires to solve two sub-systems with the matrix $A_{ff}$. It is possible to define an approximated LDU factorization of this matrix associated to each contact status [2]. This preconditioner is very efficient because it is not necessary to find an optimal parameter.

## 4 Formulation and Implementation of the Multigrid Method for Contact Problems

Classically, in the linear case, the efficiency of multigrid methods can be explained in terms of the spatial frequencies of the error between the solution and the iterates when this error is written in the basis of the eigenvectors of the discretized mapping. The smoothers (Gauss-Seidel method for example) are known to quickly reduce the high frequencies of the error, but a large number of iterations is required to reduce the low components. Multigrid methods increase the convergence on the low frequencies of the error. They correct the solution obtained after only a few smoothings on the fine grid (where the solution is needed) by a correcting error computed on a coarse grid with the current defect of equilibrium. The sequence is then composed of a few smoothings of the solution on the fine grid (large system) and of complete resolution to determine a correction error on a coarse grid (small system). Appropriate interpolation and restriction have to be carried out between the grids.
With linear problems, the problem to be solved on the coarse grid is similar to the initial one defined on the fine grid. With nonlinear problems, one of the difficulties is to write an appropriate problem on the coarse grid, which is defined on the error and associated with the defect of equilibrium computed on the fine grid and carried onto the coarse one. Here, several alternatives arise. Section 4.2.1 is devoted to the problem written on the error. In Section 4.2.2 different methods of carrying the non-linearity information needed onto the coarse grid from the fine one are proposed and tested.

## 4.1 Basic Solvers

### 4.1.1 The Discrete Problem

Finite element discretization of problems $P_1$ and $P_2$ gives the following problems $P_3$ and $P_4$ written in $\mathbb{R}^{2N_h}$, where $N_h$ is the number of fine mesh nodes. Let $\mathbb{K}_h \subset \mathbb{R}^{2N_h}$ be the cone defined by:

$$ \mathbb{K}_h = K_{h,1} \times K_{h,2} \times \dots \times K_{h,2N_h} $$

Let $I_h$ be the set of indexes of the equations relating to the normal components of the contact nodes. Let $L_h$ be the set of suffixes of the tangential components of the nodes of $\Gamma_C$.

$$ K_{h,i} = ]-\infty, d] \text{ if } i \in I_h $$

$$ K_{h,i} = \mathbb{R} \quad \text{otherwise} $$

**Problem P3:** find the sliding limit $G_h$ fixed point of the application $g_h \to \mu|F_{h,N}(u_h)|$ where $u_h$ is the solution to the problem $P_4$ depending on $g_h$.

**Problem P4:** find $u_h \in \mathbb{K}_h$ such that:

$$ J_h(u_h) \leq J_h(v_h) \quad \forall v_h \in \mathbb{K}_h \quad (21) $$

$$ \text{where:} \qquad J_h(v_h) = \frac{1}{2} v_h^T A_h v_h - v_h^T f_h + |v_h|^T g_h \quad (22) $$

with

$$ |v_h|^T = (|v_{h,1}|, |v_{h,2}|, \dots, |v_{h,2N_h}|) \quad (23) $$

$$ |v_h|^T g_h = \sum_{i \in L_h} |v_{h,i}| g_{h,i} \quad (24) $$

and

- $A_h = (a_{h,ij})$ is the matrix $(2N_h \times 2N_h)$ associated to the elasticity operator,
- $f_h = (f_{h,ij})$ is the load vector,

- $g_h$ is the vector associated to the sliding limit $g$ ($g_h \in \mathbb{R}^{\text{card}(L_h)}$).

The discrete problem $P_4$ on the fine grid has been derived from an appropriate form of the approximation problem given in [23].

The problem $P_3$ is solved using a fixed point algorithm (successive approximations). Starting with $g_h^0$ equal to zero (frictionless case), the solution to problem $P_4$ gives the displacement $u_h^\ell$ and the contact force $F_{h,N}^\ell$ at each step $\ell$ ($F_{h,N}^\ell$ is computed directly from the defect of equilibrium $A_h u_h^\ell - f_h$). We then define the iterative process as $g_h^{\ell+1} = \mu |F_{h,N}^\ell|$.

### 4.1.2 Gauss-Seidel Algorithm for the Contact

Problem $P_4$, where the sliding limit $g_h^\ell$ is given by the fixed point method, is solved by an internal multigrid process. The first step in the multigrid process is the smoothing step for which a Gauss-Seidel method is often used. Here we consider a Cryer-Christopherson method [7] adapted to deal with the nondifferential term due to friction [21, 29, 30]. It is a Successive Over-Relaxation method with Projection (SORP):

$$u_{h,i}^{k+1} = P_{K_{h,i}}((1-\omega)u_{h,i}^k + \omega m) \quad (25)$$

where

$$\left\{
\begin{array}{ll}
m = \frac{1}{a_{h,ii}} \left( f_{h,i} - \sum_{j=1}^{i-1} a_{h,ij} u_{h,j}^{k+1} - \sum_{j=i+1}^{2N_h} a_{h,ij} u_{h,j}^k - \operatorname{sgn}(m) g_{h,i}^\ell \right) & \text{if } m \in \mathbb{R}^* \\
m = 0 & \text{else}
\end{array}
\right.
\quad (26)$$

with $g_{h,i}^\ell = 0$ if $i \notin L_h$. The second alternative of (26) occurs when there is no compatibility between the calculated $m$ and its assumed sign, i.e., we obtain a negative (or positive) value of $m$ when its sign (denoted $\operatorname{sgn}(m)$) was taken to be positive (or negative).

$P_{K_{h,i}}$ is the projection operator on $K_{h,i}$ and $\omega$ is the relaxation coefficient ($\omega \in [1, 2]$). Because of the lack of theoretical estimates, an optimal value of the relaxation coefficient can be only obtained through dichotomic research and this Over-Relaxation method will be used only when the method acts as a solver on the coarse grid. When the method acts as a smoother on the fine grid, $\omega$ is set equal to 1 and the smoother is then the Symmetric Gauss-Seidel method with Projection (SGSP) [9].

### 4.1.3 Mathematical Programming Solver and Complementary Formulation

The frictional contact problem introduced in Section 2 can also be formulated, after discretization, as a complementarity problem [5, 9, 16, 29, 30]).
Problem P5: let f be given, find the displacement u and the contact force F such that:

$$
\left\{
\begin{array}{ll}
A_h^* u_h = f_h^* + F_h & \\
\text{for } i \in I_h \text{ (normal components)} & \\
\left\{
\begin{array}{l}
u_{h,i} \le 0 \\
F_{h,i} \le 0 \\
u_{h,i} F_{h,i} = 0
\end{array}
\right. & (27) \\
\text{for } i \in L_h \text{ (tangential components)} & \\
\left\{
\begin{array}{l}
u_{h,i} = \lambda_2 - \lambda_1 \\
F_{h,i} = \frac{1}{2}(\varphi_1 - \varphi_2) \\
\lambda_i \le 0 \\
\varphi_i \le 0 \quad i = 1, 2 \\
\lambda_i \varphi_i = 0
\end{array}
\right. &
\end{array}
\right.
$$

Problem P5 has been reduced to the contact variables by using a condensation pro-
cedure described in [29]. Matrix A\*_h is the full sized condensed stiffness matrix which
turns out to be nonsymmetric because of the friction, and f\*_h is the condensed load-
ing force vector. The tangential contact displacement has been separated into positive
(forward) sliding λ₁ and negative (backward) sliding λ₂ relatively to the unit vector
t. This problem is known as a linear complementary problem, the general form of
which is written:

$$
w - Mz = q \tag{28}
$$

$$
w \ge 0 \tag{29}
$$

$$
z \ge 0 \tag{30}
$$

$$
wz = 0 \tag{31}
$$

We solve this problem by using the Lemke method which is a direct pivoting al-
gorithm. Details can be found in [5, 9, 16, 29].

4.2 The Two Grid Method

4.2.1 The Algorithm

We consider a fine mesh (with N<sub>h</sub> nodes) and a coarse mesh (with N<sub>H</sub> nodes). We consider the operator (prolongation) p<sub>H</sub><sup>h</sup>: R<sup>2N<sub>H</sub></sup> &rarr; R<sup>2N<sub>h</sub></sup> and the operator (restriction) p<sub>h</sub><sup>H</sup>: R<sup>2N<sub>h</sub></sup> &rarr; R<sup>2N<sub>H</sub></sup>.

Starting with $u_h^0 \in \mathbb{K}_h$, we define an iterative process:

$$
u_h^k \rightarrow u_h^{k+1} \tag{32}
$$
### Smoothing on the fine grid

We carry out v smoothing iterates (two or three iterates) of the SGSP^v procedure (Symmetric Gauss-Seidel with Projection). We denote:

$$u_h^{k+1/2} = \text{SGSP}_\omega^v(u_h^k, f_h, A_h) \qquad (33)$$

We introduce the variable $e_h^k = u_h - u_h^{k+1/2}$, $u_h$ being the (unknown) exact solution of problem P_4. Now, in order to evaluate the correction, we have to solve the new problem:

Problem P_6: find $e_h^k$ (such that $e_h^k + u_h^{k+1/2} \in \mathbb{K}_h$) such that:

$$J_h^C(e_h^k) \leq J_h^C(w_h) \quad \forall w_h \text{ with } w_h + u_h^{k+1/2} \in \mathbb{K}_h \qquad (34)$$

where

$$J_h^C(w_h) = \frac{1}{2} w_h^T A_h w_h - w_h^T (f_h - A_h u_h^{k+1/2}) + |w_h + u_h^{k+1/2}|^T g_h^\ell \qquad (35)$$

Problem P_6 now has to be set on the coarse grid. Let $r_h^k = f_h - A_h u_h^{k+1/2}$ and $r_H^k = p_H^H r_h^k$.

### Coarse grid correction

Three difficulties occur in treating the nonlinearities of contact problems with friction:

- how to carry the nonpenetration condition from the fine onto the coarse grid,

- how to treat the absolute value on the coarse grid,

- how to transfer the sliding limit $g_h^\ell$ from the fine grid onto the coarse one?

These points will be discussed in Section 4.2.2.

We take $U_H^k$ a vector chosen in $\mathbb{R}^{2N_H}$ and we set:

$$\mathbb{K}_H^C = K_{H,1}^C \times K_{H,2}^C \times \dots \times K_{H,2N_H}^C$$

with

$$\left\{ \begin{array}{ll} K_{H,i}^C = \mathbb{R} & \text{if } i \notin I_H \\ K_{H,i}^C = ]-\infty, d - U_{H,i}^k] & \text{if } i \in I_H \end{array} \right. \qquad (36)$$

We define two vectors $G_H^\ell \in \mathbb{R}^{\text{card}(\mathbb{L}_H)}$ and $V_H \in \mathbb{R}^{2N_H}$. As with classical nonlinear multigrid methods, these vectors are chosen to approximate the values on the coarse grid. A new problem (similar to eq. (35)) is written on the coarse grid:

Problem P_7: find $e_H^k \in \mathbb{K}_H^C$ such that:
$$J_H^C(e_H^k) \le J_H^C(w_H) \quad \forall w_H \in \mathbb{K}_H^C \quad (37)$$

where $J_H^C(w_H) = \frac{1}{2}w_H^T A_H w_H - w_H^T r_H + |w_H + V_H^k|^T G_H^\ell$.

How to choose $U_H^k$, $V_H^k$ and $G_H^\ell$ will be discussed in Section 4.2.2.

Problem $P_7$ involving the error $e_H^k$ has to be solved (the loads are replaced by the defect $r_H^k$ in the second member). With linear problems, it is natural to use a direct method (Cholesky). With contact problems, the SORP method or the Lemke method (linear programming method) will be used. When the Lemke method is used, $P_7$ must be changed into a complementary problem ($P_7^b$) similar to $P_5$ [9, 29]. The solution $e_H$ is interpolated onto the fine grid:

$$e_h^k = p_h^H e_H^k \quad (38)$$

The previous solution $u_h^{k+1/2}$ is corrected (we must ensure that $u_h^{k+1}$ stays within the convex $\mathbb{K}_h$):

$$u_h^{k+1} = P_{\mathbb{K}_h}(u_h^{k+1/2} + e_h^k) \quad (39)$$

and we return to Equation (33) until a given degree of precision is achieved.

### Prolongation and restriction operators

The prolongation $p_H^h$ is a linear interpolation. Consistence of the formulation is insured by choosing $p_H^h = (p_H^h)^T$. With three-node triangles when each element is divided into four triangles, we have

$$r_H^k(x) = p_H^H r_h^k(x) = r_h^k(x) + \frac{1}{2} \sum_{y \in \mathcal{C}_x} r_h^k(y).$$

$\mathcal{C}_x$ is the set of nodes connected to $x$ [20].

### The initial condition (full multigrid)

The two grid process will be initiated by solving first the problem on the coarse grid (problem $P_4$) and then using the prolongation of this solution $u_H^0$ on the fine grid as the initial condition $u_h^0 = p_H^h u_H^0$ to start the first smoothings. This is known as the full multigrid method. Numerical tests have confirmed that the choice of the full multigrid method turns out to be essential for our problem, as generally observed for nonlinear problems [9, 19].

### 4.2.2 Transfer of the Nonlinearities

#### Convex definition

The nonpenetration condition is satisfied by the condition that the solution must be within the convex $\mathbb{K}_H^k$ in problem $P_5$. On the coarse grid, the convex $\mathbb{K}_H^k$, associated with the error $e_H^k$ depends on the solution after the smoothing: it is
characterized by the condition $(U_H^k + e_H^k)_N \le d$ on $\Gamma_C$. Therefore $U_H^k$, which is not needed in linear classical case, has to be evaluated on the coarse grid. There exist various means of evaluating $U_H^k$ from the computed solution $u_h^{k+1/2}$ on the fine grid. Four of them are considered here:

(a) $U_H^k = q_h^H u_h^{k+1/2}$ where $q_h^H$ is a restriction. We can use:

- the *canonic injection* $U_H^k = u_h^{k+1/2}$,

- the *standard operator* defined in Section 4.2.1, $U_H^k = \frac{1}{2} p_h^H u_h^{k+1/2}$. Operator $p_h^H$ has been defined for the forces: it is used here on the displacements, where a coefficient 1/2 has to be introduced,

- a more sophisticated operator which involves five nodes of the fine grid.

(b) $U_H^k$ is defined as the largest normal displacement through the node under consideration and its two neighbors: $U_H^k(x) = u_h^{k+1/2}(z)$ with $z = -\operatorname{argmax}\{u_h^{k+1/2}(y), y \in C_x\}$. This is a procedure developed by Mandel [25]. It is a very good choice because it ensures that $u_h^{k+1/2} + e_h$ stays within the convex. This is convenient for theoretical reasons.

(c) $U_H^k = u_H^0$ where $u_H^0$ is the solution initially computed on the coarse grid for the full multigrid initialization defined in Section 4.3. This is derived from the *Hackbusch* [10] choice for nonlinear problems.

(d) $U_H^k = U_H^{k-1} + e_H^{k-1}$ (with the initial condition $U_H^1 = u_H^0$). This is a variant of (c), where the computed correction $e_H^k$ is applied at each iteration (*Hackbusch with updating*).

The cases (c) and (d) need the use of the nested iteration or the full multigrid method [10]. The advantage of the full multigrid method is that it provides a good starting solution at each level. The algorithm begins at the coarsest level.

* *The undifferentiable term*

The problem of the absolute value is very important from the mathematical point of view, and the term $j(v)$ is undifferentiable. The term $|e_H + V_H^k|$ depends on the choice of vector $V_H^k$. An approximation of the displacement is needed on the coarse grid. Various choices of $V_H^k$ are given in [19].

(a) $V_H^k = q_h^H u_h^{k+1/2}$ where $q_h^H$ is a restriction. We can use:

- the *standard injection* $V_H^k = \frac{1}{2} p_h^H u_h^{k+1/2}$ (tangential components),

- the *canonic injection*.

(b) the *Hackbusch choice*: $V_H^k = u_H^0$ where $u_H^0$ is the solution initially computed on the coarse grid for the full multigrid initialization. The use of a full multigrid is necessary in that case.

(c) the *Hackbusch with updating*: $V_H^k = V_H^{k-1} + e_H^{k-1}$ (with the initial condition $V_H^1 = u_H^0$).

* *The sliding limit*

$G_H^\ell$ is an approximation of the sliding limit on the coarse grid. Two cases are proposed:
<figure><img src="image_3.png" /><figcaption>Fig. 1. Plate in contact with two rigid bodies.</figcaption></figure>

(a) the *restriction* of the fine sliding limit at each multigrid iteration,

(b) the *initial value* $G_H^0$ computed with the initial resolution on the coarse grid of the full multigrid method.

## 4.3 Multigrid Method

With very large sized problems, introducing several levels of grids can keep the coarse problem small sized. This is of course important because total resolutions are conducted on the coarse grid. It is easy to generalize the previous algorithm to a MultiGrid Method (MGM) and to implement either V-cycle or W-cycle.

# 5 Optimal Strategy: Tests on the Various Alternatives

The multigrid methods have been developed in our own finite element code EU-XENE which is a specific version of GYPTIS dedicated to the contact problems at the LMA [18].

The optimal strategy for dealing with the nonlinearities is discussed in the case of a plate which is part of a pressure vessel used in nuclear engineering [21] (see Figure 1). This is an axisymmetrical problem. The structure is in contact with two rigid solids on an oblique zone (30°) and on a horizontal zone. The plate is squeezed between the two rigid obstacles by prescribing a vertical displacement (0.05 mm) of the lower one. Primary and secondary pressures are applied on the two sides of the plate. The numerical tests have been conducted by using five meshes, with 51617 nodes for the finest one and 227 nodes for the coarsest one. Of course, the solution is always computed on the finest one. The meshes are unstructured with local refinements in the contact zone (194 contact nodes on the finest mesh).
Numerical tests have been conducted on this example in order to compare the various alternatives presented in Section 4 for the implementation of multigrid methods for that nonlinear problem. The main results are given and details of the numerical results could be found in [17].

**5.1 Convex Definition**

In this case, the tests are conducted on a frictionless case. The different possible choices for the convex definition turn out to be equivalent. For mathematical reasons, the Mandel process is suitable. This is confirmed by other numerical results obtained on other examples. The Mandel procedure will be used in the computations.

**5.2 The Undifferentiable Term**

The same example is now treated with friction. The convex is defined as chosen in the previous paragraph. We test the various ways of choosing $V_H^k$ in (P$_7$) as described in Section 4.2.2.

In this case, results have shown that the use of informations from the fine grid can be very unsatisfactory when three grids are implemented. The canonic and standard injection should be avoided. The computations are conducted with 3 pre-smoothings, 2 post-smoothings and by using V-cycle. The “Hackbusch” procedure is always efficient and will be chosen in what follows.

**5.3 The Sliding Limit**

We tested the two alternatives given in Section 4.2.2 for $G_H^\ell$ in (P$_7$). The computations are conducted with 3 pre-smoothings and 2 post-smoothings. The results have shown that the Coulomb's condition $|F_T| = \mu|F_N|$ is verified with a good accuracy when the “restriction” is used to define the sliding threshold $G_H^l$. When the “initial value” is used, the accuracy varies from 5 to 10%. The choice for the coarse sliding limit is therefore the restriction of the fine sliding limit.

**5.4 Number of Smoothings**

Taking the same example, we study the effects of the number of pre- and post-smoothings on the convergence. The results have shown that it is convenient to perform both a few pre-smoothings and a few post-smoothings. We will choose in the following 3 pre-smoothings and 2 post-smoothings which turned out to be convenient choices in the numerical tests.
<table><caption>**Table 1.** Direct and multigrid resolution for two meshes of the plate with two contact zones (Lemke solver).</caption><thead><tr><th rowspan="2">Number of d.o.f</th><th>Direct method</th><th colspan="3">2-grid method</th></tr><tr><th>CPU</th><th>Iterates</th><th>Smoothings</th><th>CPU</th></tr></thead><tbody><tr><td>26018</td><td>1'43"</td><td>6</td><td>33</td><td>15.7"</td></tr><tr><td>103234</td><td>28'00"</td><td>6</td><td>33</td><td>2'16"</td></tr></tbody></table>

## 5.5 Type of Multigrid Cycle

Taking the three grid example, we have compared the results using V-cycle and W-cycle on the same example of contact with friction. The W-cycle turns out to be slightly better than the V-cycle when considering the number of resolutions on the coarse grid (number of 3-grid iterates) and the total number of smoothings. In the following W-cycles are used for the case 3-grid or more.

## 5.6 Number of Grids

The number of grids is of course an important parameter. A large number of grids is time and memory size consuming. The tests have shown that in the case of large number of degrees of freedom (very fine mesh) the choice 3 or 4 grids is efficient. For smaller sized problems 2 grids will be sufficient.

# 6 Global Efficiency of the Multigrid Method on a Contact Problem

The efficiency of the method is presented here by comparing the multigrid method to the direct resolution on two examples (other ones could be found in [17]). The example is still that of the axisymmetrical plate in contact with two rigid obstacles (see Figure 1). Results for 26018 and 103234 d.o.f. meshes are given in Table 1. The Lemke solver is used both for the direct resolution and for the resolutions on the coarse grids. The multigrid method divides the total CPU time by 8 to 10. In Figure 2, we have plotted the evolution of the CPU time relative to the number of degrees of freedom.

It has to be noted first that control of the efficiency of the multigrid methods for that nonlinear problem was not as simple as for the linear case, because among other features it depends on the contact status. So only general and qualitative comments can be given.

It is interesting to note that in nonlinear structure mechanics such as the present frictional contact analysis, the multigrid method was found to be efficient, even with medium sized discretization instead that for linear problems they are proved to be efficient only for very large sized problems. In the case of a plate with 13009 nodes,
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** CPU time versus the number of degrees of freedom for the vessel pressure plate.</figcaption></figure>

the total CPU time was already divided by 8 (Table 1) when the Lemke solver is used.

Figure 2 shows that the rate of convergence does not seem to be strongly modified by the nonlinearity of the problem. In our case no theoretical results exist. In Figure 2, the rate of convergence for the two grid method is quasi linear.

We would like to make some comments on the influence of the computer pro-
cessor properties on the multigrid efficiency. When we implemented first the mul-
tigrid methods, we have tested them on a Dec Alpha processor and also on a CISC
processor. Of course, the CPU times are drastically reduced by using the very power-
ful new processors, but it turned out that multigrid methods were more efficient with
the CISC processor than with the recent ones when the multigrid resolution is com-
pared with the direct one. It should be noted that the new generation of processors
are highly efficient at computation but still relatively slow in performing addressings
and transfers. So simple solvers such as the Gauss-Seidel and SORP ones turn out
to be more and more fast because they only do a large number of multiplications
and only a small number of transfers. That kind of solvers becomes more and more
competitive regarding sophisticated methods such as the Lemke's one and multigrid
involving a lot of transfers. So when one refers the direct computation of the solution,
it can be observed that efficiency of multigrid decreases when the processor power
increases.

**References**

1. P. Alart and A. Curnier. A mixed formulation for frictional contact problems prone to newton like solution methods. *Comp. Meth. Appl. Mech. Eng.*, 92:353–375, 1991.
2. P. Alart and F. Lebon. Solution of frictional contact problem using ilu and coarse/fine preconditioners. *Comp. Mech.*, 16:98-105, 1995.

3. A. Boersma and P. Wriggers. An algebraic multigrid solver for finite element computations in solids mechanics. *Engng. Comp.*, 14(2):202-215, 1997.

4. A. Brandt and C. W. Cryer. Multigrid algorithms for the solutions of linear complementarity problems arising from free boundary problems. *SIAM J. Sci. Statist. Computational*, 4:655-684, 1983.

5. P. Chabrand, F. Dubois, and M. Raous. Numerical methods for solving unilateral contact problems with friction. *Math. Comp. Mod.*, 28(4-8):97-108, 1998.

6. M. Cocu, E. Pratt, and M. Raous. Formulation and approximation of quasi-static frictional contact. *Int. J. Engng. Sci.*, 34(7):783-798, 1996.

7. C. W. Cryer. The solution of a quadratic programming problem using systematic overrelaxation. *SIAM J. Control*, 9(3):385-392, 1971.

8. G. Duvaut and J. L. Lions. *Inequalities in Mechanics and Physics*. Springer-Verlag, Berlin, 1976.

9. L. Grego. *Méthodes multiniveaux pour des problèmes de contact unilatéral avec frottement*. PhD thesis, Provence University, Marseille, 1995.

10. W. Hackbusch. Multigrid method of the second kind. In *Multigrid Methods*, W. Hackbush and U. Trottenberg (Eds), Lectures Notes in Math, Vol. 960. Springer-Verlag, Berlin, 1982.

11. W. Hackbusch. *Multi-Grid Methods and Applications*. Springer, Berlin, 1985.

12. W. Hackbusch and H. D. Mittelman. On multigrid methods for variational inequalities. *Numer. Math.*, 42:65-76, 1983.

13. W. Hackbusch and U. Trottenberg (Eds). *Multigrid Methods*, Lecture Notes in Math., Vol. 960. Springer-Verlag, Berlin, 1982.

14. P. W. Hemker and P. M. De Zeeuw. Some implementations of multigrid linear system solvers. In *Multigrid Methods for Integral and Differential Equations*, D. J. Paddon and H. Holstein (Eds). Clarendon Press, 1985.

15. R. H. W. Hoppe. Multigrid algorithms for variational inequalities. *SIAM J. Numer. Anal.*, 24(5):1046-1065, 1987.

16. A. Klarbring and G. Bjorkman. The treatment of problems in contact mechanics by mathematical programming. *J. Th. Appl. Mech.*, 7(1):111-128, 1988.

17. F. Lebon M. Raous L. Grego, J.-C. Latil and I. Rosu. Multigrid methods for unilateral contact problems with friction. submitted.

18. J.-C. Latil and M. Raous. Module gyptis version 1.0: Contact unilatéral avec frottement en mécanique des structures – inéquations variationnelles. *LMA Publications*, 132, 1991. Also available in INRIA-Modulef publications.

19. F. Lebon. *Résolution numérique de problème de frottement de Coulomb. Accélération de convergence par une méthode multigrilles interne*. PhD thesis, Provence University, France, 1989.

20. F. Lebon, M. A. Ledain, M. Bellet, J. L. Chenot, and M. Raous. *Introduction aux méthodes multigrilles en mécanique des structures linéaires et non linéaires*, Vol. 3. Pluralis, 1989.

21. F. Lebon and M. Raous. Multibody contact problem including friction in structure assembling. *Comp. & Struct.*, 43(5):925-934, 1992.

22. M. A. Ledain. *Méthodes multigrilles et h-adaptation: Application au calcul de structures élastiques bidimensinnelles*. PhD thesis, ENSM, Paris, 1991.
23. C. Licht, E. Pratt, and M. Raous. Remarks on a numerical method for unilateral contact including friction. *Int. Series Numer. Math.*, 101:129–144, 1991.

24. J. Mandel. Etude algébrique d'une méthode multigrille pour quelques problemes de frontière libre. *CRAS Paris*, 298(18):469–472, 1984. Série I.

25. J. Mandel. A multilevel iterative method for symmetric, positive definite linear complementarity problems. *Appl. Math. Optim.*, 11:77–95, 1984.

26. D. J. Paddon and H. Holstein (Eds). *Multigrid Methods for Integral and Differential Equations*. Clarendon Press, 1985.

27. P. D. Panagiotopoulos. *Inequality Problems in Mechanics and Application*. Birkhauser, Boston, MA, 1985.

28. Y. Pressburger and R. Perucchio. A hierarchical two-level multigrid solver. *Comp. & Struct.*, 55(3):471–483, 1995.

29. M. Raous. Quasistatic signorini problem with coulomb friction and coupling to adhesion. In P. Wriggers and P.D. Panagiotopoulos (Eds), *New Developments on Contact Mechanics*, Vol 384, CISM, Courses and Lectures, Springer, Wien, 1999, pp. 101–178.

30. M. Raous, P. Chabrand and F. Lebon. Numerical methods for frictional contact problems and applications. *J. Th. Appl. Mech.*, 7(Suppl. to 1):111–128, 1988.

31. P. Sonneveld, P. Wesseling and P. M. De Zeeuw. Multigrid and conjugate gradient methods as convergence acceleration techniques. In *Multigrid Methods for Integral and Differential Equations*, D. J. Paddon and H. Holstein (Eds). Clarendon Press, 1985.

32. K. Stuben and U. Trottenberg. Multigrid methods: Fundamental algorithms, model problem analysis and applications. In *Multigrid Methods*, W. Hackbush and U. Trottenberg (Eds). Springer-Verlag, Berlin, 1982.