Qualitative analysis of solutions to discrete static contact problems with Coulomb friction

Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Qualitative analysis of solutions to discrete static contact problems with Coulomb friction J. Haslinger, V. Janovsky´, T. Ligursky´ ⇑ Charles University in Prague, Faculty of Mathematics and Physics, Department of Numerical Mathematics, Sokolovská 83, 186 75 Prague 8, Czech Republic

a r t i c l e

i n f o

Article history: Received 13 April 2010 Accepted 15 September 2010 Available online 29 September 2010 Keywords: Contact problem Coulomb friction Mixed finite element formulation Local Lipschitz continuous branches of solutions Piecewise smooth continuation methods

a b s t r a c t We analyze properties of solutions to discrete contact problems with Coulomb friction which are parametrized by the coefficient of friction F . Using a generalized variant of the implicit-function theorem we establish conditions under which there exists a local Lipschitz continuous branch of solutions around a reference point. Finally, a piecewise smooth continuation algorithm which allows to follow such branches of solutions is proposed.  2010 Elsevier B.V. All rights reserved.

1. Introduction Contact mechanics is a branch of solid mechanics which studies the behavior of loaded deformable structures in mutual contact. Besides non-penetration conditions one should take into account the influence of friction on contacting surfaces to get a more realistic model. Although Coulomb friction is the classical one, its mathematical treatment remained open for a long time. The first mathematically justified existence proof for the static case was done in [1]. Later, results were extended to quasistatic and dynamic problems in [2]. Typically the existence of at least one solution is shown for sufficiently small coefficients of friction F . On the other hand there is no information on the structure of solutions in a general case. Some partial results are known provided that the solution enjoys a priori given properties ([3,4]). The situation is rather different in the discrete case. In examples with a small number of degrees of freedom one can find solutions ‘‘by hand”. The structure of solutions is relatively complicated even for models with one degree of freedom ([5–7]). For models arising from finite element discretizations the situation is more involved but still we have certain knowledge on the qualitative behavior of their solutions which is not available in the continuous case. Indeed, the discrete model has at least one solution for any friction coefficient F , and this solution is unique if F is small enough. On the other hand, the bounds on F guaranteeing uniqueness of the solution are ⇑ Corresponding author. Tel.: +420 22191 3364; fax: +420 22481 1036. E-mail addresses: [email protected] (J. Haslinger), [email protected]. cuni.cz (V. Janovsky´), [email protected] (T. Ligursky´). 0045-7825/$ - see front matter  2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2010.09.010

mesh-dependent. In other words, the uniqueness property of the discrete problem depends on the size of the problem. Thus it may happen that a small change of F leads to a large change of the solution but this phenomenon appears only when sufficiently fine discretizations are used. This ‘‘catastrophic” behavior can be observed in numerical experiments ([8,9]). The first step towards better understanding of the discrete model is to establish the existence of local Lipschitz continuous branches of solutions parametrized by F . The present paper extends results from [10], where F was assumed to be constant, to the case when F depends on the spatial variable. Such results are important for the development of continuation algorithms, which is another goal of this paper. The paper is organized as follows: in Sections 2 and 3 we formulate a 3D static contact problem with orthotropic Coulomb friction and solution-dependent friction coefficients. We use the fixedpoint formulation of both, the continuous as well as the discrete problem. We show that the discrete problems have a solution for any friction matrix F belonging to an appropriate class and have a unique solution under additional assumptions on F . Then we analyze how these assumptions depend on the size of the discrete problem, i.e. on the norms of finite element partitions. The new results extend and cover the classical ones for isotropic Coulomb friction with a coefficient which does not depend on the solution ([11]). In Section 4 we formulate sufficient conditions for isotropic Coulomb friction in 2D under which there exists a Lipschitz continuous solution branch in a vicinity of a reference value of F . In Section 5 we propose a piecewise smooth variant of the Moore-Penrose continuation algorithm, which enables to follow

150

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

the solution branch. Finally in Section 6 we apply this algorithm for solving very simple model examples with a small number of contact nodes. Throughout the paper we shall use the following notation: Hk(D), k P 0 integer, stands for the standard Sobolev space of functions in the set D (H0(D) = L2(D)) equipped with the norm k kk,D and the scalar product ( , )k,D. For vectors, vector functions and matrices we use bold characters. The scalar product of two vectors x; y 2 Rp is denoted by x  y. The Euclidean norm in Rp , which is frequently used in the text, is simply written as k k. The set of all x 2 Rp such that xi P 0, xi > 0 "i = 1,. . .,p will be denoted by Rpþ and Rpþþ , respectively. Finally, the fact that a constant C depends on parameters t1,. . .,ts will be emphasized by writing C := C(t1,. . .,ts). 2. Fixed-point formulation of contact problems with orthotropic Coulomb friction law Let a deformable body be represented by a bounded domain

problem we mean any displacement vector u ¼ ðu1 ; u2 ; u3 Þ :

X ! R3 satisfying the following system of differential equations and boundary conditions (here and in what follows the summation convention will be used):

@ rij ðuÞ þ F i ¼ 0 in X; i ¼ 1; 2; 3; @xj

ð2:1Þ

where

  1 @uk @ul ; þ 2 @xl @xk

kF

1

ðx; 0ÞT t ðuÞðxÞk 6 T m ðuÞðxÞ;

ut ðxÞ–0 )

9 > > > > x 2 Cc ; > > =

F ðx;kut ðxÞkÞut ðxÞ ; ¼ T m ðuÞðxÞ kF ðx;kut ðxÞkÞut ðxÞk

x 2 Cc ;

where T t ðuÞ ¼ ðT t1 ðuÞ; T t2 ðuÞÞ; T ti ðuÞ :¼ rjk ðuÞmk t i;j , is the tangential contact stress. To give the weak and fixed-point formulation of our problem we shall need the following function spaces and sets:

K ¼ fv 2 V j v m 6 0 on Cc g; X m ¼ fu 2 L2 ðCc Þ j 9 v 2 V : u ¼ v m

on Cc g;

2

X tþ ¼ fu 2 L ðCc Þ j 9 v 2 V : u ¼ kv t k on Cc g: Further, let X 0m be the dual space to Xm with the duality pairing denoted by h , im and Km be the cone of non-negative elements of X 0m :

Km ¼ flm 2 X 0m j hlm ; v m im 6 0 8 v 2 Kg: The weak formulation of the contact problem with orthotropic Coulomb friction and the solution-dependent matrix of friction coefficients reads as follows:

9 > =

Find u 2 K such that aðu; v  uÞ þ jðkut k; T m ðuÞ; v t Þ

ð2:2Þ

rij mj ¼ Pi on CP ; i ¼ 1; 2; 3;

ð2:3Þ

ðPÞ

Z aðv ; wÞ ¼ cijkl eij ðv Þekl ðwÞ dx; v ; w 2 V; Z X Z F  v dx þ P  v ds; v 2 V; ‘ðv Þ ¼ CP

2

3

2

g 2 Km ; u 2 X tþ ; v 2 V; 3

with F 2 (L (X)) , P 2 (L (CP)) . Problem ðPÞ is an implicit variational inequality of elliptic type. One of possible ways how to prove the existence of a solution to ðPÞ is to use the fixed-point approach. Let (u, g) 2 X tþ  Km be given and define the auxiliary problem

Find u :¼ uðu; gÞ 2 K such that aðu; v  uÞ þ jðu; g; v t Þ  jðu; g; ut Þ

(static boundary conditions)

> ;

where

X

(kinematical boundary conditions)

u ¼ 0 on Cu ;

ð2:5Þ

> > > > > > ;

