Computational Mechanics 20 (1997) 370±378 Ó Springer-Verlag 1997

Solution of frictional contact problems by an EBE preconditioner P. Alart, M. Barboteu, F. Lebon

370

Abstract We present an Element by Element (EBE) procedure to solve non symmetric linear systems arising from the solution of contact with friction problems. Hybrid formulation is introduced and different types of contact elements are reviewed. To deal with large scale problems the EBE method has proved to be a strongly parallel algorithm. Numerical experiments described in this paper con®rmed the ef®ciency of this speci®c solver.

1 Introduction Frictional contact phenomena are frequent in structural analysis problems. These problems are dif®cult to formulate and even more to solve because they are governed by multivalued tribological laws and some numerical resolutions can lead to unsymmetric operators. The non symmetry becomes crucial for very large problems involving three dimensional discretization and time evolution. This paper shows how to use a simple hybrid formulation together with an ef®cient Element-By-Element (EBE) preconditioned generalized conjugate gradient algorithm when dealing with frictional contact problems. The augmented Lagrangian approach given by Alart and Curnier (1991) generates non linear and non-differentiable systems whose unknowns are node displacements and Lagrangian multipliers identi®ed with contact forces. Because the tangent matrix of the system is non symmetric, non-positive de®nite, ill-conditioned and with zeros on the diagonal, we used generalized conjugate gradient methods introduced by Sonnenveld, Wesseling and de Zeeuw (1985). Appropriate preconditioners are necessary. In this paper, a special EBE procedure for elasticity and another one for frictional contact problems are presented. In Sect. 2, we give a contact ®nite element approach: we review the necessary mathematical background and we introduce the mixed formulation for the frictional contact problem. The main contact elements are described in detail.

Communicated by S. N. Atluri, 15 February 1997

P. Alart, M. Barboteu, F. Lebon Laboratoire de MeÂcanique et GeÂnie Civil Universite Montpellier 2 Pl. E. Bataillon, F-34095 Montpellier Cedex 5 France E-mail:[email protected] Correspondence to: F. Lebon

Section 3 is devoted to conjugate gradient squared methods and to a special EBE preconditioner for only elasticity. In Sect. 4, we extend this last EBE method for frictional contact problems, characterized by a simple ``node on rigid plan'' element. Finally, in Sect. 5, classical test problems are presented such as the square test and the dovetail assembly; we describe the numerical features of the EBE preconditioner; its ef®ciency is discussed and a comparison between standard and EBE preconditioners is made.

2 Contact finite element approach 2.1 Contact and friction by hybrid formulation In the following, only a 2D discretized body is considered and the mechanical laws are written for one contact node. Unilateral contact and friction laws are given by multivalued relationships between static and kinematic variables expressed with respect to a local frame noted …t; n†. Their graphs are drawn in the Fig. 1. The normal contact distance of the node from the obstacle is denoted dn ; dt is the tangential slip increment, kn and kt are the normal and tangential components of the contact nodal force. The adopted contact and friction laws derive from pseudo- (i.e. non-differentiable) potentials, i.e. the two following similar inclusions, dn 2 oWRÿ …kn † dt 2 oWC…kn † …kt † :

…1†

