# A domain decomposition strategy for nonclassical frictional multi-contact problems

Mickaël Barboteu, Pierre Alart, Marina Vidrascu

► To cite this version:

Mickaël Barboteu, Pierre Alart, Marina Vidrascu. A domain decomposition strategy for nonclassical frictional multi-contact problems. Computer Methods in Applied Mechanics and Engineering, 2001, 190 (37-38), pp.4785-4803. 10.1016/S0045-7825(00)00347-9. hal-01393140

HAL Id: hal-01393140

https://hal.science/hal-01393140

Submitted on 6 Nov 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
# A domain decomposition strategy for nonclassical frictional multi-contact problems

M. Barboteu<sup>a</sup>, P. Alart<sup>a,*</sup>, M. Vidrascu<sup>b</sup>

In this paper we present a numerical strategy to be solve large scale frictional contact problems by domain decomposition methods which are adapted to parallel computers. The motivation is given by the study of the mechanical behavior of rolling shutters composed by many hinged slats. The numerical treatment of such nonclassical contact problems leads to very large strongly nonlinear, non-symmetric and ill-conditioned systems. Domain decomposition methods are a good alternative to overcome the difficulties of classical sequential solutions. We present a nonlinear strategy adapted to problems, called "multi-contact" problems.

**Keywords:** Frictional contact; Contact element; Generalized Newton; Domain decomposition; Schur complement; GMRES; Preconditioner; Parallelism

## 1. How does contact occur in multi-body systems (?)

The mechanics of contact takes an important place in computational structural mechanics. Indeed, to obtain more information on complex systems, it is necessary to take into account imperfect joints characterized for instance by friction or clearance. We can distinguish, in particular, cases where contact is localized from cases where contact is diffuse (many contact zones). In the last case, called multi-contact problem, we identify the structure assembly or the "multi-body" systems with frictional contact conditions between the bodies or the substructures. The rolling shutters composed by many slats jointed by a hinge with play (clearance) and friction are a typical example of a multi-contact problem [2].

The aim of this paper is to present an efficient numerical scheme for this multi-body system with deformable components. To handle this problem it is necessary to develop rigorous modelling and new numerical tools adapted to the strong nonlinearity due to the large ratio of degrees of freedom concerned by contact conditions.

In Section 2, the modelling of the contact based on a hybrid formulation is presented. We focus our attention on the modelization of the hinge between two slats characterized by play and eventually by friction (Section 2). The system obtained is nonlinear and nondifferentiable, it is solved by a generalized Newton method (GNM) [4]. In Section 2.5.2, we specify the limits of the GNM associated with classical linear solvers. To overcome these difficulties, we develop, in Section 3 a domain decomposition strategy

<footnote><sup>a</sup> Mechanics and Civil Engineering Laboratory, University of Montpellier, C.C.048, 34095 Montpellier Cedex 5, France</footnote>

<footnote><sup>b</sup> INRIA, Rocquencourt, Bt 16, BP 105, 78153 Le Chesnay Cedex, France</footnote>
for this nonclassical frictional multi-contact problem. This strategy is well adapted to coarse grain parallel computers and take advantage of the multi-processor architectures. In Section 4, the numerical performance of this strategy is illustrated. A parametric study is performed on sequential and parallel computers.

## 2. Modelling of a hinge with clearance and friction

### 2.1. Bi-unilateral contact and frictional rotative laws

In contact continuum theory undergoing finite deformations, the unilateral contact conditions between two deformable solids can be introduced as a tribological interface law rather than as boundary conditions [8,15]. Among the surface interactions, contact and friction are the more usual. But the modelling of a hinge between plates with clearance (Fig. 1) imposes to define specific frictional contact laws.

**Bi-unilateral contact law.** Indeed as shown in Fig. 1, we have to consider a double unilateral contact over the upper and lower faces of the hinge, so the contact law (Fig. 2) is more complicated than the usual case; we can call it a bi-unilateral contact law. For a play equal to $2g$ we have a upper or a lower contact for a relative deflexion $d = \pm g$. So these two unilateral condition relate the so-called gap distance $d$ to the contact shear stress $\tau_e$ (Fig. 2). According to convex analysis [13], this multi-valued relation can be derived from conjugate nondifferentiable convex potentials (in the sense of the subgradient):

$$ \tau_e \in \partial\Psi_G(d) \quad \text{or} \quad d \in \partial\Psi_G^*(\tau_e), \qquad (1) $$

where $\partial\Psi_G$ denotes the subdifferential of the indicator function $\Psi_G$ of the interval $G = [-g, g]$. $\Psi_G^*$ is the Legendre Fenchel form of $\Psi_G$.

**Threshold rotative friction law.** The second particularity is that the friction law (Fig. 3) links the relative rotation increment between two jointed slats to a frictional torque. Moreover the frictional torque depends on the contact stress according to a Coulomb type coefficient. This friction law is defined for a contact

<figure><img src="image_2.png" /><figcaption>Fig. 1. Elastic plates and guide with clearance.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Bi-unilateral contact law.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Frictional rotative law.</figcaption></figure>

shear stress $\tau_e$ which is assumed to be known and is described by three conditions: the friction criterion (Coulomb criterion); the slip rule; a complementary condition. These three relations can be written in terms of the relative angular velocity $\dot{\beta}_t$ (according to the tangent to the hinge) and the friction shear $\tau$ (decomposed in a contact shear stress $\tau_e$ and a frictional torque $\tau_t$). The graph of this multi-valued tribological law, represented in Fig. 3, summarizes the relationship between $\tau_t$ and $\dot{\beta}_t$. According to convex analysis [13], the friction law can be written as follows:

$$ \dot{\beta}_t \in \partial \Psi_{C[\tau_e]}(\tau_t) \quad \text{or} \quad \tau_t \in \partial \Psi_{C[\tau_e]}^*(\dot{\beta}_t), \qquad (2) $$

where $\Psi_{C[\tau_e]}^*(\dot{\beta}_t)$ is the Fenchel conjugate form of the indicator function $\Psi_{C[\tau_e]}$. $C[\tau_e]$ denotes the interval $[-\mu|\tau_e|, \mu|\tau_e|]$.

**Frictional contact law.** A law of frictional bi-unilateral contact can be constructed by combining the bi-unilateral contact law with the frictional rotative law. This (nonassociated) frictional contact law can be written in the form (1) plus (2), provided the dependence of the friction disc $C[\tau_e]$ upon the shear contact stress $\tau_e$ is taken into account,

$$ \begin{aligned} \tau_e &\in \partial \Psi_G(d) \quad \text{or} \quad d \in \partial \Psi_G^*(\tau_e), \\ \tau_t &\in \partial \Psi_{C[\tau_e]}^*(\dot{\beta}_t) \quad \text{or} \quad \dot{\beta}_t \in \partial \Psi_{C[\tau_e]}(\tau_t). \end{aligned} \qquad (3) $$

So this law is nonsymmetric: the contact stress is independent of the friction law but the inverse is false. We can then find the associated frictional contact status of a particle in the hinge:

$$ \begin{aligned} -g < d < g, \quad \tau_e = 0 & \quad \text{gap,} \\ d = +g, \quad \tau_e > 0; & \quad |\dot{\beta}_t| = 0, \quad |\tau_t| \le \mu|\tau_e| & \quad \text{contact+, stick,} \\ d = +g, \quad \tau_e > 0; & \quad |\dot{\beta}_t| > 0, \quad |\tau_t| = \mu|\tau_e| & \quad \text{contact+, slip,} \\ d = -g, \quad \tau_e < 0; & \quad |\dot{\beta}_t| = 0, \quad |\tau_t| \le \mu|\tau_e| & \quad \text{contact-, stick,} \\ d = -g, \quad \tau_e < 0; & \quad |\dot{\beta}_t| > 0, \quad |\tau_t| = \mu|\tau_e| & \quad \text{contact-, slip.} \end{aligned} $$