F 1 ðx; kut ðxÞkÞT t ðuÞðxÞ

jðu; g; v t Þ ¼ hg; kF ðuÞv t kim ;

i,j,k,l, = 1,2,3, is given by linear Hooke’s law and cijkl are elements of the bounded 4th order elasticity tensor C satisfying the usual symmetry and ellipticity conditions,

P ‘ðv  uÞ 8 v 2 K:

9 > = > ;

ðPðu; gÞÞ

This is a contact problem with orthotropic friction of Tresca type and the fixed matrix of friction coefficients F ðuÞ. Let W: X tþ  Km ? X tþ  Km be defined by

where m = (m1,m2,m3) is the outward unit normal vector to @ X, (non-penetration conditions)

T m ðuÞ :¼ rij ðuÞmi mj 6 0;

T m ðuÞum ¼ 0 on Cc :

ut ðxÞ ¼ 0 )

jðkut k; T m ðuÞ; ut Þ P ‘ðv  uÞ 8 v 2 K;

(equilibrium equations)

um :¼ u  m 6 0;

(orthotropic Coulomb friction law)

V ¼ fv 2 ðH1 ðXÞÞ3 j v ¼ 0 on Cu g;

X  R3 whose sufficiently smooth boundary @ X is decomposed into non-empty, non-overlapping parts Cu, CP and Cc. On Cu the body is fixed, surface tractions of density P = (P1, P2, P3) act on CP. Along Cc the body is unilaterally supported by a rigid half-space S. For simplicity of our presentation we shall suppose that Cc is flat and there is no gap between X and S. Finally, body forces of density F = (F1, F2, F3) are applied to X. Our aim is to find an equilibrium state of X taking into account effects of friction between X and S. By a classical solution of this

rij ðuÞ ¼ cijkl ekl ðuÞ; ekl ðuÞ ¼

the existence of the inverse F 1 we suppose that F 1 and F 2 do not vanish. Orthotropic Coulomb friction law with the solution-dependent matrix of friction coefficients reads as follows:

Wðu; gÞ ¼ ðkut k; T m ðuÞÞ; ð2:4Þ

To formulate orthotropic Coulomb friction law, we introduce principal orthotropic axes t1 = {t1,1,t1,2,t1,3}, t2 = {t2,1,t2,2,t2,3} on the tangent plane at a point on Cc so that {t1,t2,m} forms a local orthonormal basis. Let F 1 ; F 2 be the respective coefficients of friction in the directions t1 and t2, respectively. By ut 2 R2 we denote the vector whose components are the coordinates of u with respect to t1 and t2 on Cc, i.e. ut ¼ ðut1 ; ut2 Þ; uti ¼ u  t i . In what follows we shall suppose that both F 1 and F 2 may depend also on the Euclidean norm of ut on Cc, i.e. F i :¼ F i ðx; kut ðxÞkÞ; x 2 Cc ; i ¼ 1; 2. Finally let F ¼ diagðF 1 ; F 2 Þ be the (2  2) diagonal matrix. To guarantee

ðu; gÞ 2 X tþ  Km ;

ð2:6Þ

where u solves ðPðu; gÞÞ. Then u 2 K solves ðPÞ if and only if (kutk,  Tm(u)) is a fixed point of W. For discretization purposes we introduce the equivalent definition of the mapping W which uses the Lagrange multipliers associated with the unilateral constraint u 2 K. Instead of ðPðu; gÞÞ we consider the following problem:

Find ðu; km Þ 2 V  Km such that aðu; v  uÞ þ jðu; g; v t Þ  jðu; g; ut Þ P ‘ðv  uÞ  hkm ; v m  um im hlm  km ; um im 6 0 8 lm 2 Km :

9 > > > =

8 v 2 V; > > > ;

ðMðu; gÞÞ

151

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

It is well-known that ðMðu; gÞÞ has a unique solution (u,km) := (u(u,g), km(u,g)). In addition, u solves ðPðu; gÞÞ and km = Tm(u) on Cc. Then W defined by (2.6) takes the form

ð2:7Þ

3. Discretization of contact problems with orthotropic Coulomb friction. Existence and uniqueness results This section deals with a discretization of ðPÞ and the analysis of the resulting discrete problem. Unlike the continuous setting we prove that the discrete counterpart has a solution for any F ¼ diagðF 1 ; F 2 Þ satisfying (3.7). Moreover, it has a unique solution under the additional assumptions on F 1 and F 2 . Furthermore, we shall investigate how the uniqueness result depends on the size of the discrete problem. The discretization of ðPÞ will be based on the fixed-point formulation for an appropriate discretization of the mapping W given by (2.7). To get it we use a mixed finite element discretization of ðMðu; gÞÞ. Cc Let T X h ; T H be a partition of X and Cc into finite elements T and Cc R, respectively. The norm of T X h ; T H is denoted by h and H, respectively. We use the symbol H to point out that the partition of Cc is Cc X generally independent of T X is not h , but the case when T H ¼ T h j Cc

Cc excluded. With any T X h and T H we associate finite dimensional h H spaces V and L consisting of piecewise polynomial functions of degree less or equal k and s, respectively:

V h ¼ fv h 2 CðXÞ j v hjT 2 Pk ðTÞ 8 T 2 T Xh ; v h ¼ 0 on Cu g;

H

H

V hjC ; c

W hþ

h

H

8 v h 2 V h:

kv ht k1;Cc 6 ct kv h k1;Cc

h

¼ fu 2 W j u P 0 on Cc g:

Clearly, V and KHm serve as the natural discretizations of V and Km, respectively (observe that LH may consist also of discontinuous functions). Next we shall suppose that the couple (V h, LH) which will be used in the discretization is chosen in such a way that the following condition is satisfied: h

ðlH 2 LH & ðlH ; v hm Þ0;Cc ¼ 0 8 v h 2 V h Þ )

lH ¼ 0:

ð3:1Þ

Remark 3.1. Let us comment in short on the choice of Vh and LH satisfying (3.1). If LH: = Lh = Wh, i.e. the space LH is formed by the restrictions of functions from V h on Cc , then (3.1) is automatically satisfied. If V h, LH is the space of piecewise linear, piecewise Cc constant functions over T X h and T H , respectively, then (3.1) is satisfied provided that the ratio H/h is sufficiently large, i.e. the X c partition T C H is coarser than the partition of Cc given by T h jC (see c [12]). H h h H Let u 2 W þ and g 2 Km be given and define the problem:

9 Find ðu ; km Þ 2 V  Km such that > > > > > > aðuh ; v h  uh Þ þ jðuh ; g H ; v ht Þ > = jðuh ; g H ; uht Þ P ‘ðv h  uh Þ ðMhH ðuh ; g H ÞÞ > h > H h h h > ðkm ; v m  um Þ0;Cc 8 v 2 V ; > > > > ; H H H h H ðlm  km ; um Þ0;Cc 6 0 8 lm 2 Km ; h

H

h

ð3:3Þ

Since the function kv ht k does not belong to W hþ we have to introduce ‘‘a ‘‘return” operator” rh:H1(Cc) ? Wh possessing the following approximation and monotonicity property:

8 u 2 H1 ðCc Þ;

ku  r h uk0;Cc 6 cr hCc kuk1;Cc 1

u 2 H ðCc Þ; u P 0 on Cc ) rh u 2

ð3:4Þ

W hþ ;

ð3:5Þ

where cr is a positive constant which does not depend on the norm X hCc of T X belongs to a regular family of partih j . If k = 1 and T h j Cc

Cc

tions of Cc then (3.4) and (3.5) are satisfied by the Clément interpolation operator ([13]). We shall need also the satisfaction of the following inverse inequality for elements of V hjC : there exists a poc ð1;0Þ sitive constant cinv which does not depend on hCc such that ð1;0Þ 1

