A domain decomposition strategy for nonclassical frictional multi-contact problems

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, nonsymmetric 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.

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 diuse (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 ecient 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 nondierentiable, 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 diculties, we develop, in Section 3 a domain decomposition strategy 1 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 ®nite 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 de®ne speci®c 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 de¯exion d AEg. So these two unilateral condition relate the so-called gap distance d to the contact shear stress s 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): where oW G denotes the subdierential of the indicator function W G of the interval G Àg; g. W Ã G is the Legendre Fenchel form of W 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 coef®cient. This friction law is de®ned for a contact   shear stress s 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 b t (according to the tangent to the hinge) and the friction shear s (decomposed in a contact shear stress s e and a frictional torque s t ). The graph of this multi-valued tribological law, represented in Fig. 3, summarizes the relationship between s t and b t . According to convex analysis [13], the friction law can be written as follows: where W Ã C s e b t is the Fenchel conjugate form of the indicator function W Cse . Cs e denotes the interval Àljs e j; ljs e j.
Frictional contact law. A law of frictional bi-unilateral contact can be constructed by combining the biunilateral contact law with the frictional rotative law. This (nonassociated) frictional contact law can be written in the form (1) plus (2), provided the dependance of the friction disc Cs e upon the shear contact stress s e is taken into account, So this law is nonsymmetric: the contact stress is independent of the friction law but the inverse is false. We can then ®nd the associated frictional contact status of a particle in the hinge: À g < d < g; s e 0 gap; d g; s e > 0; j b t j 0; js t j 6 ljs e j contact; stick; d g; s e > 0; j b t j > 0; js t j ljs e j contact; slip; d Àg; s e < 0; j b t j 0; js t j 6 ljs e j contactÀ; stick; d Àg; s e < 0; j b t j > 0; js t j ljs e j contactÀ; slip: This tribological interface law can be introduced to complete the formulation of the mechanical behavior of the plates system X 1 X 2 composed by the contactor and target plates (Fig. 1). A Kirchho Love model is selected for the behavior of the two plates.