WD denotes the indicator function of a convex set D and oWD the sub-differential of this function. The set C…kn † (the section of the Coulomb's cone) is de®ned as follows with friction coef®cient l,

C…kn † ˆ fkt : kkt k  ÿlkn g

…kn  0† :

…2†

As in a previous paper, based on an augmented Lagrangian approach (Alart and Curnier 1991), the equilibrium of a discretized elastic body in frictional contact with an obstacle or of elastic bodies in frictional contact is governed by the following equation:

G…u† ‡ F…u; k† ˆ 0

…3†

where G…u† is de®ned by:



F int …u† ÿ F ext G…u† ˆ 0



;

…4†

F int …u† and F ext denote the internal and the external forces. F…u; k† de®nes the continuous frictional contact operator which may be written for several contact elements as,

Fig. 1. Unilateral contact and friction laws

F…u; k† ˆ 2 3 ‰ru dn …u†Št projRÿ …rn † ‡ ‰ru dt …u†Št projC…projR rn † …rt † ÿ 4 5;   ÿ 1r k ÿ projRÿ …rn † ‡ projC…projR rn † …rt †

 Firstly, if the target body is rigid a simple ``node on rigid'' treatment is obvious. In this contact case the target body may have a rigid ``curve'' or ``¯at'' surface and it does not need to be discretized.  If a one-to-one correspondence between the boundary nodes of the bodies can be established and maintained through out the contact duration (stick or small slip contact), then a ``node on node'' contact geometry is quite adequate.  Otherwise a ``node on facet'' contact must be used to account for initial mismatching as well as subsequent sliding (moderate slip contact).

For each case, we have a symbolic manner to represent these elements. In general, a frictional contact element is …5† composed of elastic nodes (contactor body node and target with rn and rt called `augmented' multipliers, de®ned as body node(s)) and a last node with an arbitrary position follows, whose degrees of freedom are the contact force comporn ˆ kn ‡ rdn …u†; rt ˆ kt ‡ rdt …u† ; …6† nents. In the following, we give respectively the geometry, the frictional contact operator and the tangent matrix of where r is a positive `penalty' factor. the main contact ®nite elements. In order to handle the whole contact area of the disThe ``node on rigid'' contact element contains 4 degrees of cretized body, Rÿ (respectively C…kn †† must be replaced by freedom (Fig. 2 and Fig. 3), the frictional contact operator the cartesian product Rpÿ (resp. Cp …kn ††, where p is the (5) and the tangent matrix are the following:     number of eventual contact nodes. F…u; k† rM M ; A : F…u; k† ˆ ˆ ce ÿ 1r …k ÿ F…u; k†† M 1r …M ÿ I† 2.2 ÿ

Generalized Newton Method (GNM) The algorithms usually applied to augmented Lagrangian problems in optimization (e.g. Uzawa's scheme) are stable but very slow because they are based on an alternate treatment of the primal and dual variables. We prefer to treat both variables simultaneously through Newton's method. In the Eq. (3) we distinguish two parts involving the pair x ˆ …u; k†, say a differentiable part G and a non differentiable one F. Newton's method may be extended to continuous non-differentiable equations such as (3) in the form ÿ
xi‡1 ˆ xi ÿ …Ki ‡ Ji †ÿ1 G…xi † ‡ F…xi † ; …7† Ki ˆ oG…xi †; Ji 2 oF…xi † ;

…8†

The ``node on node'' contact element needs 6 degrees of freedom (Fig. 4) whose frictional contact operator and tangent matrix are:

2

3 F…u; k† 5; ÿF…u; k† F…u; k† ˆ 4 1 ÿ r …k ÿ F…u; k†† 2 3 rM ÿrM M ÿM 5 : Ace ˆ 4 ÿrM rM 1 M ÿM r …M ÿ I†

where oF…xi † is the generalized Jacobian of F at xi (Clarke 1982). Each region of differentiability (linearity in 2D discretisation) of the operator corresponds to a de®nite contact status for each contact node (gap, stick, forward slip or backward slip). For a given local status, the generalized Jacobian is reduced to the single classical jacobian matrix. These tangent matrices Ji , noted in the following Ace (contact element matrix), are described in the next section. The GNM leads to solve successively linear sys- Fig. 2. Node on rigid contact element tems. The matrix of each system is determined by the global contact status which is a combination of p local contact status. The convergence of the GNM is discussed by Alart (1997).

2.3 Frictional contact elements From a programming standpoint, a contact ``®nite'' element is associated to each node of the elastic contactor body. In this section, we describe in detail the main conFig. 3. Node on curve contact element tact elements by analogy with ®nite elements for solids.

…9†

371

The most classical GCG method is the bi conjugate gradient method (Joly 1990). This method is ef®cient but needs two directions of descent p, q: Initialization

x0 ; y0 ; c 2 Rn r0 ˆ b ÿ Ax0 and s0 ˆ c ÿ At y0

Fig. 4. Node on node contact element

372

p0 ˆ r0

The ``node on facet'' contact element needs 8 degrees of freedom (Fig. 5) characterized by the frictional contact operator and tangent matrix as follows:

2

3 F…u; k† 6 ÿ…1 ÿ a†F…u; k† 7 7 F…u; k† ˆ 6 4 ÿaF…u; k† 5 ; ÿ 1r …k ÿ F…u; k†† 2

rM …a ÿ 1†rM 6 …a ÿ 1†rM …a ÿ 1†2 rM Ace ˆ 6 4 ÿarM ÿa…a ÿ 1†rM M …a ÿ 1†M

ÿarM ÿa…a ÿ 1†rM a2 rM ÿaM

…10† 3

M …a ÿ 1†M 7 7: ÿaM 5 1 r …M ÿ I†

…11† In the previous expressions a is the barycentric coef®cient of the impact (contactor node) over the deformable facet. The submatrix M is a rank-one matrix which depends on the local contact element. Moreover, for each frictional contact status the operator F…u; k† is de®ned as follows:

8 < 0 gap status F…u; k† ˆ rn n ‡ rt t ˆ r stick status : : r …n ÿ elt† slip status n 3 Solving non symmetric systems

q0 ˆ s0 Iterations

For k ˆ 0; 1;


ÿ
ÿ
ak ˆ rk ; sk = Apk ; qk xk‡1 ˆ xk ‡ ak pk and rk‡1 ˆ rk ÿ ak Apk yk‡1 ˆ yk ‡ ak qk and sk‡1 ˆ sk ÿ ak At qk ÿ
bk‡1 ˆ rk‡1 ; sk‡1 =…rk ; sk † pk‡1 ˆ rk‡1 ‡ bk‡1 pk and qk‡1 ˆ sk‡1 ‡ bk‡1 qk Remark 1 It is possible to show that if the coef®cient (Apk ; qk † is always different to zero then the method converges. Then the number of iterations is smaller than the order of the matrix. For our problem, we have only shown that if the friction coef®cient is suf®ciently small, the two ®rst coef®cients are different to zero. The idea of the conjugate gradient squared method (CGS), de®ned by Sonnenveld et al. Zeeuw (1985) to improve the previous algorithm, is to construct a polynomial Sn minimizing a special norm ……Sn …A†r 0 †T Aÿ1 …Sn …A†r 0 ††1=2 under the condition Sn …0† ˆ 1. The resulting algorithm veri®es the condition

rk‡1 ˆ …Sk‡1 †2 …A†r 0 :

…13†

Preconditioning techniques are necessary to accelerate these kind of method. Different preconditioners have been proposed (see for example Dupont et al. 1968; Concus et al. 1982 or Wallis 1983; Hughes et al. 1987; Joly 1990; Alart 3.1.1 and Lebon 1995). Traditionally, these methods are based Algorithm with preconditioner on Jacobi, block Jacobi, SOR iterations, incomplete facFor large and sparse matrices A, it is advisable to use generalized conjugate gradient method (GCG) to solve the torizations or coarse/®ne decomposition. If C denotes the preconditioner matrix, the precondiEq. (7). A general form of this kind of method is given by: tioned CGS method is written: k‡1 k k k k 0 0 0

3.1 Conjugate gradient squared method (CGS)

x

ˆ x ‡ a p and p ˆ Tk …A†r with r ˆ b ÿ Ax : …12† Initialization Tn is a polynomial of degree n in A. We have x0 2 Rn rk‡1 ˆ b ÿ Axk‡1 ˆ Sk‡1 …A†r0 where Sn is a polynomial of 0 0 degree n. Note that the polynomial Sn veri®es Sn …0† ˆ 1. r ˆ b ÿ Ax z 0 ˆ Cÿ1 r 0 p0 ˆ r 0 q0 ˆ r0 Iterations Fig. 5. Node on facet contact element

For k ˆ 0; 1; . . .

yk ˆ Cÿ1 Aqk ÿ
ak ˆ z0 ; zk †=…z0 ; yk ÿ
xk‡1 ˆ xk ‡ ak pk ‡ pk ÿ ak yk and ÿ
rk‡1 ˆ rk ‡ ak A pk ‡ pk ÿ ak yk

tracting the gathering operation out of the round brackets of the Eq. (15). For this purpose, let C1 be the preconditioner matrix obtained after this ®rst approximation:

zk‡1 ˆ Cÿ1 rk‡1 ÿ
ÿ
bk‡1 ˆ z0 ; zk‡1 = z0 ; zk ÿ
pk‡1 ˆ zk‡1 ‡ bk‡1 pk ÿ ak yk and ÿ
qk‡1 ˆ pk ‡ bk‡1 pk ÿ ak yk ‡ bk‡1 qk

with

C1 ˆ W1=2

Nel Y eˆ1

!

 e Š W1=2 ; ‰A

 e ˆ Ie ‡ Wÿ1=2 …Ae ÿ We †Wÿ1=2 : A ge ge

…16†

…17†

Wge is the restriction of the diagonal global matrix W to the degrees of freedom of the eth element and it must be  e is the distinguished from the elementary diagonal We
A Remark 2 If the bi conjugate gradient method converges Winget regularized element matrix. This regularization  e is positive-de®nite, hence a Cholesky factorthen the CGS method converges twice as fast. As for the ensures A previous method the coef®cient …z0 ; yk † is different to zero ization LU is well de®ned. So the preconditioner may be for the two ®rst iterations. written C1 as follows:

3.2 EBE preconditioner for elasticity Hughes et al. (1983) have introduced an element-by-element preconditioner to overcome the excessive cost of the preconditioning matrix C, especially when the scale of the system is too large. In this paper, we generalize the algorithm for our nonsymmetric systems. But, it is interesting to present in detail a EBE preconditioner for a symmetric system involving a LU factorization instead of the LDLt factorization in order to generalize it for the non symmetric case. Moreover some modi®cations and detailed information are introduced in comparison to the original method of Hughes and Winget (1985). Element-by-element preconditioners are motivated by a desire to maintain the element based data structure of current ®nite element codes. In order to arrive at the Cholesky EBE preconditioner, we start from the standard ®nite element method assembly operation. If we denote the eth element contribution to the global matrix A by Ae , we have: Aˆ

Nel Y eˆ1

Ae ;

…14†

C1 ˆ W1=2

Nel Y eˆ1

!

‰Le :Ue Š W1=2 ;

with

 e Pe Š†Pe and Ue ˆ Pe …U  e Pe Š†Pe ; e ‰Pe A  e ‰Pe A Le ˆ Pe …L …19†  e Pe Š and U  e Pe Š are given by the Cholesky e ‰Pe A  e ‰Pe A where L  e Pe
Pe is a perfactorization of the product matrix Pe A mutation matrix whose role is to interchange rows and  e such that they are consistent with the local columns of A nodal ordering. It is convenient to express Le and Ue by using a decomposition into upper-triangular (F†, lower-triangular (E†, and diagonal (D† parts, with respect to the global ordering: Le ˆ Ee ‡ Ie ;

Ue ˆ De ‡ Fe :

eˆ1

ˆW

ÿ1=2

I‡W

Nel Y eˆ1

…21†

ˆ W1=2 …E ‡ D ‡ F†W1=2 :

eˆ1

1=2

…20†

Consequently, the matrices Ee ; Fe are not necessarily triangular according the local element nodes ordering. From a computational point of view, the matrices Le and Ue are stored in a single array as a matrix; the identity matrix is implicit. The gathering operation deals with Ee ; Fe ; De as a single standard matrix, hence the global splitting E; F; D is relevant:

! and Nel is the number of elements. This global matrix can Nel Y 1=2 be decomposed as follows: C2 ˆ W ‰Ee ‡ De ‡ Fe Š W1=2 Nel Y A ˆ W ‡ ‰Ae ÿ We Š

…18†

! ÿ1=2

‰Ae ÿ We ŠW

1=2

W

; …15†

where W (resp. We ) is the diagonal part of A (resp. Ae ). The last factorization with the diagonal matrix W has a regularizing effect. Indeed, for diagonal scaling preconditioner, this regularization nondimentionalizes the equation system, ensuring that the residual norm used to monitor convergence is well-de®ned. The ®rst approximation for reaching this preconditioner consists in ex-

Since the sum decomposition E ‡ D ‡ F corresponds to the computational storage of the LU decomposition, we can approximate the preconditioner C2 by the new one noted C3 :

    C3 ˆ W1=2 E ‡ …D ‡ F†W1=2

…22†

The last sum decomposition (22) may be identi®ed with an approximation based upon the global Cholesky factorization. This allows us to de®ne a EBE preconditioner C4 :

   C4 ˆ W1=2 L UW1=2 ˆ L
U :

…23†

373

!

Nel Some experiments show that a ®nal modi®cation improves Y ~ 1=2 J ‡ W ~ 1=2 ; ~ ÿ1=2 ‰Ae ÿ W  e ŠW ~ ÿ1=2 W the behaviour of the preconditioner: the diagonal part of ˆ W eˆ1 U, noted D, is replaced by the initial diagonal part of A (noted W). This last approximation is a main advance added to the standard EBE preconditioner. Then we pos- with tulate the ®nal EBE preconditioner C5 as: 

1 if i 2 ce ~ ii ˆ ; W = ce Wii if i 2   It seems that this modi®ed preconditioner takes advantage 0 if i 2 ce 0 if e ˆ ce  of the classical EBE preconditioner and the standard diand We ˆ Jii ˆ ; 1 if i 2 = ce We if e ˆ ee agonal one (Barboteu et al. 1997 in preparation). Indeed the convergence of the residual norm in the conjugate gradient squared algorithm is more stable, i.e. this norm where i represents a degree of freedom from a node and ee (resp. ce) represents the ®nite elastic (resp. contact) eledecrease is monotonous. ment. Remark 3 In elasticity, the stiffness matrix A is symmetFor the contact and friction problems, we can take the ric. It is then appropriate to use a Crout LDLt decompo- same structure as the elastic preconditioner algorithm sition. The Crout and Cholesky EBE preconditioner construction. So the ®rst approximation in (16) is modialgorithms are obviously quite similar in structure, except ®ed as follows: the computational storage. The matrix C3 (resp. C4 ) de®ned in (15) (resp. in (18)) must be replaced by the folNel h i Y ~ 1=2 ; ~ 1=2 ‡ ~ ÿ1=2  ~ ÿ1=2 W lowing expression: C1 ˆ W Je ‡ W ge …Ae ÿ We †Wge     eˆ1 C3 ˆ W1=2 E ‡D‡ Et W1=2 …resp. C4 ˆ L
D
Lt † : …28† C5 ˆ L…U ÿ D ‡ W† :

374

…27†

…24†

…25† Cholesky EBE requires more operations per iterations than Crout EBE because the elements diagonals Dii cannot be multiplied into a global diagonal scaling matrix W1=2 (see Hughes et al. 1987). The substitution in (24) amounts to replace D (the diagonal part of L
D
Lt decomposition) by W, then we end up with the following preconditioner C5 :

C5 ˆ L
W
Lt :

…26†

with

Je ˆ

n

0 I

if e ˆ ce if e ˆ ee

and

~ ge ˆ W



I Wge

if e ˆ ce : if e ˆ ee

According to the section 3.2 we know that the Winget  ee , as mentioned in the relation (30) leads regularization A to a LU decomposition of Aee . Moreover we can see in the next section that the factorization is made directly on the elementary contact matrices Ace . Then we get the next formula:

! ! 4 Nel Nel Y Y e W ~ 1=2 ˆ W ~ 1=2 ~ 1=2 ~ 1=2 ; Extension to hybrid contact elements (node on rigid) A C1 ˆ W ‰Le
Ue Š W To deal with non symmetric systems issued from our eˆ1 eˆ1 formulation of frictional contact problems, it is necessary …29† to modify the previous global scheme. Moreover the LU where factorization of the elementary contact matrices is not obvious and it must be discussed for each contact status in  e ˆ Le
Ue A a second section. (  ee ˆ Ie ‡ W ~ ÿ1=2 ~ ÿ1=2 ˆ Lee :Uee A ge …Ae ÿ We †W ge ; ˆ 4.1 Ace ˆ Lce :Uce EBE preconditioner for hybrid element …30† The ®rst dif®culty comes from the eventual zeros on the diagonal part of the matrix A due to the hybrid elements. Consequently the matrix W cannot be used as previously. where Le and Ue are the Cholesky decomposition terms  e. Moreover the Winget regularization is not appropriate for of A We use a similar method as the one described for the elementary contact matrices; indeed, some of hybrid maEBE preconditioner computation in the previous section. trices may be singular, according to the contact status (they do not have the appropriate properties as the ®nite The last modi®cation concerns the substitution of some diagonal terms in order to de®ne the appropriate candielement matrix). The decomposition of the matrix A is date C5 . Note that we only replace the diagonal contribu~  performed by introducing the matrices W; J and We tion which is not related to the contact and friction: instead of W; I and We : ÿ
Nel ~ ‡W  ; Y C5 ˆ L
U ÿ D …31† ~ ‡ ‰Ae ÿ W  eŠ AW eˆ1 with

Table 1. Correspondence between the standard and frictional contact procedure Global matrices Standard procedure Frictional contact procedure

~ ii ˆ D



0 Uii

if i 2 ce ; if i 2 = ce

I J

 ii ˆ W

Local matrices

W ~ W  W;



0 Aii

L L

if i 2 ce : if i 2 = ce

We can summarize in Table 1 the main modi®cations in the procedure in order to obtain a speci®c frictional contact preconditioner in comparison with the standard elastic one.

D ~ D

M11 6ˆ 0

U U

We e W

Ie Je

Wge ~ ge W

8 d1 ˆ M11 > > < d ˆ det M ˆ 0 2 M11 : 21 > a1 ˆ M > M11 : a2 ˆ M12

…37† 375

 If M11 is equal to zero and since M is a rank-one matrix, we can distinguish two cases.  For M21 is zero, M is an upper triangular matrix and an obvious Cholesky decomposition may be written with 4.2 L ˆ I; U ˆ M. LU factorization of elementary contact matrices  If M12 vanishes, M is lower triangular, then we postulate The LU factorization of the elementary contact matrices is for an approximated L U decomposition as follows: not standard and is presented for the ``node on rigid'' L ˆ M; U ˆ I. Indeed this factorization is not standard element (see Sect. 3.2). The elementary matrix of such because the diagonal part of L is not the identity. element has the form below: Consequently, we perform so a new Cholesky approxi  rM M mation of the matrix Ace without modifying imple; …32† Ace ˆ mentation, because the diagonal part of L is implicit as M 1r …M ÿ I† mentioned above. We can summarize these two special where the rank-one matrix M depends on the frictional cases: contact status. Although two of these matrices Ace are    I 0 rM M singular, it is possible to ®nd a Cholesky LU decomposi; M11 ˆ 0; and M21 ˆ 0 Aec ˆ I I 0 ÿ rI tion. For the gap and stick status, a factorization of the r elementary matrix Ace is easily obtained: …38†

gap status …i.e., M ˆ 0†; " #  # " 0 0 I 0 0 0 ˆ ˆ Lce
Uce ; Ace ˆ 0 ÿ rI 0 I 0 ÿ rI

…33† stick status …i.e., M ˆ I†; #" #   " I 0 rI I rI I ˆ Lce
Uce : ˆ I Ace ˆ I 0 ÿ rI I 0 r …34† t For the slip status (i.e., M ˆ …n ÿ elt†n †, the situation is less simple. If we ®nd a Cholesky decomposition of the submatrix M (noted LU), it is possible to determine a factorization of the frictional contact matrix Ace as follows:      L 0 rU U rM M ˆ : Ace ˆ L 0 ÿ rI M 1r …M ÿ I† I r …35† Then the problem amounts to study the existence of the LU factorization of M, as it is presented below. The submatrices L and U are determined by four coef®cients a1 ; a2 ; d1 ; d2 :  L ˆ I ‡ a1 e2 et1 : …36† U ˆ d1 e1 et1 ‡ d2 e2 et2 ‡ a2 e1 et2  If the ®rst component of M, noted M11 , does not vanish, the coef®cients of L and U are clearly identi®ed:



M11 ˆ 0; and M12 ˆ 0

M 0 Aec ˆ M I r



rI I 0 ÿ rI



:

…39† 4.3 Some extensions and topics A straightforward generalization concerns the elementary matrices issued from the ``node of facet'' contact element (see Sect. 2.3), which have the following form: Aec ˆ 2 rM 6 …a ÿ 1†rM 6 6 4 ÿarM M

…a ÿ 1†rM …a ÿ 1†2 rM ÿa…a ÿ 1†rM …a ÿ 1†M

ÿarM

M

3

ÿa…a ÿ 1†rM …a ÿ 1†M 7 7 7 ÿaM 5 a2 rM 1 ÿaM r …M ÿ I† …40†

It is easy to verify that the following factorization exists by assuming the decomposition of M,

Aec ˆ 2 32 rU …a ÿ 1†rU ÿarU L 6 …a ÿ 1†L I 76 0 0 6 76 6 76 4 ÿaL 5 4 0 0 I L 0 0 I r

3 U 0 7 7 7 0 5 ÿ rI …41†

5 Numerical results

376

5.1 Parametric study on a classical benchmark We have tested the previous frictional contact EBE algorithm on the elastostatic square test. This example is very elementary but with complicated contact status on the contact area according to the loadings and the frictional coef®cient (see Raous et al. 1988; Alart and Lebon 1995 and Curnier and Alart 1988). This application involves a long bar pressed on a plane. The coarse mesh contains 16 contact nodes. The ®ner meshes are constructed by a global re®nement, i.e. each ®nite element is divided into four new elements. So we obtained meshes with 16, 32, 64, 128, 256 contact nodes and respectively 220, 780, 2620, 7980, 23 200 degrees of freedom (see Fig. 6). Consequently, the proportion of the degrees of freedom associated to a contact node or a multiplier node (therefore concerned by the non symmetry) decreases when the mesh in ®ner: 30%, 16%, 10%, 6%, 4%. To appreciate the ef®ciency of the EBE preconditioner we compare the EBE preconditioner with the Diagonal and ILU method in terms of CGS iterations number and CPU time. Usually, the EBE preconditioner requires the same number of iterations as the ILU method (Figs. 7, 8 and 9). The gap is moderate and decreases when the friction coef®cient increases (see Fig. 7). It means that the EBE method is less sensitive to the non symmetry emphasized by the frictional coef®cient. Figure 10 con®rms this interpretation because, for a same number of degrees of freedom, the EBE preconditioner is more ef®cient if contact occurs. On the other hand, for a suf®ciently large number of degrees of freedom (>104 ), the ILU preconditioner requires more number iterations when there are contact elements (see Fig. 10). The EBE method is more sensitive to the penalty factor (Fig. 8). It seems that an optimal value of this factor exists, around 0.1 multiplied by the Young modulus, for our example.

Fig. 6. Square test (®ne mesh)

Fig. 7. In¯uence of the friction coef®cient on preconditioner/CGS convergence behaviour (square test with 64 contact nodes and r ˆ 0:1)

Fig. 8. In¯uence of the penality factor on preconditioner/CGS convergence behaviour (square test with 64 contact nodes and l ˆ 0:1)

Fig. 9. Comparison between preconditioner/CGS convergence behaviour according to the degrees of freedom (square test with r ˆ 0:1 and l ˆ 0:1)

377

Fig. 10. In¯uence of the frictional contact from the elasticity case Fig. 11. Dovetail assembly (®ne mesh) on EBE and ILU/CGS convergence behaviour according to the degrees of freedom (r ˆ 0:1 and l ˆ 0:1)

effected with a coarse mesh (2362 nodes with 128 contact nodes) and a ®ne mesh (14 354 nodes with 512 contact These comments are con®rmed by the study of the CPU nodes) and for a frictional coef®cient (resp. the penalty times (see Tables 2 and 3). factor) equal to 0.3 (resp. 0.1). The previous comments are not modi®ed by the oblique character of the contact area for which the elementary 5.2 factorization is different (see Sect. 4.2). This example leads Dovetail assembly This example, which concerns an industrial problem, is a to another interesting conclusion: the EBE preconditioner two body contact and friction problem with two oblique is more ef®cient than ILU preconditioner for a ®ne mesh contact zones (see Fig. 11). The numerical tests have been (Table 5 versus Table 4). Table 2. Average CPU times per newton iterations according to the problem size

EBE 16 contact nodes Average CPU times per newton iterations (sequential)

32 64 128 256

± ± ± ±

Table 4. Numerical results with the dovetail assembly (128 contact nodes)

Table 5. Numerical results with the dovetail assembly (512 contact nodes)

CGS average iterations CPU times (sequential)

CGS average iterations CPU times (sequential)

1.8

0.73

DIAG 1.3

7.6 33.3 164 414

4.4 25.6 123 343.5

8.2 50 425 no conv

EBE

ILU

DIAG

16 contact nodes

6.7

2.3

3.7

32 64 128 256

22 93 455 1113

12 63 340 978

22 142 1068 no conv

Table 3. Total CPU times according to the problem size Total CPU times (sequential)

ILU

± ± ± ±

per newton iterations total

per newton iterations total

EBE

ILU

DIAG

287 49 235.2

268 46.8 223.6

no conv ± ±

EBE

ILU

DIAG

740.4 763.6 4313

884.2 1045.6 5854

no conv ± ±

378

6 Conclusion Our results indicate that the EBE preconditioner can be competitive with the ILU preconditioner in terms of CGS iteration number and CPU time. It seems to be less sensitive to the non symmetry due to contact elements. Moreover the EBE procedure has been introduced for vectorized computation (see for example Muller 1985; Nour-Omid and Parlett 1985 and Hughes et al. 1987) and, presently, may be implemented on parallel computers which use the aggregate power of many processors. This algorithm is particularly well adapted to computers with a large number of processors, at best one processor per element. For a few processors, the mesh may be divided into sub-structures; but in this case, the EBE procedure will have to be compared with domain decomposition methods which are known as very ef®cient algorithms. References

Alart, P. (1997): Methode de Newton geÂneÂraliseÂe en meÂcanique du contact. Journal de MatheÂmatiques Pures et AppliqueÂes (to appear) Alart, P.; Curnier, A. (1991): A mixed formulation for frictional contact problems prone to Newton like solution methods. Comp. Meth. Appl. Mech. Eng. 92, 353±375 Alart, P.; Lebon, F. (1995): Solutions of frictional contact problems using ILU and coarse/®ne preconditioners. Comp. Mech. 16, 98±105 Barboteu, M.; Alart, P.; Lebon, F. (1997): An improved EBE preconditioner for elastostatics, in preparation Behie, A.; Forsyth, P. A. (1983): Comparison of fast iterative methods for symmetric systems. Siam. J. Num. Anal. 3, 41±63 Clarke, F. H. (1982): Optimization and nonsmooth analysis. Wiley Ed., New-York Concus, P.; Golub, G.; Meurant, G. (1982): Block preconditioning for the conjugate gradient method. Report LBL 14856, University of California

Curnier, A.; Alart, P. (1988): A generalized Newton method for contact problems with friction. J. Theo. Appl. Mech., Special edition: Number. Methods Mech. Contact Involving Friction, 7, Suppl. No. 1 Dupont, T.; Kendall, R. P.; Rachford, H. H. (1968): An approximate factorization procedure for solving self-adjoint elliptic difference equations. Siam, J. Num. Anal. 9, 559±573 Hughes, T. J. R.; Ferencz, R. M.; Hallquist, J.O. (1987): Largescale vectorized implicit calculations in solid mechanics on a Cray X-MP/48 utilizing EBE preconditioned conjugate gradients. Comp. Meth. Appl. Mech. Eng. 61, 215±248 Hughes,T. J. R.; Levit, I.; Winget, J. M. (1983): An element by element solution algorithm for problems of structural and solid mechanics. Comp. Meth. Appl. Mech. Eng. 36, 241±254 Hughes, T. J. R.; Winget, J. M. (1985): Solution algorithms for non linear transient heat conduction analysis employing element by element iteratives strategies. Comp. Meth. Appl. Mech. Eng. 52, 711±815 Joly, P. (1990): Mise en oeuvre de la meÂthode des eÂleÂments ®nis. Ellipse, Paris (in French) Muller, A. M. (1985): Element by element iterative procedures in structural ®nite analysis. Thesis, Standford University, Standford, California Nour-Omid, B.; Parlett, B. N. (1985): Element preconditioning using splitting techniques. SIAM J. Scien. Stat. Comp. 6, 761± 770 Raous, M.; Chabrand, P.; Lebon, F. (1988): Numerical methods for frictional contact problems and applications. J. Theo. Appl. Mech., Special issue, 7, Suppl. No. 1 Sonnenveld, P.; Wesseling, P.; de Zeeuw, P. M. (1985): Multigrid and conjugate gradient methods as convergence acceleration techniques. Multigrid Meth. Integr. different. Clarendon Press, Oxford, 117±167 Wallis, J. R. (1983): Incomplete Gaussian elimination as a preconditioning for generalized conjugate gradient acceleration. In: Proc. 7th Symposium on numerical simulation of reservoir performance of the society of petroleum engineers of AIME, San Francisco