kv h k1;Cc 6 cinv hCc kv h k0;Cc

8 v h 2 V h:

ð3:6Þ

Let us recall that (3.6) is satisfied provided that T X hj

belongs to a

Cc

regular family of partitions of Cc satisfying the so-called inverse assumption (see [14, Th. 3.2.6]). Finally, we shall suppose that the coefficients of friction F 1 and F 2 are continuous and bounded:

9 > =

ð3:7Þ

>

; R1þ ;

where F min ; F max are given positive numbers. Let WhH : W hþ  KHm ! W hþ  KHm be defined by

H

h

h

and there exists a constant ct > 0 independent of v 2 V and h such that

i ¼ 1; 2; 8 ðx; nÞ 2 Cc  h 3

V ¼ ðV Þ ; Km ¼ fl 2 L j l P 0 on Cc g; W ¼

h

0 < F min 6 F i ðx; nÞ 6 F max ;

and set

ð3:2Þ

2

F 1 ; F 2 2 CðCc  R1þ Þ;

LH ¼ flH 2 L2 ðCc Þ j lHjR 2 Ps ðRÞ 8 R 2 T CHc g

h

v ht ¼ ðv ht ; v ht Þ 2 ðH1 ðCc ÞÞ2 8 v h 2 V h 1

Wðu; gÞ ¼ ðkut k; km Þ:

h

x ´ (t1(x), t2(x)) associating any x 2 Cc with the principal orthotropic axes is sufficiently smooth so that

H

where the meaning of all symbols has been introduced in Section 2. Because of (3.1), this problem has a unique solution ðuh ; kHm Þ :¼ ðuh ðuh ; g H Þ; kHm ðuh ; g H ÞÞ for any ðuh ; g H Þ 2 W hþ  KHm . Before we give the definition of the discrete contact problem with orthotropic Coulomb friction we formulate assumptions needed in what follows. First we shall suppose that the vector field

WhH ðuh ; g H Þ ¼ ðrh kuht k; kHm Þ; ðuh ; g H Þ 2 W hþ  KHm ; h

H

ð3:8Þ h

H

where ðu ; km Þ is the unique solution of ðMhH ðu ; g ÞÞ. This definition is meaningful since (3.2) implies that kuht k 2 H1 ðCc Þ. The mapping WhH serves as the approximation of W. Analogously to the continuous setting we say that uh 2 Vh is a solution of the discrete contact problem with orthotropic Coulomb friction and the solution-dependent coefficients of friction iff the couple ðr h kuht k; kHm Þ is a fixed point of WhH. Next we shall examine the existence, eventually the uniqueness of the fixed points of WhH. To this end we introduce the following norm in Wh  LH:

kðuh ; lH ÞkW h LH ¼ kuh k0;Cc þ klH k1=2;h ; ðuh ; lH Þ 2 W h  LH ; ð3:9Þ where

klH k1=2;h ¼ sup

0–v h 2V h

ðlH ; v hm Þ0;Cc : kv h k1;X

ð3:10Þ

In view of (3.1), the formula (3.10) defines a mesh-dependent norm in LH. To prove the existence of a fixed point we use Brouwer’s fixedpoint theorem. Since the mapping WhH is continuous in virtue of (3.7), it remains to show that it maps a closed convex set into itself. This property follows from the next lemma. Lemma 3.1. Let (3.1)–(3.6) be satisfied. Then there exist positive constants R1, R2 which do not depend on h and H such that WhH maps W hþ  KH m \ B into itself, where B is the ball

B ¼ fðuh ; lH Þ 2 W h  LH j kuh k0;Cc 6 R1 ; klH k1=2;h 6 R2 g:

152

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

Proof. Inserting vh = 0, 2uh into the first inequality in ðMhH ðuh ; g H ÞÞ and using the V-ellipticity of a with the constant of ellipticity a > 0 we easily obtain that uh is bounded:

e :¼ RðkFk e 9R 0;X ; kPk0;CP Þ > 0 :

e akuh k1;X 6 R:

ð3:11Þ

Furthermore,

  kr h kuht kk0;Cc  r h kuht k  kuht k0;Cc þ kkuht kk0;Cc ð3:4Þ

 cr hCc kuht k1;Cc þ kuht k0;Cc

ð3:3Þ

 cr ct hCc kuh k1;Cc þ kuh k0;Cc

ð3:6Þ

ð1;0Þ

 ðcr ct cinv þ 1Þctr kuh k1;X ;

where ctr is the norm of the trace mapping from (H1(X))3 into (L2(Cc))3. From this and (3.11) the existence of R1 :¼ R1 ðkFk0;X ; kPk0;CP ; aÞ > 0 bounding kr h kuht kk0;Cc follows. Restricting ourselves to functions from Vh such that v ht ¼ 0 on Cc we easily get:

  e6R e kak þ 1 :¼ R2 ; kkHm k1=2;h 6 kakkuh k1;X þ R

a

e is from (3.11) and kak is the norm of a. where R

4. Existence of local Lipschitz continuous branches of solutions h

Remark 3.2. It is worth mentioning that R1 and R2 are independent of F . A natural question arises, namely under which conditions the fixed point of WhH is unique. To this end let us suppose that the friction coefficients F 1 , F 2 are smooth enough and denote

(

 ) @F i ðx; nÞ   L ¼ max sup  i¼1;2 @n  x2Cc ; n>0 and

jðF Þ ¼ sup x2Cc ; n>0

maxfF 1 ðx; nÞ; F 2 ðx; nÞg : minfF 1 ðx; nÞ; F 2 ðx; nÞg

The quantity jðF Þ can be considered as a ‘‘measure” of orthotropy (for the isotropic friction jðF Þ ¼ 1). It can be shown (see [15]) that WhH is Lipschitz continuous in W hþ  KHm \ B, i.e. there exists C > 0 such that

 h ; gH ÞkW h LH 6 Ckðuh ; g H Þ  ðu

ð3:12Þ

 h ; gH Þ 2 W hþ  KHm \ B, where B is the holds for every ðuh ; g H Þ; ðu same as in Lemma 3.1. Moreover, the Lipschitz constant C in (3.12) is of the form C ¼ maxfC 1 ðF max ; HÞ; C 2 ðL; jðF Þ; H; hCc Þg, where the constants C1, C2 enjoy the following properties: (a) C 1 ðF max ; HÞ ! 0 if F max ! 0þ for any H > 0 fixed, C 2 ðL; jðF Þ; H; hCc Þ ! 0 if L ? 0+ for any H; hCc fixed and jðF Þ bounded; (b) if F does not depend on kuht k, i.e. L  0, then C 2 ð0; jðF Þ; H; hCc Þ ¼ 0 for any F and any H; hCc > 0; (c) if F is fixed then C 1 ðF max ; HÞ; C 2 ðL; jðF Þ; H; hCc Þ behaves as H1/2 and ðHhCc Þ1=2 , respectively, for H; hCc ! 0þ provided that the Babuška-Brezzi condition for {Vh, LH} is satisfied:

0–v

h 2V h

ðlH ; v hm Þ0;Cc P bklH k;Cc kv h k1;X

8 lH 2 LH ;

where b > 0 does not depend on h and H and

klH k;Cc

hl ; v m im ¼ sup : 0–v 2V kv k1;X H