Hybrid formulation
We focus our attention on the incremental problem which may be formulated as a minimization problem: I sys x; s e : inf where I sys is the total potential energy of this two jointed plates system. V 1 and V 2 denote the virtual¯exion displacement spaces of contactor and target plates.x (resp.x Ã ) is the real (resp. virtual) de¯exion of the global system X 1 X 2 and is de®ned byx x 1 ; x 2 (resp.x Ã x Ã 1 ; x Ã 2 . I e1 and I e2 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 eects in the hinge resulting from the jointed plates reciprocal interactions. I e and I c are de®ned by where D is the stiffness tensor and j is the tensor of changes of curvature of plates. In the expression of I c , the ®rst 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 s e which depends on the solutionx. For this reason, the minimization problem (4) is considered as a quasi-variational problem. To overcome these dif®culties, 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 quasiaugmented Lagrangian is postulated [1,4,7] for the frictional contact problem; it gives an unconstrained formulation: In the context of contact mechanics, the Lagrange multipliers m, c (m Ã ; c Ã for the virtual variables) can be interpreted as the nominal contact shear stress and the nominal frictional torque, respectively. The terms l q e and l r t represent the regularization of the functions W G d and W Ã Csex c, where q 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Ĉ for the original one Cs e : Cm; dx : Cprox W Ã qG m q dx: This approach of``quasi''-augmented Lagrangian permits to satisfy exactly the contact constraints and friction criteria contrary to penalty techniques. The dierentiability regions on the augmented Lagrangien are presented in Fig. 4.
A solution x; k, k m; c of the quasi-minimax problem is characterized by vanishing of the ®rst variations of L r;q sys DI e x; dx DL q;r c x; k;x; k; dx; dk 0; Vdx; Vdk; 7 where the contact part has the following expression, dC: The calculation of the gradient gives The discretization of frictional contact interface between the two plates is performed by a speci®c hybrid ®nite 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 ®nite 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 ®ctitious node numbered q 1 whose degrees of freedom are the frictional contact force components.

5
We denote by e u (or u if no ambiguity arises) the generalized displacements of the eth element and e k (or k) the generalized contact forces: e u u u 1 ; u 2 ; . . . ; u q ; e k \ e u q1 " k: We have to de®ne a nonclassical contact element to model the discretization of the hinge with clearance and rotative friction. The discretization of the rolling shutters uses elastic ®nite plates (DKTP elements) which have three nodes with three degrees of freedom for each node, the de¯exion x, the two rotations b x and b y with respect to the axis x and y. Similarly to the DKTP elastic ®nite element, this hinge element involves three nodes (the contactor, the target node and the contact stress node). We denote byx the vector concatenating the de¯exion x and the two rotations b x and b y , k the vector concatenating the contact force m and the two torque components c x and c y : where e denotes the normal vector to the plane of the plate and t the tangent vector to the hinge in the plane of the plate (normal to the vector n).
The kinematic contact variables are the gap dx and the relative rotation increment dbx between two loading increments i and i 1, By using the expression of the gradients given by (8), we get the form of the elementary frictional contact operator Fx 1 ; x 2 ; k: where r m q dx e T k qx 1 À x 2 ; s c rdb: r and s are the augmented Lagrange multipliers associated, respectively, to the contact force and the rotative frictional torque (q and r are the penalty parameters). The function f q; r is equal to q for the de¯exion degrees of freedom and to r for the rotations. Moreover with the gradients of the gap and the relative rotation increment (rxg e x e T ; Àe T and rxdb t tt T ; Àtt T ), the elementary frictional contact operator has the following expression: 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 F x 1 ; x 2 ; k corresponding to these contact status, takes the following form: (e equals to 1 or À1 with respect to the slip direction ( or À)). The tangent matrix A ec r x 1 ;x 2 ;k Fx 1 ; x 2 ; k associated to this operator has the following form: with M 0 gap status; I stick status for the play g or À g; e À elte T slip status for the play g or À g: V X 13

Generalized Newton method (GNM)
In the augmented Lagrangian approach, the equilibrium of a discretized contact bodies system is characterized by the ®rst variation of L r sys (7) which takes the form, Gx Fx 0: 14 So we distinguish two parts involving the pair x u; k, a dierentiable part Gx Gx and a nondierentiable one Fx. 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 nondierentiability of Eq. (14), Newton's method may be extended to the following iterative form [4,6]: The matrix K i is the usual elastic stiness matrix and A i c represents a tangent matrix issued from the assembly of the tangent matrix A ec (de®ned 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: The linear systems may be solved by an iterative method as the square preconditioner conjuguate 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 particular speci®c 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¯exion 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 ®t into each other. To facilitate the rolling of the shutters at the opening, the pro®le 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.

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 in¯uence 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 de¯exion passes from 15.24 cm for the classical hinge to 17.02 cm for our speci®c hinge model. For more mechanical considerations (study of the in¯uence of the friction or the play) the reader can refer to [1,7]. To carry out the mechanical study, we have met some diculties 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 speci®c strategies well suited to the solution of large scale problems with a large number of contact zones. Moreover to overcome the diculties 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.

Nonlinear domain decomposition strategies for frictional multi-contact problems
Domain decomposition methods (substructuring techniques) are ecient 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 ecient solvers in a classical monoprocessor 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 X into smaller subdomains X 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,

Ku F 17
with K a square, symmetric positive de®nite matrix arising from ®nite element discretization of a linear, elliptic, self-adjoint boundary value problem on a domain X. The ®rst step of the Schur complement method consists then in splitting the domain into small local nonoverlapping subdomains with interfaces de®ned as follows (see Fig. 10): With this partition we denote by K n and F n the subdomain stiness 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 K n ; u n and F n take this form: 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: where S N n1 R n t S n R n ; S n K n ll À K n li K n ii À1 K n il ; B N n1 R n t F n l À K n li K n ii À1 F n i ; B n F n l À K n li K n ii À1 F n i : 20 R n is the restriction operator which goes from C to C n . The matrix S is the Schur complement matrix; S n are the local Schur complement matrices.

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- 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 (condS O1=H 2 1 H =h, where H is the subdomain size and h is the mesh size) than the initial system (con-dK O1=h 2 . Hereafter, we use the multi-level Neumann Neumann preconditioner. This iterative technique never requires the explicit calculation of the matrix S. We have just to form the matrix vector products Sp k and Mr k1 by solving independent auxiliary Dirichlet and Neumann problems on the local subdomains. We will brie¯y 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 S N n1 R n t S n R n is approximated by a weighted sum of the matrices S n À1 , M N n1 R n t P n S n À1 P n R n : 21 P n is a diagonal weighting matrix which must verify N n1 R n t P n R n Idj C . A generic choice is: P n K n ii = ndi n1 K n ii 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), S n À1 denotes an arbitary regularized inverse. When this preconditioner is used the condition number is (condMS O1=H 2 1 logH =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¯oating 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 M and of the residual gradient r k1 have the following form: where S n À1 S n À1 if X n is not a floating subdomain; the inverse of the projection of the image of S n ; & and Gc 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 C in order to minimize the residual: Its minimum is attained for the function Gc which vanishes its gradient; so we have to solve the following coarse problem, where c is the unknown, R n t P n S n À1 P n R n r k1 : 23 For the balancing method, the condition number is asymptotically optimal (condM S O1 logH =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].

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 dierence compared to a standard Newton procedure consists of speci®c adaptations to take into account new rigid body motions in¯oating subdomains with contact zones and the nonsymmetry of the coarse problem (Section 3.2.2).

Newton Schur algorithm: general purpose
The algorithm. To extend domain decomposition method to frictional contact problem we have to overcome many dif®culties: the nonlinearity, the non differentialibility and the nonsymmetry. The nonlinearity 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 perconditioner (with and without coarse spaces) is carried to overcome the poor conditioning and the treatment of¯oating subdomains; but we will see in the Section 3.2.2 that the balancing method is sensitive to the nonsymmetry and requires speci®c adaptations. The Newton Schur algorithm takes the form described in Table 2.
Subtracturation stratgies. 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 ef®ciency.

Newton Schur algorithm: speci®c adaptations
The balancing method requires two speci®c adaptations due to frictional contact conditions

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 S n À1 is solved by subdomain. As already mentioned, S n À1 denotes an arbitrary regularized inverse. The actual choice, is to consider the exact inverse inverse for non¯oating subdomains. For¯oating subdomains the construction of a regularized inverse involves the use of a set of elements which de®ne the null space. This choice has an in¯uence on the size of the coarse Table 2 Newton Schur algorithm for frictional contact problems 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 eciency 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¯oating 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.

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 de®ned for nonsymmetric problems. The ®rst idea is to replace the matrix S by the symmetrized matrix S s S s S S t . 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 l 0. We consider the interface reduced matrix S Ã with a zero friction coecient (S Ã S l0 ) to evaluate the norm of the dierence between M and S À1 . Then the minimization problem takes the following form: This minimum is reached for the function Gc which vanishes its gradient, The solution of this equation veri®es the following equality: R n P n S n À1 P n R n r k1 ; 25 which de®nes the coarse problem adapted to the nonsymmetric case.

Numerical results
In this section, we study the numerical and parallel behavior of the Newton Schur strategy. In the ®rst section we present a parametric study to analyse the convergence behavior of this algorithm and speci®c 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 eciency of the parallel implementation of our algorithm.

Parametric study
We consider two decompositions for the rolling shutters: a band decomposition without¯oating subdomains and a square decomposition with¯oating subdomains which permits to validate our speci®c adaptations.

Band decomposition
The in¯uence of three parameters is studied: the number of hinges, the friction and the number of subdomains.
In¯uence of the number of hinges. The in¯uence 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 14 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.
In¯uence 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 (l 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.
In¯uence 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.

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 coecient and the number of subdomains.
In¯uence of the friction coef®cient. Table 5 underlines the contribution of the balancing method when oating 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 magni®es the degradation of the algorithm due to¯oating subdomains. In order to understand the in¯uence of the nonsymmetry on the algorithm and the advantages of the speci®c improvements presented in Section 3.2.2, we study friction coef®cient varying from 0 to 2 for a given example with 16 slats and 15 Â 2 subdomains (26¯oating subdomains) (Fig. 14). If the friction coef®cient is zero, the problem is symmetric and the two balancing method are equivalent and ef®cient. As soon as we introduce a small friction coef®cient, the number of GMRes iterations increases quickly because nearly all contact elements have a slip status which introduce the nonsymmetry in the matrices. For l 0:2, we have the higher ratio of slip status. For l > 0:2, the ratio of sticking status increases and the nonsymmetry decreases. For l > 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 coef®cient is weak. Notice that the speci®c balancing method is less sensitive to the nonsymmetry: it requires twice less iterations than the standard one.
In¯uence of the number of subdomains. The good behavior of the speci®c preconditioner is con®rmed when the number of¯oating 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 ®rst conclusion, the algorithm has the desired numerical scalability. Thus the strategy developed above is ef®cient for large scale frictional contact problems.

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- 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 de®ned 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 dierent; 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 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 de®nition 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.

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 diusion 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 speci®c diculties. 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 dicult to appreciate the eciency of the dierent parallel strategies. But from a mechanical point of view, the use of parallel computations is necessary to handle a realistic sample.