This tribological interface law can be introduced to complete the formulation of the mechanical behavior of the plates system $\Omega_1 \cup \Omega_2$ composed by the contactor and target plates (Fig. 1). A Kirchhoff Love model is selected for the behavior of the two plates.

## 2.2. Hybrid formulation

We focus our attention on the incremental problem which may be formulated as a minimization problem:

$$ I_{\text{sys}}(\hat{\omega}; \tau_e) := \inf_{\hat{\omega}^* \in V_1 \times V_2} \left\{ I_{e1}(\omega_1^*) + I_{e2}(\omega_2^*) + I_e(\hat{\omega}^*; \tau_e) \right\}, \qquad (4) $$
where $I_{\text{sys}}$ is the total potential energy of this two jointed plates system. $V_1$ and $V_2$ denote the virtual flexion displacement spaces of contactor and target plates. $\hat{\omega}$ (resp. $\hat{\omega}^*$) is the real (resp. virtual) deflexion of the global system $\Omega_1 \cup \Omega_2$ and is defined by $\hat{\omega} = [\omega_1, \omega_2]$ (resp. $\hat{\omega}^* = [\omega_1^*, \omega_2^*]$. $I_{c1}$ and $I_{c2}$ are the sum of the internal elastic energies stored in each plate minus the external potential. $I_c$ represents a quasi-potential induced by the frictional contact effects in the hinge resulting from the jointed plates reciprocal interactions. $I_c$ and $I_c$ are defined by

$$ 
\begin{aligned}
I_c(\omega(x, t)) &= \frac{1}{2} \int_{\Omega} D_{\alpha\beta\gamma\delta} \kappa_{\gamma\delta}(\omega(\mathbf{x})) \kappa_{\alpha\beta}(\omega(\mathbf{x})) d\Omega - \int_{\partial_f \Omega} \mathbf{f} \cdot \mathbf{u}(\mathbf{x}, t) d\partial_f \Omega, \\
I_c(\hat{\omega}^*, \tau_c) &= \int_{\Gamma^c} \Psi_G(d(\hat{\omega}^*)) d\gamma + \int_{\Gamma^c} \Psi_{C(\tau_c(\hat{\omega}))}^*(\hat{\beta}_t(\hat{\omega}^*)) dt d\gamma,
\end{aligned}
\quad (5) $$

where $D$ is the stiffness tensor and $\kappa$ is the tensor of changes of curvature of plates. In the expression of $I_c$, the first integral represents the constraint (from a mathematical point of view) due to the bi-unilateral contact. The second one provides a supplementary energy to the system induced by the rotative friction effects in the hinge. The problem (4) is not a standard minimization problem. Indeed the Coulomb disk is function of the contact shear stress $\tau_c$ which depends on the solution $\hat{\omega}$. For this reason, the minimization problem (4) is considered as a quasi-variational problem. To overcome these difficulties, Alart and Curnier [4] developed a mixed penalty duality formulation of frictional contact problem inspired from an augmented Lagrangian approach to regularize the nondifferentiable contact and friction terms. A quasi-augmented Lagrangian is postulated [1,4,7] for the frictional contact problem; it gives an unconstrained formulation:

$$ L_{\text{sys}}^{r, \rho}(\hat{\omega}, v, \gamma; \tau_c) = \inf_{\hat{\omega}^*} \sup_{v^*} \sup_{\gamma^*} \left\{ I_{c1}(\omega_1^*) + I_{c2}(\omega_2^*) + L_c^{r, \rho}(\hat{\omega}^*, v^*, \gamma^*; \tau_c) \right\}, $$

with

$$ L_c^{r, \rho}(\hat{\omega}^*, v^*, \gamma^*; \tau_c) = \int_{\Gamma^c} l_c^\rho(\hat{\omega}^*, v^*) d\gamma^* + \int_{\Gamma^c} l_t^r(\hat{\omega}^*, \gamma^*; \tau_c) d\gamma^* $$

In the context of contact mechanics, the Lagrange multipliers $v$, $\gamma$ ($v^*$, $\gamma^*$ for the virtual variables) can be interpreted as the nominal contact shear stress and the nominal frictional torque, respectively. The terms $l_c^\rho$ and $l_t^r$ represent the regularization of the functions $\Psi_G(d)$ and $\Psi_{C(\tau_c(\hat{\omega}))}^*(\gamma)$,

$$ 
\begin{aligned}
l_c^\rho(\hat{\omega}^*, v^*) &= -\frac{1}{2\rho}|v^*|^2 + \frac{1}{2\rho}\text{dist}_{\rho G}^2\{v^* + \rho d(\hat{\omega}^*)\}, \\
l_t^r(\hat{\omega}^*, \gamma^*; \tau_c) &= \hat{\beta}_t(\hat{\omega}^*)dt\gamma^* + \frac{r}{2}|\hat{\beta}_t(\hat{\omega}^*)dt|^2 - \frac{1}{2r}\text{dist}_C^2\{\gamma^* + r\hat{\beta}_t(\hat{\omega}^*)dt\},
\end{aligned}
\quad (6) $$

where $\rho$ and $r$ are the positive penalty factors related respectively to bi-unilateral contact and rotative friction. According to the Alart and Curnier approach [4], we introduce a particular form of the “quasi” Lagrangian by substituting an augmented convex set $\hat{C}$ for the original one $C[\tau_c]$:

$$ \hat{C}(v, d(\hat{\omega})) := C(\text{prox}_{\rho G} (v + \rho d(\hat{\omega}))). $$

This approach of “quasi”-augmented Lagrangian permits to satisfy exactly the contact constraints and friction criteria contrary to penalty techniques. The differentiability regions on the augmented Lagrangian are presented in Fig. 4.

A solution $(\hat{\omega}, \lambda)$, $(\lambda = (v, \gamma))$ of the quasi-minimax problem is characterized by vanishing of the first variations of $L_{\text{sys}}^{r, \rho}$

$$ D I_c(\hat{\omega}, \delta\hat{\omega}) + D L_c^{r, \rho}((\hat{\omega}, \lambda; \hat{\omega}, \lambda); (\delta\hat{\omega}, \delta\lambda)) = 0, \quad \forall \delta\hat{\omega}, \forall \delta\lambda, \quad (7) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Augmented frictional contact laws.</figcaption></figure>

where the contact part has the following expression,

$$DL_{e}^{\rho,r}(\hat{\omega}, \lambda; (\delta\hat{\omega}, \delta\lambda)) = \int_{\Gamma^{c}} \left( \frac{\delta\hat{\omega}}{\delta\lambda} \right) \left( \begin{array}{l} \nabla_{\hat{\omega}^{*}} l_{e}^{\rho}(\hat{\omega}, v) + \nabla_{\hat{\omega}^{*}} l_{t}^{r}(\hat{\omega}, \gamma; v) \\ \nabla_{\lambda^{*}} l_{e}^{\rho}(\hat{\omega}, v) + \nabla_{\lambda^{*}} l_{t}^{r}(\hat{\omega}, \gamma; v) \end{array} \right) d\Gamma.$$

The calculation of the gradient gives

$$\begin{align}
\nabla_{\hat{\omega}^*} (l_e^\rho(\hat{\omega}, v) + l_t^r(\hat{\omega}, \gamma)) &= [\nabla_{\hat{\omega}^*} d(\hat{\omega})]^T \text{prox}_{\psi_{\hat{\rho}}^*} (v + \rho d(\hat{\omega})) + [\nabla_{\hat{\omega}^*} \hat{\beta}_t(\hat{\omega})]^T \text{proj}_{\hat{C}}(\gamma + r \hat{\beta}_t(\hat{\omega}) dt), \nonumber \\
\nabla_{\hat{\omega}^*} (l_e^\rho(\hat{\omega}, v) + l_t^r(\hat{\omega}, \gamma)) &= -\frac{1}{\rho} \{v - \text{prox}_{\psi_{\hat{\rho}}^*} (v + \rho d(\hat{\omega}))\} [\nabla_{\hat{\omega}^*} d(\hat{\omega})] - \frac{1}{r} \{\gamma - \text{proj}_{\hat{C}}(\gamma + r \hat{\beta}_t(\hat{\omega}) dt)\}. \tag{8}
\end{align}$$

## 2.3. Hinge contact elements (Fig. 5)

The discretization of frictional contact interface between the two plates is performed by a specific hybrid finite contact element which gives an elementary contribution of the frictional contact stress and an associated elementary tangent matrix (see Eq. (12)). This description requires the use of a GNM to solve the resulting nonlinear equation (see Section 2.4).

Using the finite element method to discretize the solids, the discrete contact interface depends on the discretization and the motions of the bodies. Generally a contact element is composed by

* one node of the contactor surface, numbered 1;

* a local geometry of the target surface composed by nodes (numbered 2, 3, ..., q);

* a “Lagrange multiplier” or fictitious node numbered q + 1 whose degrees of freedom are the frictional contact force components.

<figure><img src="image_6.png" /><figcaption>Fig. 5. “Hinge” contact element.</figcaption></figure>
We denote by $^e\mathbf{u}$ (or $\mathbf{u}$ if no ambiguity arises) the generalized displacements of the $e$th element and $^e\boldsymbol{\lambda}$ (or $\boldsymbol{\lambda}$) the generalized contact forces:

$$ ^e\mathbf{u} = \mathbf{u} = (\mathbf{u}_1, \mathbf{u}_2, \dots, \mathbf{u}_q); \quad ^e\boldsymbol{\lambda} = "^e\mathbf{u}_{q+1}" = \boldsymbol{\lambda}. $$

We have to define a nonclassical contact element to model the discretization of the hinge with clearance and rotative friction. The discretization of the rolling shutters uses elastic finite plates (DKTP elements) which have three nodes with three degrees of freedom for each node, the deflexion $\omega$, the two rotations $\beta_x$ and $\beta_y$ with respect to the axis $x$ and $y$. Similarly to the DKTP elastic finite element, this hinge element involves three nodes (the contactor, the target node and the contact stress node). We denote by $\tilde{\omega}$ the vector concatenating the deflexion $\omega$ and the two rotations $\beta_x$ and $\beta_y$, $\boldsymbol{\lambda}$ the vector concatenating the contact force $\mathbf{v}$ and the two torque components $\gamma_x$ and $\gamma_y$:

$$ \tilde{\omega} = \omega\mathbf{e} + \beta_x\mathbf{t}_x + \beta_y\mathbf{t}_y = \omega\mathbf{e} + \boldsymbol{\beta}, $$

$$ \boldsymbol{\lambda} = \mathbf{v}\mathbf{e} + \gamma_x\mathbf{t}_x + \gamma_y\mathbf{t}_y = \mathbf{v}\mathbf{e} + \boldsymbol{\gamma}, $$

where $\mathbf{e}$ denotes the normal vector to the plane of the plate and $\mathbf{t}$ the tangent vector to the hinge in the plane of the plate (normal to the vector $\mathbf{n}$).

The kinematic contact variables are the gap $d(\tilde{\omega})$ and the relative rotation increment $\delta\beta(\tilde{\omega})$ between two loading increments $i$ and $i+1$,

$$ d(\tilde{\omega}) = \omega_1 - \omega_2, \quad \delta\beta(\tilde{\omega}) = (\beta_1^{i+1} - \beta_1^i) - (\beta_2^{i+1} - \beta_2^i) = \delta\beta_1 - \delta\beta_2. $$

By using the expression of the gradients given by (8), we get the form of the elementary frictional contact operator $\mathcal{F}(\omega_1, \omega_2, \boldsymbol{\lambda})$:

$$ \mathcal{F}(\tilde{\omega}, \boldsymbol{\lambda}) = \left[ \begin{array}{c} [\nabla_{\tilde{\omega}} d]^T \text{prox}_{\psi_{\rho G}^*}(\sigma) + [\nabla_{\tilde{\omega}} \delta \beta_t]^T \text{proj}_{C(\text{prox}_{\psi_{\rho G}^*}(\sigma))}(\tau_t)\mathbf{t} \\ -\frac{1}{f(\rho,r)} \left\{ \boldsymbol{\lambda} - [\nabla_{\tilde{\omega}} d]^T \text{prox}_{\psi_{\rho G}^*}(\sigma) - [\nabla_{\tilde{\omega}} \delta \beta_t]^T \text{proj}_{C(\text{prox}_{\psi_{\rho G}^*}(\sigma))}(\tau_t)\mathbf{t} \right\} \end{array} \right], \quad (9) $$

where

$$ \sigma = v + \rho d(\tilde{\omega}) = \mathbf{e}^T(\boldsymbol{\lambda} + \rho(\omega_1 - \omega_2)), $$

$$ \tau = \gamma + r\delta\boldsymbol{\beta}. $$

$\sigma$ and $\tau$ are the augmented Lagrange multipliers associated, respectively, to the contact force and the rotative frictional torque ($\rho$ and $r$ are the penalty parameters). The function $f(\rho, r)$ is equal to $\rho$ for the deflexion degrees of freedom and to $r$ for the rotations. Moreover with the gradients of the gap and the relative rotation increment ($\nabla_{\tilde{\omega}} g_e(\tilde{\omega}) = (\mathbf{e}^T, -\mathbf{e}^T)$ and $\nabla_{\tilde{\omega}}(\delta\beta_t) = (\mathbf{t}\mathbf{t}^T, -\mathbf{t}\mathbf{t}^T)$), the elementary frictional contact operator has the following expression:

$$ \mathcal{F}(\omega_1, \omega_2, \boldsymbol{\lambda}) = \left[ \begin{array}{c} \mathbf{F}(\omega_1, \omega_2, \boldsymbol{\lambda}) \\ -\mathbf{F}(\omega_1, \omega_2, \boldsymbol{\lambda}) \\ \frac{-1}{f(\rho,r)}(\boldsymbol{\lambda} - \mathbf{F}(\omega_1, \omega_2, \boldsymbol{\lambda})) \end{array} \right] \quad (10) $$

with

$$ \mathbf{F}(\omega_1, \omega_2, \boldsymbol{\lambda}) = \mathbf{e} \text{ prox}_{\psi_{\rho G}^*}(\sigma) + \text{proj}_{C(\text{prox}_{\psi_{\rho G}^*}(\sigma))}(\tau). $$

Seven frictional contact status are associated to this hinge contact interactions (laws (1) and (2)) according to the slat touches the upper ($d = +g$) or lower ($d = -g$) part of the hinge. The contact operator $\mathcal{F}(\omega_1, \omega_2, \boldsymbol{\lambda})$ corresponding to these contact status, takes the following form:
$$
F(\omega_1, \omega_2, \boldsymbol{\lambda}) =
\begin{cases}
\mathbf{0} & |\sigma| \le \rho g \quad \text{gap,} \\
(\sigma - \rho g)\mathbf{e} + \boldsymbol{\tau} & \begin{cases} \sigma > \rho g \\ \|\boldsymbol{\tau}\| < \mu(\sigma - \rho g) \end{cases} \quad \text{upper stick,} \\
(\sigma - \rho g)(\mathbf{e} - \varepsilon\mu\mathbf{t}) & \begin{cases} \sigma > \rho g \\ \|\boldsymbol{\tau}\| \ge \mu(\sigma - \rho g) \end{cases} \quad \text{upper } \varepsilon \text{ slip,} \\
(\sigma + \rho g)\mathbf{e} + \boldsymbol{\tau} & \begin{cases} \sigma < -\rho g \\ \|\boldsymbol{\tau}\| < -\mu(\sigma + \rho g) \end{cases} \quad \text{lower stick,} \\
(\sigma + \rho g)(\mathbf{e} - \varepsilon\mu\mathbf{t}) & \begin{cases} \sigma < -\rho g \\ \|\boldsymbol{\tau}\| \ge -\mu(\sigma + \rho g) \end{cases} \quad \text{lower } \varepsilon \text{ slip,}
\end{cases}
\tag{11}
$$

($\epsilon$ equals to +1 or -1 with respect to the slip direction (+ or -)). The tangent matrix $\mathbf{A}_{\text{ec}} = \nabla_{\omega_1,\omega_2,\boldsymbol{\lambda}} \mathcal{F}(\omega_1, \omega_2, \boldsymbol{\lambda})$ associated to this operator has the following form:

$$
\mathbf{A}_{\text{ec}} = \nabla_{\omega_1, \omega_2, \boldsymbol{\lambda}} \mathcal{F}(\omega_1, \omega_2, \boldsymbol{\lambda}) = \begin{bmatrix} f(\rho, r)\mathbf{M} & -f(\rho, r)\mathbf{M} & \mathbf{M} \\ -f(\rho, r)\mathbf{M} & f(\rho, r)\mathbf{M} & -\mathbf{M} \\ \mathbf{M} & -\mathbf{M} & \frac{1}{f(\rho, r)}(\mathbf{M} - \mathbf{I}) \end{bmatrix} \quad (12)
$$

with

$$
\mathbf{M} = \begin{cases} \mathbf{0} & \text{gap status,} \\ \mathbf{I} & \text{stick status (for the play } g \text{ or } -g\text{),} \\ (\mathbf{e} - \varepsilon\mu\mathbf{t})\mathbf{e}^{\text{T}} & \text{slip status (for the play } g \text{ or } -g\text{).} \end{cases} \tag{13}
$$

2.4. Generalized Newton method (GNM)

In the augmented Lagrangian approach, the equilibrium of a discretized contact bodies system is
characterized by the first variation of L<sub>sys</sub><sup>r</sup> (7) which takes the form,

$$
\mathbf{G}(\mathbf{x}) + \mathcal{F}(\mathbf{x}) = 0. \tag{14}
$$

So we distinguish two parts involving the pair **x** = (**u**, **λ**), a differentiable part **G**(**x**)(= **G**(**x**)) and a non-
differentiable one **F**(**x**). **G** represents the elastic part and **F** denotes the global frictional contact operator
which is obtained by assembling the elementary contributions (11). We treat both variables simultaneously
through Newton’s method. To overcome the nondifferentiability of Eq. (14), Newton’s method may be
extended to the following iterative form [4,6]:

$$
\begin{align}
\mathbf{x}^{i+1} &= \mathbf{x}^i - (\mathbf{K}^i + \mathbf{A}_c^i)^{-1} (\mathbf{G}(\mathbf{x}^i) + \mathcal{F}(\mathbf{x}^i)), \tag{15} \\
\mathbf{K}^i &= \partial \mathbf{G}(\mathbf{x}^i), \quad \mathbf{A}_c^i \in \partial \mathcal{F}(\mathbf{x}^i). \nonumber
\end{align}
$$

The matrix **K**<sup>*i*</sup> is the usual elastic stiffness matrix and **A**<sub>c</sub><sup>*i*</sup> represents a tangent matrix issued from the as-
sembly of the tangent matrix **A**<sub>ec</sub> (defined in (12) for example). For a given local status, the generalized
Jacobian is reduced to the single classical Jacobian matrix. The GNM leads us to solve at each iteration *i*
the following linear system:

$$
(\mathbf{K}^i + \mathbf{A}_c^i) \Delta \mathbf{x}^i = -(\mathbf{G}(\mathbf{x}^i) + \mathcal{F}(\mathbf{x}^i)), \quad \text{where} \quad \Delta \mathbf{x}^i = \mathbf{x}^{i+1} - \mathbf{x}^i. \tag{16}
$$

The linear systems may be solved by an iterative method as the square preconditioner conjugate gradient method (denoted by SPCG) [5].

2.5. A benchmark: the rolling shutters

2.5.1. The rolling shutters

The aim of the study is to simulate the quasi-static behavior of such shutters submitted to strong winds (the
ultimate aim is to study the dynamic behavior but is beyond the aim of this paper). A rolling shutter is a
<figure><img src="image_7.png" /><figcaption>Fig. 6. Rolling shutters composed by hinged slats.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Rolling shutters composed by 16 jointed slats.</figcaption></figure>

particular specific case of multi-contact structure. The rolling shutters for shops, stores and hangars are formed by a succession of slats jointed by a hinge (Fig. 6). Such a structure is then composed by an assembly of elastic structures (plates in flexion and torsion) which leads to consider a large number of contact zones. The edges of the slates are designed in such a way that the slats fit into each other. To facilitate the rolling of the shutters at the opening, the profile of the slat requires a clearance or a play in the hinge. So we use the the modelization developed in Section 2 to model hinges with clearance and eventually rotative friction between slats.

### 2.5.2. The limits of the classical numerical treatment

The example chosen (Fig. 7) concerns a square shutter of 4 m length, 5 mm of thickness with a play of 3 mm embedded on its boundary and submitted to its own weight. To analyse the influence of our frictional contact modelization we compare two models of hinge: a hinge with play and friction (Fig. 8); a perfect hinge without play and friction (Fig. 9). We notice that the maximal deflexion passes from 15.24 cm for the classical hinge to 17.02 cm for our specific hinge model. For more mechanical considerations (study of the influence of the friction or the play) the reader can refer to [1,7]. To carry out the mechanical study, we have met some difficulties due to the multiplicity of the contact interfaces. In Table 1 we show the behavior of the classical algorithm Newton-SPCG according to the number of slats. The multi-contact structure composed by the rolling shutters leads to drastically nonlinear large-scale ill-conditioned problem with a large ratio of contact degrees of freedom. The contact ratio increases with the slats number (Table 1). Simultaneously, the number of Newton iterations increases moderately. In addition the CPU time blows up and the square conjugate gradient method does not converge for the shutter composed by 16 and 32 slats.

In conclusion, a high contact ratio spoils the behavior of classical methods and we need specific strategies well suited to the solution of large scale problems with a large number of contact zones. Moreover to overcome the difficulties of classical sequential solutions the domain decomposition methods are a good alternative. In the next section, we present a nonlinear domain decomposition strategy well suited for large scale nonlinear nonsymmetric and ill-conditioned problems.
<figure><img src="image_9.png" /><figcaption>Fig. 8. Deformation with hinge (with play and friction, $\mu = 0.2$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Deformation with perfect hinge (without play and friction).</figcaption></figure>

<table><caption>Table 1
Numerical behavior of the classical algorithm according to the contact ratio</caption><thead><tr><th rowspan="2">Description of the rolling shutter</th><th colspan="4">Number of slats</th></tr><tr><th>4</th><th>8</th><th>16</th><th>32</th></tr></thead><tbody><tr><td>Total dof</td><td>6924</td><td>7704</td><td>9294</td><td>12 384</td></tr><tr><td>Contact element</td><td>96</td><td>224</td><td>480</td><td>992</td></tr><tr><td>Contact ratio</td><td>12.4%</td><td>26.1%</td><td>46.5%</td><td>72.1%</td></tr><tr><td>Newton ite.</td><td>9</td><td>12</td><td>?</td><td>?</td></tr><tr><td>CPU time (s)</td><td>1824</td><td>11 534</td><td>No conv.</td><td>No conv.</td></tr></tbody></table>

### 3. Nonlinear domain decomposition strategies for frictional multi-contact problems

Domain decomposition methods (substructuring techniques) are efficient because they allow to reduce memory storage and calculation time. Moreover these methods take advantage of the new multi-processor generation of computers as they exhibit an intrinsic parallelism with a high granularity. The main component of the domain decomposition algorithm is a numerical solver based on the solution of local independent subproblems on subdomains. In addition, these methods are efficient solvers in a classical mono-processor environment as well. First we investigate domain decomposition methods introduced for linear systems. The method used hereafter is the primal Schur complement method which consists in imposing the displacement continuity on the interfaces and in controlling the normal stress gap. Next we develop a strategy to solve a large scale multi-contact problem by using the GNM combined to the domain decomposition method.
### 3.1. Schur complement method in elasticity

#### 3.1.1. The principle: reduction to an interface problem

In substructuring methods, the parallel solution of a structural problem is achieved by splitting the original domain of computation $\Omega$ into smaller subdomains $\Omega^n$ ($n = 1, N$) and by reducing the initial problem to an interface system. So we recall the notations and the formulation. We consider a system of linear algebraic equations,

$$ \mathbf{Ku} = \mathbf{F} \tag{17} $$

with $\mathbf{K}$ a square, symmetric positive definite matrix arising from finite element discretization of a linear, elliptic, self-adjoint boundary value problem on a domain $\Omega$.

The first step of the Schur complement method consists then in splitting the domain into small local nonoverlapping subdomains with interfaces defined as follows (see Fig. 10):

$$ \Gamma^n = \partial\Omega^n \cap \left( \bigcup_{p=1, p \neq n}^{N} \partial\Omega^p \right) - \partial\Omega. $$

With this partition we denote by $\mathbf{K}^n$ and $\mathbf{F}^n$ the subdomain stiffness matrix and right hand side. For each subdomain we can distinguish the internal degrees of freedom (denoted by $i$) from the interface degrees of freedom (denoted by $l$). So the expressions of $\mathbf{K}^n$, $\mathbf{u}^n$ and $\mathbf{F}^n$ take this form:

$$ \mathbf{K}^n = \begin{bmatrix} \mathbf{K}_{ii}^n & \mathbf{K}_{il}^n \\ \mathbf{K}_{li}^n & \mathbf{K}_{ll}^n \end{bmatrix}, \quad \mathbf{u}^n = \begin{Bmatrix} \mathbf{u}_i^n \\ \mathbf{u}_l^n \end{Bmatrix}, \quad \mathbf{F}^n = \begin{Bmatrix} \mathbf{F}_i^n \\ \mathbf{F}_l^n \end{Bmatrix}. \tag{18} $$

The second step consists in reducing the global system to the interface problem by a block Gaussian elimination of the internal degrees of freedom. Finally, the interface problem to be solved is written as:

$$ \text{find } \boldsymbol{\mu} \text{ such that } S\boldsymbol{\mu} = \mathbf{B}, \tag{19} $$

where

$$ \begin{align} \mathbf{S} &= \sum_{n=1}^{N} (\mathbf{R}^n)^t \mathbf{S}^n \mathbf{R}^n, \quad \mathbf{S}^n = \mathbf{K}_{ll}^n - \mathbf{K}_{li}^n (\mathbf{K}_{ii}^n)^{-1} \mathbf{K}_{il}^n, \notag \\ \mathbf{B} &= \sum_{n=1}^{N} (\mathbf{R}^n)^t (\mathbf{F}_i^n - \mathbf{K}_{li}^n (\mathbf{K}_{ii}^n)^{-1} \mathbf{F}_i^n), \quad \mathbf{B}^n = \mathbf{F}_i^n - \mathbf{K}_{li}^n (\mathbf{K}_{ii}^n)^{-1} \mathbf{F}_i^n. \tag{20} \end{align} $$

$\mathbf{R}^n$ is the restriction operator which goes from $\Gamma$ to $\Gamma^n$. The matrix $\mathbf{S}$ is the Schur complement matrix; $\mathbf{S}^n$ are the local Schur complement matrices.

#### 3.1.2. Interface problem solution

Modern domain decomposition methods solve the interface problem by an iterative preconditioned conjugate gradient type algorithm. The desirable preconditioner has nice parallel properties and its per-

<figure><img src="image_11.png" /><figcaption>Fig. 10. Substructuraction of the domain $\Omega$ in $N$ subdomains.</figcaption></figure>
formance must be insensitive to the discretization step and the number of subdomains. Notice that even without preconditioning the interface system has a better condition number (cond($\mathbf{S}$) = $\mathcal{O}((1/H^2)(1+H/h))$, where $H$ is the subdomain size and $h$ is the mesh size) than the initial system (cond($\mathbf{K}$) = $\mathcal{O}(1/h^2)$). Hereafter, we use the multi-level Neumann Neumann preconditioner. This iterative technique never requires the explicit calculation of the matrix $\mathbf{S}$. We have just to form the matrix vector products $\mathbf{S}\mathbf{p}_k$ and $\mathbf{M}\mathbf{r}_{k+1}$ by solving independent auxiliary Dirichlet and Neumann problems on the local subdomains. We will briefly present the multi-level Neumann Neumann preconditioners: the original one and the balancing method introduced by Mandel [12].

*Original Neumann Neumann preconditioner.* The inverse of the sum $\mathbf{S} = \sum_{n=1}^{N} (\mathbf{R}^n)^t \mathbf{S}^n \mathbf{R}^n$ is approximated by a weighted sum of the matrices $(\mathbf{S}^n)^{-1}$,

$$ \mathbf{M} = \sum_{n=1}^{N} (\mathbf{R}^n)^t \mathbf{P}^n (\bar{\mathbf{S}}^n)^{-1} \mathbf{P}^n \mathbf{R}^n. \qquad (21) $$

$\mathbf{P}^n$ is a diagonal weighting matrix which must verify $\sum_{n=1}^{N} (\mathbf{R}^n)^t \mathbf{P}^n \mathbf{R}^n = \text{Id}|_R$. A generic choice is: $\mathbf{P}^n = \mathbf{K}_{ii}^n / \sum_{n=1}^{nd_i} \mathbf{K}_{ii}^n$ with $nd_i$ the subdomains which contain the degrees of freedom $i$. When the local Neumann problems are not well posed (no Dirichlet boundary conditions), $(\bar{\mathbf{S}}^n)^{-1}$ denotes an arbitrary regularized inverse. When this preconditioner is used the condition number is $(\text{cond}(\mathbf{MS}) = \mathcal{O}((1/H^2)(1 + [\log(H/h)]^2)))$. This preconditioner is not optimal as it does not scale well with the number of subdomains.

*Balancing method.* Mandel [12] has generalized the Neumann Neumann preconditioner by accounting for the rigid body motions of the floating subdomains. The idea of Mandel is to add to the local independent subproblems on subdomains a global coarse problem with few unknowns for each subdomains. The product of the preconditioner $\mathbf{M}$ and of the residual gradient $\mathbf{r}_{k+1}$ have the following form:

$$ \mathbf{M}\mathbf{r}_{k+1} = \sum_{n=1}^{N} \left\{ (\mathbf{R}^n)^t \mathbf{P}^n (\bar{\mathbf{S}}^n)^{-1} \mathbf{P}^n \mathbf{R}^n \mathbf{r}_{k+1} \right\} + \mathbf{G}_{\gamma}, $$

where

$$ (\bar{\mathbf{S}}^n)^{-1} = \begin{cases} (\mathbf{S}^n)^{-1} & \text{if } \Omega^n \text{ is not a floating subdomain,} \\ \text{the inverse of the projection of the image of } \mathbf{S}^n, & \end{cases} $$

and $\mathbf{G}_{\gamma}$ is the projection of the set of the linear combination of the rigid body motion of the subdomains over the interface. In practice, we solve a global optimization problem [9,10] over the interface $\Gamma$ in order to minimize the residual:

$$ \min_{\mathbf{G}_{\gamma}} \|\mathbf{r}_{k+1}\|_{\Gamma}^{2} := \min_{\mathbf{G}_{\gamma}} \{(\mathbf{S}(\mathbf{M} - \mathbf{S}^{-1})\mathbf{r}_{k+1})^{t}(\mathbf{M} - \mathbf{S}^{-1})\mathbf{r}_{k+1}\}. \qquad (22) $$

Its minimum is attained for the function $\mathbf{G}_{\gamma}$ which vanishes its gradient; so we have to solve the following coarse problem, where $\gamma$ is the unknown,

$$ (\mathbf{G}'\mathbf{S}\mathbf{G})\gamma = -\mathbf{G}'\mathbf{S}\sum_{n=1}^{N} (\mathbf{R}^n)^t \mathbf{P}^n (\bar{\mathbf{S}}^n)^{-1} \mathbf{P}^n \mathbf{R}^n \mathbf{r}_{k+1}. \qquad (23) $$

For the balancing method, the condition number is asymptotically optimal $(\text{cond}(\mathbf{MS}) = \mathcal{O}(1 + (\log(H/h))^2))$; the preconditioner is then insensitive to the number of subdomains.

This preconditioner can also be described and analyzed in a unique abstract framework. Indeed, the Neumann Neumann algorithm is a standard additive Schwarz algorithm [11,14].

## 3.2. Domain decomposition method, Newton solver and frictional contact

The strategy to solve the nonlinear problem is the following: consider the GNM as the nonlinear standard solver and at each Newton iteration solve the linearized problem by an iterative Schur
complement solver. We call this algorithm “Newton Schur” (Section 3.2.1). The main difference compared to a standard Newton procedure consists of specific adaptations to take into account new rigid body motions in floating subdomains with contact zones and the nonsymmetry of the coarse problem (Section 3.2.2).

### 3.2.1. Newton Schur algorithm: general purpose

**The algorithm.** To extend domain decomposition method to frictional contact problem we have to overcome many difficulties: the nonlinearity, the non differentialibility and the nonsymmetry. The non-linearity and the nondifferentiability are treated by the generalized algorithm as for classical methods (Newton-PSGC). The Schur complement method is applied to the linearized system at each Newton iteration. The nonsymmetry led us to choose the GMRes method to solve the interface problem. Finally the multi-level Neumann–Neumann preconditioner (with and without coarse spaces) is carried to overcome the poor conditioning and the treatment of floating subdomains; but we will see in the Section 3.2.2 that the balancing method is sensitive to the nonsymmetry and requires specific adaptations. The Newton Schur algorithm takes the form described in Table 2.

**Subtracturation strategies.** The main feature of this nonlinear domain decomposition strategy consists in distinguishing the physical contact interfaces from the numerical subdomain interfaces. Contrary to current approaches we suggest to treat the physical contact interfaces inside the subdomains. For the rolling shutters, the contact interfaces (here the hinges) must be inside the subdomains and do not constitute interfaces as presented in the Fig. 11 for a band and square decomposition. The decomposition of the mechanical system is not forced to respect the geometry of its components. It makes it possible to balance the size of the subdomains and get an optimal decomposition for parallel efficiency.

### 3.2.2. Newton Schur algorithm: specific adaptations

The balancing method requires two specific adaptations due to frictional contact conditions

#### 3.2.2.1. Treatment of the rigid body motions.
The convergence of the proposed domain decomposition method is optimal for the choice of the preconditioner described in Section 3.1.2. At each iteration of the algorithm a local problem involving the matrix $(\tilde{S}^n)^{-1}$ is solved by subdomain. As already mentioned, $(\tilde{S}^n)^{-1}$ denotes an arbitrary regularized inverse. The actual choice, is to consider the exact inverse inverse for nonfloating subdomains. For floating subdomains the construction of a regularized inverse involves the use of a set of elements which define the null space. This choice has an influence on the size of the coarse

<table><caption>Table 2
Newton Schur algorithm for frictional contact problems</caption><tr><td><b>Newton iterations:</b></td><td>i = 1, ...</td></tr><tr><td><b>Schur complement method:</b></td><td><div style="text-align: center;"><span style="display: inline-block; border: 1px solid black; padding: 0 2px; line-height: 1em;">Linearization of the global problem:</span><br><b>G</b>(<b>u</b>) + <i>F</i>(<b>u</b>, <b>λ</b>) = 0 &rarr; <b>K</b><sub>i</sub>Δ<b>x</b> = <b>b</b><sub>i</sub></div></td></tr><tr><td><b>Reduction to a nonsymmetric Γ-interface problem:</b></td><td><div style="text-align: center;"><span style="display: inline-block; border: 1px solid black; padding: 0 2px; line-height: 1em;">K</b><sub>i</sub>Δ<b>x</b> = <b>b</b><sub>i</sub> &rarr; <b>S</b><sub>i</sub><i>χ</i> = <b>B</b><sub>i</sub></div></td></tr><tr><td><b>Interface solution of <i>S</i><sub>i</sub><i>χ</i> = <b>B</b><sub>i</sub> by a GMRes algorithm with multilevel Neumann–Neumann preconditioner</b></td><td>until convergence.</td></tr></table>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Band and square decomposition substructuraction strategies.</figcaption></figure>

problem, and therefore on the cost of the algorithm. The use of rigid body motions allows to minimize the size of the coarse problem and thus improve the numerical efficiency of the algorithm.

Since the contact interfaces are inside the subdomains we have to consider more rigid body motions than in the classical case. The actual number of more rigid body motions depends on the global status of the contact interfaces. For instance, with a square decomposition a floating subdomain with a single contact interface may have six rigid body motions if the contact status in all contact elements is gap, four rigid body motions if this status is slip and three rigid body motions if it is stick (usual case). The explanation of the number of rigid body motions according to the contact conditions is shown in Fig. 12.

### 3.2.2.2. An adapted balancing preconditioner.
The numerical experiments [7] show that the behavior of the iterative Schur a complement solver (GMRes algorithm) is strongly perturbed when friction occurs, i.e., when nonsymmetry is introduced in the tangent matrices. Indeed the minimization problem (25) is not well defined for nonsymmetric problems. The first idea is to replace the matrix S by the symmetrized matrix S<sup>s</sup> (S<sup>s</sup> = S + S<sup>t</sup>). But, as we will see later, a better choice is to use asymmetric matrix which has a mechanical meaning: the Schur complement matrix associated to a frictionless contact status by imposing μ = 0. We consider the interface reduced matrix S<sup>*</sup> with a zero friction coefficient (S<sup>*</sup> = S<sub>μ=0</sub>) to evaluate the norm of the difference between M and S<sup>-1</sup>. Then the minimization problem takes the following form:

$$ \min_{G_{\gamma}} |||\mathbf{r}_{k+1}|||_{I}^{2} := \min_{G_{\gamma}} \{(\mathbf{S}^{*}(\mathbf{M} - \mathbf{S}^{-1})\mathbf{r}_{k+1})^{t}(\mathbf{M} - \mathbf{S}^{-1})\mathbf{r}_{k+1}\}. \quad (24) $$

This minimum is reached for the function G<sub>γ</sub> which vanishes its gradient,

$$ \nabla_{G_{\gamma}} \left\{ 2 \left\langle \sum_{n=1}^{N} \mathbf{H}_{k+1}^{n} \cdot \mathbf{S}^{*} \mathbf{G}_{\gamma} \right\rangle + \langle \mathbf{S}^{*} \mathbf{G}_{\gamma}, \mathbf{G}_{\gamma} \rangle \right\} = 0 $$

with

$$ \mathbf{H}_{k+1}^{n} = ( (\mathbf{R}^{n})^{t} \mathbf{P}^{n} (\bar{\mathbf{S}}^{n})^{-1} \mathbf{P}^{n} \mathbf{R}^{n} ) \mathbf{r}_{k+1}. $$

The solution of this equation verifies the following equality:

$$ (\mathbf{G}^{t} \mathbf{S}^{*} \mathbf{G})_{\gamma} = - \mathbf{G}^{t} \mathbf{S}^{*} \sum_{n=1}^{N} (\mathbf{R}^{n} \mathbf{P}^{n} (\bar{\mathbf{S}}^{n})^{-1} \mathbf{P}^{n} \mathbf{R}^{n}) \mathbf{r}_{k+1}, \quad (25) $$

which defines the coarse problem adapted to the nonsymmetric case.
<figure><img src="image_13.png" /><figcaption>Fig. 12. Detection of rigid body motion for a subdomain with a contact interface.</figcaption></figure>

# 4. Numerical results

In this section, we study the numerical and parallel behavior of the Newton Schur strategy. In the first section we present a parametric study to analyse the convergence behavior of this algorithm and specific adaptations with respect to the decomposition, the number of hinges with and without friction and the number of subdomains. In a second part, we analyse the efficiency of the parallel implementation of our algorithm.

## 4.1. Parametric study

We consider two decompositions for the rolling shutters: a band decomposition without floating subdomains and a square decomposition with floating subdomains which permits to validate our specific adaptations.

### 4.1.1. Band decomposition

The influence of three parameters is studied: the number of hinges, the friction and the number of subdomains.

*Influence of the number of hinges.* The influence on an increasing number of degrees of freedom in the contact zones due to the multiplicity of contact zones on classical solvers is the main motivation of the domain decomposition development (see Section 2.5.2). Table 3 summarizes mesh parameters associated to different rolling shutters with a varying number of slats (i.e. hinges). Fig. 13 shows that for a decomposition
<table><caption>Table 3
Discretization of the rolling shutter according to the number of slats</caption><thead><tr><th rowspan="2">Description of the rolling shutter</th><th colspan="4">Number of slats</th></tr><tr><th>4</th><th>8</th><th>16</th><th>32</th></tr></thead><tbody><tr><td>Total dof</td><td>26 124</td><td>27 672</td><td>30 768</td><td>36 960</td></tr><tr><td>Contact element</td><td>192</td><td>448</td><td>960</td><td>1984</td></tr><tr><td>Contact dof</td><td>1728</td><td>4032</td><td>8640</td><td>17 856</td></tr><tr><td>Contact ratio</td><td>6.6%</td><td>14.5%</td><td>28.1%</td><td>48.3%</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Newton iterations and GMRes average iterations according to the number of slats.</figcaption></figure>

into 31 subdomains, the number of Newton iterations increases moderately but, contrary to the classical method, the linear solver has a good behavior. Indeed the average number of GMRes iterations decreases slowly when the contact ratio increases.

*Influence of the friction.* On an example with 16 slats and 30 768 degrees of freedom, the convergence of the Newton Schur algorithm is not perturbed by the friction: 22 Newton iterations (resp. 52 GMRes average iterations) without friction and 23 (resp. 56 GMRes average iterations) with friction ($\mu = 0.2$). On the other hand memory requirements increases strongly when friction occurs due to the nonsymmetry (from 21.5 Mo without friction, to 37.5 Mo with friction). This enlargement becomes catastrophic for direct solvers (from 104 to 209 Mo) and the algorithm fails for lack of memory.

*Influence of the number of subdomains.* In Table 4, we consider a rolling shutter with 16 slats decomposed in 3, 7, 15 and 31 subdomains. We compare our strategy with the classical algorithm in terms of CPU time. For the example with 9264 degrees of freedom, the CPU time decreases with respect to the number of subdomains until 15 which seems to be an optimal number in this case. With this optimal decomposition, when 9264 degrees of freedom are considered the Newton Schur strategy converges 85 times better than the classical method. When the number of subdomains increases the band decomposition does not respect the requirement of nice aspect ratio for one subdomain. This explains why 15 subdomains are here an optimum.

<table><caption>Table 4
Newton Schur CPU time (and number of GMRes average iterations) according to the number of subdomains comparison with the classical Newton direct method</caption><thead><tr><th rowspan="2">Shutter<br>16 slats</th><th colspan="4">Number of subdomains</th><th rowspan="2">Classical<br>algorithm</th></tr><tr><th>3</th><th>7</th><th>15</th><th>31</th></tr></thead><tbody><tr><td>9264 dof</td><td>2788 (12)</td><td>814 (22)</td><td>520 (58)</td><td>—</td><td>44 265</td></tr><tr><td>30 768 dof</td><td>*</td><td>15 752 (28)</td><td>8050 (70)</td><td>12 281 (188)</td><td>*</td></tr></tbody></table>

<footnote>* Not enough memory.</footnote>
### 4.1.2. Square decomposition

We now carry out a square decomposition of the rolling shutters in 15 × 2, 15 × 3 and 15 × 4 subdomains. Two parameters have to be studied: the friction coefficient and the number of subdomains.

*Influence of the friction coefficient.* Table 5 underlines the contribution of the balancing method when floating subdomains appear in the decomposition. Without friction and without balancing method, the behavior of the GMRes algorithm is spoiled (307 iterations versus 108). By introducing the standard balancing method, the convergence rate is improved (46 iterations versus 73), but the friction cancels this advantage and the standard balancing method magnifies the degradation of the algorithm due to floating subdomains. In order to understand the influence of the nonsymmetry on the algorithm and the advantages of the specific improvements presented in Section 3.2.2, we study friction coefficient varying from 0 to 2 for a given example with 16 slats and 15 × 2 subdomains (26 floating subdomains) (Fig. 14). If the friction coefficient is zero, the problem is symmetric and the two balancing method are equivalent and efficient. As soon as we introduce a small friction coefficient, the number of GMRes iterations increases quickly because nearly all contact elements have a slip status which introduce the nonsymmetry in the matrices. For $\mu = 0.2$, we have the higher ratio of slip status. For $\mu > 0.2$, the ratio of sticking status increases and the nonsymmetry decreases. For $\mu > 2$, we have only stick status and we recover a symmetric problem. The algorithm with the two methods behaves like the evolution of the ratio of slip status, even if the friction coefficient is weak. Notice that the specific balancing method is less sensitive to the nonsymmetry: it requires twice less iterations than the standard one.

*Influence of the number of subdomains.* The good behavior of the specific preconditioner is confirmed when the number of floating subdomains increases (Fig. 15). Beyond 30 subdomains, the adapted balancing method does not depend on the number of subdomains; it behaves almost like the standard method on a frictionless problem. As a first conclusion, the algorithm has the desired numerical scalability. Thus the strategy developed above is efficient for large scale frictional contact problems.

### 4.2. Parallel behavior

In this section, we investigate the parallel scalability of the actual implementation of the Newton Schur strategy. A nonlinear, nonsymmetric problem (the rolling shutters with hinges contact interfaces) is com-

<table><caption>Table 5
Influence of the friction on the GMRes average iterations with multi-level Neumann-Neumann preconditioners ($\mu = 0.2$)</caption><thead><tr><th>16 slats and 9264 dof</th><th colspan="2">Original Neumann</th><th colspan="2">Standard balancing method</th></tr><tr><th>Partitionning</th><th>15×2</th><th>15×4</th><th>15×2</th><th>15×4</th></tr></thead><tbody><tr><td>Hinge without friction</td><td>108</td><td>307</td><td>73</td><td>46</td></tr><tr><td>Hinge with friction</td><td>111</td><td>318</td><td>241</td><td>483</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Influence of the friction coefficient on preconditioned GMRes convergence behavior for a decomposition of 15 × 2 subdomains.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15. Influence of the number of subdomains on preconditioned GMRes convergence behavior.</figcaption></figure>

pared with a linear symmetric example (a rolling shutter hinged without play and friction). For each experiment, we map one subdomain on each processor of the parallel computer (an IBM SP2 with 207 nodes). We use MPI as message passing library. Moreover for each experiment, we try to appreciate the gain of the parallel solution with respect to the sequential one, consequently the problems presented are not very large. In Table 6, the results with the rolling shutters hinged without play and friction (linear symmetric problem) are summarized. The rolling shutter is decomposed in 8, 16 and 32 subdomains, then the parallelization is carried out with 8, 16 and 32 processors. The gain is here defined as the ratio between sequential CPU time and parallel CPU time. For this test the best gain is obtained for a decomposition into 16 subdomains and decreases sensitively for 32 subdomains. For this medium size problem, the poor performance is due to the large size (Table 7) of the interface with respect to the subdomain size. The results of the parallel implementation for a realist example with a rolling shutter with around 400 000 degrees of freedom are presented in Table 6. On the other hand, a similar analysis with the problem of the rolling shutter with play and friction (nonlinear problem) is presented in Table 8. The rolling shutter is decomposed into 7 and 15 subdomains. Now we can compare performances for the nonlinear case with respect to the linear one. The gains for the nonlinear problem are better than for the linear problem. These gains get closer to the maximal gain which is the number of processors (5.6 for 7 processors and 11.8 for 15 processors).

Notice that the problem of how to measure the parallel performance is rather complex. It is obvious that the initial decomposition of the domain has to be well balanced. This may change during the computations, depending on the status of contact zones. In addition, if the number of subdomains changes the whole algorithm is different; when the number of subdomains increases the size of local problems is reduced and the size of the coarse problem grows. An optimal has to be found such that the granularity of local

<table><caption>Table 6
CPU time and gains of the Schur primal complement method</caption><thead><tr><th>Processors used</th><th>8</th><th>16</th><th>32</th></tr></thead><tbody><tr><td>99 075 dof sequential solution (T<sup>1</sup><sub>seq</sub>)</td><td>316.3</td><td>290.2</td><td>429</td></tr><tr><td>99 075 dof parallel solution (T<sup>1</sup><sub>para</sub>)</td><td>58.6</td><td>47.1</td><td>126.1</td></tr><tr><td>Gain = T<sup>1</sup><sub>seq</sub>/T<sub>para.ex1</sub></td><td>5.4</td><td>6.2</td><td>3.4</td></tr><tr><td>394 755 dof parallel solution (T<sup>2</sup><sub>para</sub>)</td><td>528.4</td><td>298.7</td><td>356.7</td></tr></tbody></table>

<table><caption>Table 7
Sizes of the subdomains and the decomposition interfaces</caption><thead><tr><th>Nb. of subdomains</th><th>8</th><th>16</th><th>32</th></tr></thead><tbody><tr><td>Subdomains dof</td><td>12 723</td><td>6555</td><td>3571</td></tr><tr><td>Interface dof</td><td>2709</td><td>5805</td><td>11 997</td></tr></tbody></table>
<table><caption>Table 8
CPU time and gains of the Newton Schur algorithm</caption><thead><tr><th rowspan="2">Rolling shutter with 16 slats</th><th colspan="2">Number of processors/subdomains</th></tr><tr><th>7</th><th>15</th></tr></thead><tbody><tr><td>30 768 dof solution sequential (T<sup>1</sup><sub>seq</sub>)</td><td>15 752</td><td>8050</td></tr><tr><td>30 768 dof solution parallel (T<sup>1</sup><sub>para</sub>)</td><td>2810</td><td>680</td></tr><tr><td>Gain = T<sup>1</sup><sub>seq</sub>/T<sup>1</sup><sub>para</sub></td><td>5.6</td><td>11.8</td></tr><tr><td>110 595 dof solution parallel (T<sup>2</sup><sub>para</sub>)</td><td>*</td><td>7680</td></tr></tbody></table>

problems is preserved and that the time spend to solve the coarse problem remains reasonable compared to the one used for local problems. On the other hand the definition used for the gain (ratio between sequential and parallel time) is not very appropriate because for realistic target problems it is not possible to run on one processor. Nevertheless the results presented show that the use of domain decomposition makes it possible to solve very challenging problems.

## 5. Conclusions

The method presented in this paper and the numerical results which illustrate its potential show that it is now possible to solve real life large contact problems using a domain decomposition approach. This algorithm has nice scalability properties for both numerical and parallel point of view. In addition it is more robust than standard iterative solvers and is reliable on parallel computers. Further investigations are needed to improve the preconditioners which takes into account the fact that the operator is nonsymmetric. Such a preconditioner was studied for advection diffusion problems; a Robin Robin preconditioner [3] which considers the non symmetry of the problem was derived and the results obtained are very promising.

The present approach has been tested on an original problem, the rolling shutters. But is was used on more classical multi-contact problems as granular media with deformable grains or cellular materials involving self contact between the walls of the cells. The advantages of the rolling shutters are the periodicity of the structure and the simplicity of the contact geometry allowing an easy well balanced decomposition. On the other hand, the granular and cellular media present specific difficulties. For cellular material, the large deformations have to be accounted for. For cellular material, the large deformations have to be accounted for. For granular media, the contact elements changes during the incremental or dynamical process in such a way that it is difficult to appreciate the efficiency of the different parallel strategies. But from a mechanical point of view, the use of parallel computations is necessary to handle a realistic sample.

## References

[1] K. Ach, Etude du comportement de structures polyarticulées avec ou sans jeu: techniques d'homogénéisation et modélisation numérique, Ph.D. thesis, Université Montpellier II, 2000.

[2] K. Ach, P. Alart, M. Barboteu, F. Lebon, B.M. Bodji, Parallel frictional contact algorithms and industrial application, Comput. Meth. Appl. Mech. Eng. 177 (1999) 169 181.

[3] Y. Achdou, P. Le Tallec, F. Nataf, M. Vidrascu, A domain decomposition preconditioner for an advection diffusion problem, Comput. Meth. Appl. Mech. Eng. 184 (2000) 145 170.

[4] P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Comput. Meth. Appl. Mech. Eng. 92 (1991) 353 375.

[5] P. Alart, M. Barboteu, F. Lebon, Solution of frictional contact problems by an EBE preconditioner, Comput. Mech. 20 (1997) 370 378.

[6] P. Alart, Méthode de Newton généralisée en mécanique du contact, J. Math. Pure Appl. 76 (1997) 83 108.

[7] M. Barboteu, Contact, frottement et techniques de calcul parallele, Ph.D. thesis, Université Montpellier II, 1999.

[8] A. Curnier, Q.C. He, A. Klarbring, Continuum mechanics modelling of large deformation contact with friction, in: Contact Mechanics, Plenum Press, New York, 1995.

<footnote>* Not enough memory for one node.</footnote>
[9] P. Le Tallec, Domain decomposition methods in computational mechanics, Comput. Mech. Adv. 1 (1994) 121-220.

[10] P. Le Tallec, M. Vidrascu, Méthodes de décomposition de domaines en calcul des structures, Premier Colloque National en Calcul des Structures, Giens, 1993.

[11] P. Le Tallec, M. Vidrascu, Generalized Neumann-Neumann preconditioners for iterative substructuring, in: Proceedings of the Ninth International Symposium on Domain Decomposition Methods for Partial Differential Equations, Bergen, June 1996, 1998, pp. 413-425.

[12] J. Mandel, Balancing domain decomposition, Commun. Appl. Num. Meth. 9 (1993) 233-241.

[13] J.J. Moreau, On unilateral constraints, friction and plasticity, in: G. Capriz, G. Stampacchia (Eds.), CIME II, New Variational Technol. Math. Phys. (1973) (175-322).

[14] P. Le Tallec, J. Mandel, M. Vidrascu, A Neumann-Neumann domain decomposition algorithm for solving plate and shell problems, SIAM J. Numer. Math. 35 (1998) 836-867.

[15] G. Pietrzak, A. Curnier, Continuum mechanics modelling and augmented Lagrangian formulation of multibody, large deformation frictional contact problems. in: CIMNE, Comput. Plast. (1997) 878-973.