In this section we shall consider solutions to discrete contact problems to be a function of friction coefficients. We restrict ourselves to 2D-case and isotropic Coulomb friction with the coefficient F which depends on the spatial variable x 2 Cc. From the previous section we know that if F 6 F crit then there exists a unique solution of our problem. Moreover, it can be shown ([6]) that for such F the mapping uh : F #uh ðF Þ is Lipschitz continuous if (3.12) is satisfied. In what follows we shall investigate the case when F > F crit . More precisely, if F 0 > F crit is a reference point, under which conditions there exists a neighborhood U d ðF 0 Þ such that the mapping uh : F #uh ðF Þ; F 2 U d ðF 0 Þ, has a Lipschitz continuous branch, implying (among others) local uniqueness of the respective solution. To this end we use the algebraic formulation of the discrete contact problem with Coulomb friction which involves two Lagrange multipliers: one releasing the unilateral constraint and the other regularizing the non-smooth frictional term. It reads as follows:

9 Find ðu; km ; kt Þ 2 Rn  Km  Kt ðkm Þ such that > > > T T = Au ¼ f  B k  B Fk ; m m

t

t

ðlm  km Þ  Bm u þ Fðlt  kt Þ  Bt u 6 0

 h ; gH ÞkW h LH kWhH ðuh ; g H Þ  WhH ðu

sup

Suppose that jðF Þ is bounded. From (a) it follows that for any H; hCc fixed one can find F crit :¼ F crit ðHÞ > 0 and Lcrit :¼ Lcrit ðH; hCc Þ > 0 such that if F max 6 F crit and L 6 Lcrit then the mapping WhH is contractive. Therefore there exists a unique fixed point of WhH and the method of successive approximations converges. Recall that each iterative step of the method of successive approximations is represented by a contact problem with Tresca friction and in the case of solution-dependent friction coefficients we update not only the slip threshold but also the matrix F . If F does not depend on the solution, i.e. L  0, then using (b) and (c) we recover the classical result (see [11]). If both H and hCc tend to zero then to preserve the contractivity of WhH, the parameters F max ; L have to decay at least as fast as H1/2 and ðHhCc Þ1=2 , respectively. This is a consequence of (c). The mesh dependency of discrete models can be read in two ways: either (i) the matrix F is fixed, then passing from coarser to finer meshes we may loose unicity of the approximate solution or (ii) finite element meshes are fixed, then setting F f ¼ fF ; f P 0, one can find fcrit > 0 such that the discrete model has a unique solution for f 6 fcrit and eventually multiple solutions if f > fcrit. This behavior has been observed in computations (see [9]).

8 ðlm ; lt Þ 2 Km  Kt ðkm Þ;

ð4:1Þ

> > > ;

where A is an (n  n) stiffness matrix, u 2 Rn is the nodal displacement vector, Bm, Bt are (p  n) matrices representing the linear mappings u ´ um and ut, respectively, and p is the number of the contact nodes. We shall suppose that

BTm lm þ BTt lt ¼ 0 () ðlm ; lt Þ ¼ ð0; 0Þ 2 Rp  Rp : Further, f 2 Rn is the load vector, F ¼ diagðF 1 ; . . . ; F p Þ is the diagonal matrix with F i being the value of F at the ithcontact node, Km ¼ Rpþ and

Kt ðgÞ ¼ fl 2 Rp j jli j 6 g i 8 i ¼ 1; . . . ; pg;

g 2 Km :

It is well-known that km and Fkt is the opposite of the discrete normal and tangential contact stress, respectively. Let f 2 Rn be fixed. Using tools of convex analysis, problem (4.1) can be equivalently written as the following system of generalized equations:

Find y 2 Rnþ2p such that 0 2 C f ðF ; yÞ þ Q ðyÞ; Rpþþ

nþ2p

nþ2p

nþ2p

ð4:2Þ nþ2p

where C f : R !R and Q : R R is the singlevalued and the set-valued mapping, respectively, defined by

153

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

0

A B C f ðF ; yÞ ¼ @ Bm

10

0 1 u f BTt F CB C B C 0 A@ km A  @ 0 A;

BTm 0

FBt 0 0 1 0 B C Q ðyÞ ¼ @ N Km ðkm Þ A; N Kt ðkm Þ ðkt Þ F :¼ ðF 1 ; . . . ; F p Þ 2

0

Rpþþ ;

1

kt

S F ðgf Þ ¼ fðu; km ; kt Þg;

where {(u,km,kt)} denotes the set of all solutions to the discrete contact problem with Coulomb friction with the coefficient F and the load vector gf, we get the following result.

0

y :¼ ðu; km ; kt Þ 2 R

nþ2p

;

with N Km ðlÞ; N Kt ðkm Þ ðlÞ denoting the normal cones of Km and Kt(km) at l 2 Rp , respectively. The system (4.2), in which F plays the role of the perturbation parameter, can be analyzed by using the generalized implicit-function theorem ([16]). Let S f : Rpþþ Rnþ2p be the solution (multi-valued) mapping of (4.2):

S f ðF Þ ¼ fy 2 Rnþ2p j 0 2 C f ðF ; yÞ þ Q ðyÞg;

F 2 Rpþþ ;

and ðF 0 ; y0 Þ 2 Rpþþ  Rnþ2p be a reference point such that y0 2 S f ðF 0 Þ. Theorem 2.1 in [16] states that S f has a local Lipschitz continuous branch passing through y0 at F 0 if the so-called strong regularity condition is satisfied at ðF 0 ; y0 Þ. This condition pertains to the properties of the multi-valued mapping Rf : Rnþ2p Rnþ2p , where

Rf ðgÞ ¼ fy 2 Rnþ2p j

g 2 C f ðF 0 ; y0 Þ þ $y C f ðF 0 ; y0 Þðy  y0 Þ þ Q ðyÞg; g 2 Rnþ2p ; i.e. Rf (g) is the solution set of the system of generalized equations arising from the partial linearization of the smooth part Cf in (4.2) at ðF 0 ; y0 Þ with respect to the second variable. Taking g :¼ ðgu ; gm ; gt Þ 2 Rnþ2p ; Rf ðgÞ consists of all y satisfying

9 0 ¼ Au þ BTm km þ BTt F 0 kt  f  gu ; > = 0 2 Bm u  gm þ N Km ðkm Þ;

0 2 F 0 Bt u  gt þ N Kt ðkm Þ ðkt Þ;

> ;

ð4:3Þ

with F 0 ¼ diagðF 01 ; . . . ; F 0p Þ. Setting

^ :¼ u þ u



Bm 0

F Bt

þ 



gm ; gt

 þ Bm where stands for the Moore-Penrose pseudo-inverse of 0 F Bt   Bm , we have F 0 Bt



        Bm Bm Bm þ gm ^ u ¼ u þ ¼ gt F 0 Bt F 0 Bt F 0 Bt F 0 Bt Bm

Bm u þ gm

!

F 0 Bt u þ gt

9 ^ þ B m km þ 0 ¼ Au > > >  þ   > > gm Bm >  gu ; = f  A 0 gt F Bt > > > ^ þ N Km ðkm Þ; > 0 2 Bm u > > ; 0 ^ þ N Kt ðkm Þ ðkt Þ: 0 2 F Bt u

Theorem 4.1. Let us suppose that S F 0 has a local Lipschitz continuous branch containing y0 in a vicinity of f 2 Rn , i.e. there exist: a singlevalued Lipschitz continuous function /F 0 from a neighborhood O of f b of y0 such that into Rnþ2p and a neighborhood Y

BTt F 0 kt

b & /F 0 ðgf Þ ¼ S F 0 ðgf Þ \ Y

/F 0 ðf Þ ¼ y0

8 gf 2 O:

Then there are neighborhoods U, Y of F and y0, respectively, and a single-valued Lipschitz continuous function rf: U ? Y satisfying 0

rf ðF 0 Þ ¼ y0 & rf ðF Þ ¼ S f ðF Þ \ Y

8 F 2 U:

This theorem says that locally the dependence of a solution on F can be deduced from the dependence of the solution on the load vector keeping F fixed. It is worth mentioning that the latter one is much simpler since the dependence on the load vector is piecewise linear. The local behavior of the set-valued mapping f #S F ðf Þ; f 2 Rn , for F 2 Rpþþ fixed can be studied in more details by writing (4.1) as a system of piecewise smooth equations (see (5.1)) and applying the implicit-function theorem for such a class of functions (Theorem 4.2.2 in [17]). From it one concludes that the existence of local Lipschitz continuous branches is not guaranteed only around points y 2 S F ðf Þ for f belonging to the union of subspaces of dimension strictly lower than n. (For more details see [6].)

5. Numerical continuation of solution curves In this section we shall propose the algorithm for tracing the (local) Lipschitz continuous branches of solutions, the existence of which was established in the previous section. More precisely, taking a smooth path a 2 I #F ðaÞ ¼ ðF 1 ðaÞ; . . . ; F p ðaÞÞ 2 Rpþ ; I  R1 open, we shall approximate the curve consisting of all solutions to the discrete contact problem with Coulomb friction and the coefficient F ðaÞ for a running over I. Next we use the equivalent formulation of (4.1) as the system of non-smooth equations involving the projection mappings P Km : Rp ! Km and P Kt ðF km Þ : Rp ! Kt ðF km Þ, where Km ¼ Rpþ and

Kt ðF gÞ ¼ fl 2 Rp j jli j 6 F i g i 8 i ¼ 1; . . . ; pg;

g 2 Km :

Taking F :¼ F ðaÞ; a 2 I, this system becomes:

Find x 2 Rnþ2p  I such that HðxÞ ¼ 0; where H : R

nþ2p

0 B HðxÞ ¼ @

and (4.3) transforms into T

gf 2 Rn ;

nþ2p

I !R T

Au þ Bm km þ

ð5:1Þ

is defined by

BTt kt

f

km  P Km ðkm þ rBm uÞ

1 C A;

x :¼ ðu; km ; kt ; aÞ 2 Rnþ2p  I:

kt  P Kt ðF ðaÞkm Þ ðkt þ rBt uÞ ð4:4Þ

^ ; km ; kt Þ Comparing this with (4.2), it is readily seen that the triplet ðu satisfies (4.4) iff it is a solution to the discrete contact problem with Coulomb friction with the coefficient F 0 and the perturbed load vector

 þ   gm B þ gu : gf :¼ f þ A 0 m gt F Bt Hence, introducing the set-valued mapping S F : Rn Rn  Rp  Rp for F 2 Rpþþ fixed by

Here r > 0 is a parameter and the components of P Km ; P Kt ðF ðaÞkm Þ are given by

ðP Km Þi ðlÞ ¼ P½0;þ1Þ ðli Þ; ( P½F i ðaÞkm;i ;F i ðaÞkm;i ðli Þ if km;i P 0; ðP Kt ðF ðaÞkm Þ Þi ðlÞ ¼ P ½F i ðaÞkm;i ;F i ðaÞkm;i ðli Þ if km;i < 0; i ¼ 1; . . . ; p; l 2 Rp ; with P[0,+1), P[a,b] being the projections of R1 onto [0, + 1) and [a,b],  1 < a 6 b < +1, respectively. The meaning of the remaining symbols in the definition of H is the same as in (4.1). Note that H is a piecewise smooth function, i.e. for every  2 Rnþ2p  I there exists an open neighborhood O  Rnþ2p  I; x

154

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

 2 O, x and a finite number of smooth functions HðiÞ : O ! Rnþ2p ; i ¼ 1; . . . l, (the so-called selection functions) such that HðxÞ 2 fHð1Þ ðxÞ; . . . ; HðlÞ ðxÞg for every x 2 O. The set

Þ  fi 2 f1; . . . ; lg j Hðx Þ ¼ HðiÞ ðx Þg IH ðx  and the functions HðiÞ ; i 2 IH ðx Þ, is known as the active index set at x . are termed the active selection functions for H at x Classical continuation techniques require H in (5.1) to be smooth. Next, we shall show how these techniques can be adapted to our non-smooth case. In particular, we shall modify the MoorePenrose continuation, which is presented e.g. in [18]. This continuation procedure consists of two steps. In the correction, we use the piecewise smooth Newton method (7.2.14 Algorithm in [19]) instead of the smooth one when the gradient $H is replaced by the gradient of one of its active selection functions if necessary. On the other hand, a generalization of the prediction is more complicated. Indeed, if one tries to use the tangential direction of the form

H0 ðx j ; s j Þ ¼ 0;

ks j k ¼ 1;

j ¼ 1; 2; . . . ;

where H0 ðx j ; s j Þ denotes the directional derivative of H at x j in the direction s j, the continuation may fail when approaching a point of non-differentiability on the solution curve. It is caused by the fact that the Newton corrections are only locally convergent and one has to take a suitable tangential vector to reach their zone of convergence (see Fig. 1 for illustration). This is why we have to develop a special approach to pass through such points. Clearly, the non-differentiability of H is caused by the functions:



x #km  P Km ðkm þ rBm uÞ;

ð5:2Þ

x #kt  P Kt ðF ðaÞkm Þ ðkt þ rBt uÞ:

x #  rðBm uÞi

and

x #kt;i þ F i ðaÞkm;i ;

x #  rðBt uÞi ;

x #kt;i  F i ðaÞkm;i ;

x #ð2kt þ rBt uÞi ;

i = 1,. . .,p, respectively, representing all possible values of the mappings in (5.2). Furthermore, we define the so-called test functions hk ¼ ðhk1 ; . . . ; hkp Þ : Rnþ2p  I ! Rp ; k ¼ 1; 2; 3, by

h1i ðxÞ ¼ ðkm þ rBm uÞi ;

$Hði2 Þ ðx j Þs j ¼ 0;

h2i ðxÞ ¼ ðkt þ rBt uÞi  F i ðaÞkm;i ;

h3i ðxÞ ¼ ðkt þ rBt uÞi þ F i ðaÞkm;i ; i ¼ 1; . . . ; p; x 2 Rnþ2p  I: From their definition it immediately follows that there is the one-toone correspondence between the signs of the components of h1(x), h2(x) and h3(x) and the selection functions for H which are active at x. (Possible zero components indicate that more than one selection function is active.) Suppose that x 2 Rnþ2p  I is a point where only

ks j k ¼ 1:

To select the direction of this vector, the notion of orientation is adapted from the theory of smooth continuations. Definition 5.1. Let H be smooth at a point x 2 Rnþ2p  I. The tangential vector s 2 Rnþ2pþ1 satisfying

$HðxÞs ¼ 0;

ksk ¼ 1;

is called positively oriented iff



With (5.2) the following selection functions will be associated:

x #km;i ;

one selection function is active, i.e. all components of the test functions are nonzero there. Assembling the signs of the test functions into a (p  3) array in such a way that the ith column corresponds to hi(x), i = 1,2,3, we see that every selection function for H can be represented by a (p  3) array and this representation is unique.  of non-difLet x j be a current point which is close to a point x ferentiability of H as illustrated in Fig. 1. Assume that exactly  and there extwo selection functions Hði1 Þ and Hði2 Þ are active at x  which ists a piecewise smooth curve of solutions passing through x consists of two smooth branches belonging to the solution sets to Hði1 Þ ðxÞ ¼ 0 and Hði2 Þ ðxÞ ¼ 0. Supposing that x j is a root of Hði1 Þ , we shall describe how to reach the unknown smooth branch of the curve corresponding to Hði2 Þ ðxÞ ¼ 0. Clearly, one of the test functions, say hk, has a zero component , say the mth one, and this component changes its sign when at x . Moreover, continuity of hk ensures that hkm ðx j Þ passing through x is close to zero. If the (p  3) array represents Hði1 Þ then changing the sign of hkm we obtain the representative of the selection function Hði2 Þ . This leads us to the following choice of the tangential vector s j at x j:

det

$HðxÞ

 > 0:

sT

In the opposite case it is called negatively oriented. In our non-smooth case, we determine the direction of s j to preserve the orientation in the following sense:

det

$Hði1 Þ ðx j1 Þ

!

ðs j1 ÞT

det

x

xj x j−1

Fig. 1. Necessity of a good tangential prediction.

! > 0:

ðs j ÞT

Let us note that the expounded procedure can be also applied when  is met exactly, i.e. x j ¼ x . Neverthe point of non-differentiability x theless, this situation is highly improbable. On the basis of the above considerations we propose the following algorithm. Algorithm 5.1 (Piecewise smooth variant of the Moore-Penrose continuation). 0

Data: e,e > 0, hmax P hinit P hmin > 0, hinc > 1 > hdec > 0, kmax P kthr > 0 and x0 2 Rnþ2p  I; s0 2 Rnþ2pþ1 satisfying:

H0 ðx0 ; s0 Þ ¼ 0;

kHðx0 Þk < e; h jτ j

$Hði2 Þ ðx j Þ

ks0 k ¼ 1:

Step 1: Set h0: = hinit, j: = 0. Step 2: Set mdec: = 0. Step 3 (prediction): Set X0: = xj + hjsj, T0: = sj, k: = 0. k Step 4 (correction): Select ! an index  ik in IH ðX Þ and set: !

B :¼

$Hðik Þ ðX k Þ ðT k ÞT

e :¼ B1 R; T

;

T kþ1 :¼

X kþ1 :¼ X k  B1 Q :

R :¼ e T ; ek kT

0 1

;

Q :¼

HðX k Þ 0

;

155

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161 0

Step 5: If kHðX kþ1 Þk < e and kXk+1  Xkk < e , set xj+1: = Xk+1, sj+1: = Tk+1 and go to Step 9. Step 6: If k < kmax, set k: = k + 1 and go to Step 4. Step 7: If hj > hmin, set hj: = max{hdechj,hmin}, mdec: = mdec + 1 and go to Step 3. Step 8: According to a component of h1(xj), h2(xj) or h3(xj) close to 0 select a function HðiÞ which is likely to be active in a vicinity of xj and compute a vector sj satisfying

$HðiÞ ðxj Þsj ¼ 0;

ksj k ¼ 1

In the forthcoming computations, we will use the following parameter setting: e = e0 = 106, hmin = 105, 0.1 6 hmax 6 5, hinit = 0.05, hinc = 1.3, hdec = 0.5, kmax = 10, kthr = 4. 6.1. One contact node This subsection deals with an elementary problem involving a single linear triangular finite element shown in Fig. 2. This problem was suggested and analyzed in [10]: let f = (fm, ft) be the load vector, x = (um, ut,km,kt,a) be the solution of the piecewise smooth system

and preserving the orientation. Set hj: = hinit and go to Step 2. Step 9: Set

hjþ1 :¼



minfhinc hj ; hmax g if k < kthr and mdec ¼ 0 hj

otherwise

and j: = j + 1; go to Step 2. 0

Here e and e are convergence tolerances, hmin, hmax and hinit is the minimal, maximal and initial step length, respectively, and hinc, hdec are the scale factors for adjustment of the step length. Further, kmax stands for the maximal number of corrections allowed and mdec denotes the number of the step length reductions of hj. In Step 7 the current step length is shortened in case of nonconvergence of the corrections. Step 9 defines the step length for the prediction in the next iteration. The new step length hj+1 can be larger than hj only if the number of corrections (Step 4) does not exceed kthr given a priori and mdec = 0. These parts of the routine together with the prediction and the corrections are taken directly from the classical Moore-Penrose continuation. Step 8 is added to handle the situation when the corrections do not converge even for h = hmin. Then we determine a new tangential prediction, making use of our test functions, and return to the classical part of the procedure.

Fig. 3. Geometry of the second elementary problem.

1.6 1.4 1.2

branch 2

1 0.8 0.6 0.4 0.2

Remark 5.1

0

(i) This algorithm can be also used to pass through a point where more than two selection functions for H are active. In this case, however, one has more possibilities to choose a new selection function when ‘‘switching” between smooth branches. (ii) Control of changes of the signs of components of h1, h2 and h3 during the continuation routine enables us to detect points of non-differentiability with an arbitrarily chosen accuracy.

branch 1 −0.2

0

10

20

30

40

50

0

10

20

30

40

50

1.6 1.4 1.2 1

6. Numerical examples 0.8

The last section is devoted to the application of Algorithm 5.1. We shall restrict ourselves to very simple model examples with one and two contact nodes (see Figs. 2 and 3). zero Dirichlet condition

0.6 0.4 0.2 0 −0.2

finite element rigid foundation Fig. 2. Geometry of the first elementary problem.

Fig. 4. Example 6.1. The line interpretations (at the top): solid (no contact), dashed (contact-stick), dash-dotted (contact-slip). At the bottom: An illustration of the adaptive step length refinement. The asterisks show the initial points of the continuation, whereas the arrows indicate the positively oriented tangents at the initial points.

156

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

(5.1) with n ¼ 2; p ¼ 1; F ðaÞ ¼ a > 0 and k,l > 0 be the Lamé coefficients. Then (j) if ((k + 3l)fm + (k + l)ft 6 0 & fm 6 0) _ ((k + 3l)fm + (k + l)ft > 0) then there exists just one solution branch; (jj) if ((k + 3l)fm + (k + l)ft < 0 & fm > 0) then there exist two separated solution branches; (jjj) if ((k + 3l)fm + (k + l)ft = 0 & fm > 0) then there exist two solution branches which meet one another in a bifurcation point. The respective solution branches are given explicitly in [10]. For comparison, we will compute them numerically. The following three examples pertain to the cases (jj), (j) and (jjj) with k = l = 1.

and in what follows by solid, dashed and dash-dotted curves, respectively. Along these curves, just one selection function is active. The points marked by the circle refer to the transition points between the contact modes: at such points at least two selection functions are active and Step 8 of Algorithm 5.1 has to be executed. 1.2

1

branch 2 0.8

0.6

Example 6.1. Data: f = (1.5,  4). There exist two disjoint solution branches (1 and 2), see Fig. 4 at the top. We always plot only the relation a ´ km, since the components um, ut and kt are uniquely determined by km. The performance of Algorithm 5.1 is illustrated at the bottom of Fig. 4: the continuation starts at two points on each branch marked by the asterisks. We follow the solution in both positive and negative directions (the arrows always mark the positive directions). The solutions obtained by Algorithm 5.1 are represented by the dots . The solution points are classified as no contact (meaning that um < 0), contact-stick ðkm > 0 & jkt j < F km ) and contact-slip points (km > 0 & ut – 0). The sets of such points will be denoted here

0.4

0.2

branch 1 0

−0.2

0

5

10

15

20

0

5

10

15

20

0

5

10

15

20

−6

5

x 10

5.5 4 5 3

4.5 4

2

3.5

1

3

0

2.5

−1

2 −2 1.5 −3 1 −4

0.5

−5

0 0

5

10

15

20 1.2

5.5 5

1

4.5 0.8

4 3.5

0.6

3 0.4

2.5 2

0.2 1.5 1

0

0.5 −0.2

0 0

5

10

Fig. 5. Example 6.2.

15

20 Fig. 6. Example 6.3. The dotted line at the top denotes solutions corresponding to grazing contact.

157

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

Example 6.2. Data: f = (1.5,7). The solution set consists of just one branch, see Fig. 5. Example 6.3. Data: f = (1,  2). The solution set is depicted at the top of Fig. 6: for a = (k + 3l)/(k + l) = 2 branch 1 bifurcates and the corresponding contact problem possesses continuum of solutions represented by the vertical segment. The actual computed paths when starting on branches 1 and 2 are depicted in Fig. 6 in the center and at the bottom, respectively. Although the whole branch 1 consists of the grazing contact points, i.e. um = km = 0, and two selection functions are active on it, the algorithm copes even with this situation. Observe that the direction of continuation during the transition between branches 1 and 2 is determined in such a way that the orientation is preserved. Nevertheless, if one changes the orientation in the corresponding transition point, the other part of branch 1 is revealed. It is worth mentioning that an arbitrarily small perturbation of f destroys the bifurcation. For example, if we set f = (0.98,  2) then we arrive at (jj): the solution set consists of two distinct branches, see Fig. 7. 6.2. Two contact nodes Next we shall consider a discrete contact problem with two contact nodes (see Fig. 3), i.e. n = 4, p = 2 in (5.1). We set

F ðaÞ ¼ ða; aÞ; a > 0; f ¼ ðf 1 ; f 2 Þ; f i ¼ ðfm;i ; ft;i Þ; i ¼ 1; 2, and k = l = 1. Unlike to the previous case, this time the exact solutions are not at our disposal. The following figures are constructed on the basis of the solutions computed by Algorithm 5.1. Example 6.4. Data: f1 = (1,  3), f2 = (1,  1). We found numerically three branches depicted in Figs. 8 and 9 (one has to be aware of different scaling of the figures). Again we plot only the relation a ´ km,1,km,2, i.e. the normal stresses in the first and the second contact node. Solutions of the system (5.1) can be deduced from the figures in the following way: for a fixed a they consist of pairs (km,1,km,2) whose components are the intersection of the vertical line at a with the same branch i associated with the first and the second contact node. Nevertheless, when coupling the components km,1 and km,2 from the same branch, one has to take into account also the order in which they were obtained during the continuation. Let us take a = 5, for example. Then (5.1) has (at least) five solutions — one on branch 1, two on branch 2 (the upper intersection in the first component with the upper intersection in the second component and the lower one with the lower one) and two on branch 3 (the lower intersection in the first component with the upper intersection in the second component and vice-versa). The performance of Algorithm 5.1 is illustrated in Fig. 10: the continuation is initialized at the points marked by the asterisk corresponding to the value a = 5 and each point is computed with a 0 guaranteed precision given by the tolerances e = e = 106.

1.2 1 1 0.8

0.8

branch 2

branch 2 0.6

0.6

0.4 0.4 0.2 0.2 0

branch 1

branch 1 −0.2

0

5

10

15

0 2

20

3

4

5

6

7

8

9

10

11

7

8

9

10

11

1.2 1 1 0.8

0.8

branch 2 0.6

branch 2

0.6

0.4 0.4 0.2 0.2 0

branch 1

branch 1 −0.2

0

5

10

15

Fig. 7. A small perturbation destroys the bifurcation.

0 20

2

3

4

5

6

Fig. 8. Example 6.4, branch 1 and 2: a ´ km,1 (at the top), a ´ km,2 (at the bottom).

158

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

The crucial problem is to find starting points x0 for the continuation routine. In general, we can fix a > 0 and compute a corresponding solution to (5.1) by the piecewise smooth Newton method, see e.g. [19].

Remark 6.1. In order to find all roots, we proceed as follows: let us fix a > 0. We determine all particular selection functions for H. For each selection function, which is linear, we solve (5.1). Finally, we sort out the classes of the roots. Obviously, such an algorithm has the exponential complexity.

0.6

Example 6.5. Data: f1 = (0.4,  1.4), f2 = (1.5,  1.3). We found three branches depicted in Figs. 11 and 12. It is worth mentioning that the points on branch 3 in the graph a ´ km,1 has to be viewed twice, i.e. there are two solutions with the same km,1 but different km,2. The numerical performance is illustrated in Fig. 13. Next, we give an example of a bifurcation, similar to that one shown in Fig. 6.

0.5

0.4

branch 3 0.3

0.2

0.1

0 2

3

4

5

6

7

8

9

10

11

Example 6.6. Data: f1 = (2,  3), f2 = (1,  1). In Fig. 14, we observe the bifurcation phenomenon: branch 1 bifurcates at a = 1.3333 (for better understanding, the transition points are numbered). Note that the abscissa between transition points 3 and 4 in the graph a ´ km,1 is vertical and it represents together with the corresponding singleton in the graph a ´ km,2 continuum of solutions of the respective contact problem. Furthermore, the whole branch 1 consists of points corresponding to grazing contact of the first contact node and no contact of the second contact node. The

0.2

1.6 0.15 1.4 1.2

0.1

1

branch 3

0.05

0.8 0.6

branch 2

0

0.4 0.2

−0.05

2

3

4

5

6

7

8

9

10

11

Fig. 9. Example 6.4, branch 3: a ´ km,1 (at the top), a ´ km,2 (at the bottom).

branch 1

0 −0.2

0

2

4

6

8

10

1.6 1 1.4 0.8

1.2

branch 2

branch 2

1 0.6 0.8 0.6

0.4

0.4 0.2 0.2

branch 3

branch 1

0

0 2

3

4

5

6

7

8

9

10

11

Fig. 10. Illustration of the adaptive step length control in Example 6.4: branch 2 and 3 for a ´ km,2.

−0.2

0

2

4

6

Fig. 11. Example 6.5, branch 1 and 2.

8

10

159

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

second contact node comes into contact in transition point 3 on branch 2. The output of the numerical continuation of branch 2 when starting from a = 5 is shown in Fig. 15.

1

2 1.8

2 1.6

3

0.6

1.4 1.2

0.5

branch 2

1 0.4

0.8 0.6

0.3

0.4 0.2

0.2

branch 1

4 0

branch 3

0.1

0

1

2

3

4

5

6

7

0

−0.1

1

1 0

2

4

6

8

10 0.8

0.6

2

0.5

0.6

0.4

0.4

branch 2

branch 3

0.3

0.2

0.2

branch 1

3, 4 0

0.1

0

0

−0.1

1

2

3

4

5

6

7

Fig. 14. Example 6.6. The abscissa connecting transition points 3 and 4 at the top corresponds to a single point at the bottom. The dotted line at the top denotes solutions corresponding to grazing contact.

0

2

4

6

8

10

Fig. 12. Example 6.5, branch 3. Continuing in the negative direction, a goes back and forth while km,1 remains the same.

Example 6.7. Data: F ¼ ð4; 4Þ; f ðaÞ ¼ ðf 1 ðaÞ; f 2 ðaÞÞ with f1(a) = (0.1a + 0.4, 1.1a + 0.2), f2(a) = (0.2a + 1.8,0.8a  0.1),  3 < a < 3. The linear loading path results in the piecewise linear solution branch (see Fig. 16). Again, it suffices to plot the relation a ´ km,1, km,2. Note that we encounter at most five solutions of the corresponding contact problem for a fixed. For example, there are three intersections of the branch with the vertical line at a = 1.6 at the first contact node and five intersections at the second contact node. However, any point between transition points 1 and 2 of the first contact node has to be counted three times as seen from Fig. 16.

1.6 1.4

Remark 6.2. It is worth mentioning that a simple adaptation of Algorithm 5.1 enables us to continue along a loading path a 2 I ´ f(a) keeping F fixed, see [20].

branch 2

1.2 1 0.8 0.6 0.4

7. Conclusions

branch 3 0.2 0

2

3

4

5

6

7

8

9

Fig. 13. Illustration of the numerical continuation in Example 6.5.

10

Mathematical models of contact problems with Coulomb friction in elastostatics have at least one solution provided that the coefficient of friction F is sufficiently small. On the other hand there is no information on the structure of solutions (uniqueness, multiplicity) in a general case. The situation is somewhat different

160

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161

0.8

1

2 1.8

0.7

2

1.6

4

0.6

3 1.4

0.5

1.2

branch 2

0.4

1

5

0.8

0.3

0.6

0.2

0.4

0.1

0.2 0

4

0 0

21

1

2

3

4

5

6

7

−3

−2

3 −1

0

1

2

3

2

3

3

1

1

−0.1

2.5 0.8

2 2 0.6

5

branch 2 1.5

4

0.4

3

1 0.2 0.5

2 3, 4

0

1

0 0

1

2

3

4

5

6

7

Fig. 15. Performance of the continuation routine in Example 6.6: branch 2.

for discrete versions of these problems. It is well-known that appropriate discretizations have a solution for any F belonging to a suitable class of coefficients and this solution is unique for F small enough. The aim of this paper was twofold: a) to show how the uniqueness property depends on the number of degrees of freedom, i.e. on the norm of finite element meshes used to get the discrete model; b) how (eventually) multiple solutions to the discrete problems can be captured. Sections 2 and 3 of this paper were devoted to a). We considered the orthotropic Coulomb friction law with friction coefficients which may depend on the solution itself. The discrete solutions are represented by fixed points of the mapping WhH acting on the contact zone. This mapping was defined by means of the mixed finite element formulation of auxiliary contact problems with Tresca friction. We presented conditions under which WhH is contractive, implying not only the uniqueness of the fixed point but also convergence of the method of successive approximations. Unfortunately, these conditions are mesh-dependent. To fulfill them, the coefficient F has to decay in an appropriate rate depending on the mesh norms. The rest of the paper deals with b). To understand better the structure of discrete solutions, we analyzed conditions ensuring the existence of local Lipschitz continuous branches of solutions assumed to be a function of F . This was done in Section 4 by using a variant of the implicit-function theorem for generalized equations. In a standard way of solving contact problems with Coulomb friction one finds some solution (which depends for example on the choice of the initial

−3

−2

−1

0

1

Fig. 16. Example 6.7. Between transition points 1 and 2 at the top, a goes forth, back and forth while km, 1 remains the same.

approximation) but without any further information on its position within other potential solutions. This is why we proposed a piecewise smooth variant of the Moore-Penrose continuation in Section 5 which allows us to follow branches of solutions parametrized by the friction coefficient F , the load vector f, etc. Unlike the classical Moore-Penrose continuation for smooth (differentiable) problems, some modifications in the prediction step ensuring the transition through points of non-differentiability have to be done. Using this technique, various graphs of solutions were obtained in the model examples presented in Section 6. Therefore, the proposed method seems to be promising for solving larger models arising from finite element discretizations of contact problems. Acknowledgements This research was supported under the grant No. 201/07/0294 of the Grant Agency of the Czech Republic and MSM 0021620839. T. Ligursky´ thanks also the Necˇas Center for Mathematical Modeling for its support. References [1] J. Necˇas, J. Jarušek, J. Haslinger, On the solution of the variational inequality to the Signorini problem with small friction, Boll. Unione Mat. Ital. V. Ser., B 17 (1980) 796–811.

J. Haslinger et al. / Comput. Methods Appl. Mech. Engrg. 205–208 (2012) 149–161 [2] C. Eck, J. Jarušek, M. Krbec, Unilateral contact problems. Variational methods and existence theorems, Vol. 270, Pure and Applied Mathematics (Boca Raton), Chapman & Hall/CRC, Boca Raton, FL, 2005. [3] P. Hild, Non-unique slipping in the Coulomb friction model in twodimensional linear elasticity, Q.J. Mech. Appl. Math. 57 (2) (2004) 225–235. [4] Y. Renard, A uniqueness criterion for the Signorini problem with Coulomb friction, SIAM J. Math. Anal. 38 (2) (2006) 452–467. [5] V. Janovsky´, Catastrophic features of Coulomb friction model, in: The mathematics of finite elements and applications IV, MAFELAP 1981, Proc. Conf., Uxbridge/Middlesex 1981 (1982), pp. 259–264. [6] T. Ligursky´, Theoretical analysis of discrete contact problems with Coulomb friction, submitted to Appl. Math. (2009). [7] S. Basseville, A. Leger, E. Pratt, Investigation of the equilibrium states and their stability for a simple model with unilateral contact and Coulomb friction, Arch. Appl. Mech. 73 (5–6) (2003) 409–420. [8] M. Raous, Quasistatic Signorini problem with Coulomb friction and coupling to adhesion, in: Wriggers, Peter, et al. (Eds.), New developments in contact problems, CISM Courses Lect, Vol. 384, Springer, Wien, 1999, pp. 101–178. [9] J. Haslinger, R. Kucˇera, O. Vlach, Bifurcations in contact problems with local Coulomb friction, in: K. Kunisch et al. (Eds.), Numerical Mathematics and Advanced Applications. Proceedings of ENUMATH 2007, Springer, Berlin, 2008, pp. 811–818. [10] P. Hild, Y. Renard, Local uniqueness and continuation of solutions for the discrete Coulomb friction problem in elastostatics, Quart. Appl. Math. 63 (2005) 553–573. [11] J. Haslinger, Approximation of the Signorini problem with friction, obeying the Coulomb law, Math. Methods Appl. Sci. 5 (1983) 422–437.

161

[12] J. Haslinger, I. Hlavácˇek, J. Necˇas, Numerical Methods for Unilateral Problems in Solid Mechanics, in: P.G. Ciarlet et al. (Eds.), Handbook of Numerical Analysis, Vol. IV, North-Holland, Amsterdam, 1995, pp. 313–485. [13] P. Clément, Approximation by finite element functions using local regularization, Rev. Franc. Automat. Inform. Rech. Operat. 9 (R-2) (1975) 77– 84. [14] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, Studies in Mathematics and its Applications, Vol. 4, North-Holland Publishing Company, Amsterdam - New York - Oxford, 1978. [15] J. Haslinger, R. Kucˇera, T. Ligursky´, Discretization and numerical realization of 3D elastostatic contact problems with orthotropic Coulomb friction and solution-dependent coefficients of friction, J. Comput. Appl. Math., (2010), submitted for publication. [16] S.M. Robinson, Strongly regular generalized equations, Math. Oper. Res. 5 (1980) 43–62. [17] S. Scholtes, Introduction to piecewise differentiable equations, Preprint No. 53/ 1994, Institut für Statistik und Mathematische Wirtschaftstheorie, Universität Karlsruhe, 1994. [18] A. Dhooge, W. Govaerts, Y.A. Kuznetsov, MATCONT: A MATLAB package for numerical bifurcation analysis of ODEs, ACM Trans. Math. Softw. 29 (2) (2003) 141–164. [19] F. Facchinei, J.-S. Pang, Finite-dimensional variational inequalities and complementarity problems, Springer Series in Operations Research, Vol. II, Springer, New York, NY, 2003. [20] V. Janovsky´, T. Ligursky´, Computing non unique solutions of the Coulomb friction problem, submitted to Math. Comput. Simul. (2009).