- Email: [email protected]
A FEM–BEM coupling method for a nonlinear transmission problem modelling Coulomb friction contact
Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466 www.elsevier.com/locate/cma
A FEM–BEM coupling method for a nonlinear transmission problem modelling Coulomb friction contact Matthias Maischak *, Ernst P. Stephan Institut fu¨r Angewandte Mathematik, Universita¨t Hannover, Welfengarten 1, D-30167 Hannover, Germany Received 12 December 2003; received in revised form 22 March 2004; accepted 22 March 2004
Abstract A nonlinear transmission problem modelling elastoplastic interface problems with contact and Coulomb friction is reduced to a boundary/domain variational inequality. We present a corresponding FEM/BEM coupling procedure and derive for its Galerkin solution an a priori error estimate. Furthermore we reformulate the problem as an equivalent saddle point problem whose discretization can be solved by the Uzawa algorithm. Convergence of the FEM/BEM coupling method is proved and numerical results are given. Ó 2004 Elsevier B.V. All rights reserved. Keywords: Transmission problem; Boundary elements; Finite elements; Steklov–Poincare´ operator; Friction
1. Introduction In this paper, we consider a nonlinear transmission problem which models unilateral contact with Coulomb friction. The problem under consideration can be seen as a simplified model of the displacement field of a nonlinear elastic body occupying a bounded domain X with smooth boundary C which is partly in unilateral frictional contact with a linear elastic body. We introduce a boundary/domain formulation (P) with a variational inequality. Using boundary integral operators the original transmission problem is reformulated as a coupling problem consisting of a weak formulation of a nonlinear differential operator in the bounded domain X and the Steklov–Poincare´ operator on the boundary of X. The use of boundary integral operators in formulations with variational inequalities goes back to the work of [3,10–12]. We present a FEM/BEM coupling method generalizing the approach in [9] which deals with a linear problem and
*
Corresponding author. E-mail addresses: [email protected] (M. Maischak), [email protected] (E.P. Stephan).
0045-7825/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.03.018
454
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
therefore uses a boundary integral operator formulation only. For another pure boundary element formulation, which uses regularization and therefore yields variational equations, see [6]. With a discrete version of the Steklov–Poincare´ operator we derive an a priori error estimate for the FEM/BEM coupling solution. Furthermore we present an equivalent saddle point formulation whose discretization can be solved with an Uzawa algorithm. Numerical experiments underline the usefulness and convergence of our approach. 2. Boundary/domain variational inequality Let X Rd be a bounded domain with Lipschitz boundary C. Let C ¼ Ct [ Cs where Ct and Cs are nonempty, disjoint and open in C. In the interior part X we consider a nonlinear partial differential equation, whereas in the exterior part Xc :¼ Rd n X we consider the Laplace equation. Then the friction contact problem under consideration reads: Find u such that divð.ðjrujÞ ruÞ ¼ f in X;
ð1Þ
Du ¼ 0 in Xc ;
ð2Þ
with the radiation condition uðxÞ ¼ a þ oð1Þ if d ¼ 2; uðxÞ ¼ Oðjxj2d Þ if dP3;
ðjxj ! 1Þ;
ð3Þ
where a 2 R is constant for any u, but varying with u. Here .: [0, 1) ! [0, 1) is a C1[0, 1) function with t Æ .(t) being monotonously increasing with t, .(t) 6 .0, (t Æ .(t)) 0 6 .1 and let .(t) + t Æ min{0, .(t)} P a > 0. 1 2 Writing u1 :¼ ujX and u2 :¼ ujXc , the tractions on C are given by .ðjru1 jÞ ou and ou with normal vector on on n pointing into Xc. We consider transmission conditions on Ct u1 jCt u2 jCt ¼ u0 jCt
and
.ðjru1 jÞ
ou1 ou2 j j ¼ t0 jCt on Ct on Ct
ð4Þ
and Coulomb friction conditions with given friction force g > 0 on Cs, ou1 ou2 ou1 6 g on Cs ; .ðjru1 jÞ j j ¼ t0 jCs ; .ðjru1 jÞ on Cs on Cs on .ðjru1 jÞ
ou1 ðu0 þ u2 u1 Þ þ gju0 þ u2 u1 j ¼ 0 on Cs : on
ð5Þ
ð6Þ
1 We remark that the above boundary conditions imply: if j.ðjru1 jÞ ou j < g then u0 + u2 u1 = 0, if on ou1 ou1 .ðjru1 jÞ on ¼ g then u0 + u2 u1 P 0, and if .ðjru1 jÞ on ¼ g then u0 + u2u1 6 0. These conditions reflect 1 the Coulomb law of friction. Here .ðjru1 jÞ ou corresponds to the tangential component of the traction in on 1 linear elasticity. The first implication means that the body sticks if the absolute value of .ðjru1 jÞ ou is less on than the friction coefficient. The remaining implications mean that if this bound is attained, the body slides. In order to derive a variational formulation of (1)–(6) we assume that the given data satisfy
f 2 L2 ðXÞ;
u0 2 H 1=2 ðCÞ;
together with Z f dx þ ht0 ; 1i ¼ 0;
t0 2 H 1=2 ðCÞ;
if d ¼ 2;
X
and we look for u 2 H 1 ðXÞ [ H 1loc ðXc Þ.
g 2 L1 ðCs Þ;
gP0;
ð7Þ
ð8Þ
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
455
To this end we introduce the following operators and forms. Let the function q be given by . through Z t s .ðsÞ ds: q : ½0; 1Þ ! ½0; 1Þ; t 7! qðtÞ ¼ 0
From the assumptions on .(t) it follows that 06qðtÞ6 12 .0 t2 and hence Z GðuÞ :¼ qðjrujÞ dx X
is finite for any u 2 H1(X) and strictly convex. The Fre´chet derivative of G is Z T DGðu; rÞ :¼ .ðjrujÞðruÞ rr dx 8u; r 2 H 1 ðXÞ
ð9Þ
X
and DG is uniformly monotone with respect to the semi-norm j jH 1 ðXÞ ¼ kr kL2 ðXÞ in H1(X), i.e., there exists a constant c > 0 such that cju rj2H 1 ðXÞ 6 DGðu; u rÞ DGðr; u rÞ
8u; r 2 H 1 ðXÞ
(see [3, Proposition 2.1]). The Coulomb friction is described by the following functional Z jðvÞ :¼ g jvj ds:
ð10Þ
ð11Þ
Cs
Next we define the boundary integral operators Z Z o V /ðzÞ :¼ 2 /ðxÞ Gðz; xÞ dsx ; K/ðzÞ :¼ 2 /ðxÞ
Gðz; xÞ dsx ; on x C C Z o o K 0 /ðzÞ :¼ 2 /ðxÞ
Gðz; xÞ dsx ; W /ðzÞ :¼ K/ðzÞ; onz onz C with the fundamental solution 8 1 > > log jx yj if d ¼ 2; > < 2p Gðx; yÞ ¼ d 2 > C > > 2 :
jx yj2d if dP3: 4pd=2 The Steklov–Poincare´ operator for the exterior problem S ¼ 12ðW þ ðK 0 IÞV 1 ðK IÞÞ;
ð12Þ
maps H1/2(C) continuously into its dual space H1/2(C) and the bilinear form hSv, vi is positive definite for ~ 1=2 ðCs Þ :¼ all v 2 H1/2(C), where hÆ, Æi is understood as L2-duality on C. Additionally, we define H s g. fv 2 H 1=2 ðCÞ : supp v C ~ 1=2 ðCs Þ we set Finally, for all ðu; vÞ; ðr; sÞ 2 X :¼ H 1 ðXÞ H Aðu; v; r; sÞ :¼ DGðu; rÞ þ hSðujC þ vÞ; rjC þ si;
ð13Þ
and we introduce for all (u, v) 2 X the linear functional k Z kðu; vÞ :¼ f u dx þ ht0 þ Su0 ; ujC þ vi:
ð14Þ
X
456
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
Then the variational formulation (P) of (1)–(6) reads as boundary/domain variational inequality: Find ð^ u; ^vÞ 2 X such that Að^ u; ^v; u ^ u; v ^vÞ þ jðvÞ jð^vÞPkðu ^ u; v ^vÞ 8ðu; vÞ 2 X :
ð15Þ
In order to describe the connection between (P) and the original friction problem we define the auxiliary problem (U). We consider the functional Z Z Z 1 2 Uðu1 ; u2 Þ :¼ Gðu1 Þ þ jru2 j dx f u1 t0 u2 jC ds: ð16Þ 2 Xc X C Then problem (U) consists in finding the minimizer ð^u1 ; ^u2 Þ in C of Uðu1 ; u2 Þ þ jððu2 u1 þ u0 ÞjCs Þ, where the convex set C is given by C :¼ fðu1 ; u2 Þ 2 H 1 ðXÞ H 1loc ðXc Þju1 jCt u2 jCt ¼ u0 jCt and u2 2 L2 g;
ð17Þ
L2 :¼ fv 2 H 1loc ðXc ÞjDv ¼ 0 in H 1 ðXc Þ and ðfor d ¼ 2 there exists a 2 R such thatÞ v satisfies (3)g: ð18Þ Theorem 1. Problem (U) is equivalent to (1)–(6) in the sense of distributions. Proof. C is convex and closed in H 1 ðXÞ H 1loc ðXc Þ (c.f. [3, Remark 4]). Due to standard arguments (see [3,4]) (U) is equivalent to: Find ð^ u1 ; ^ u2 Þ 2 C such that u2 Þ; ðu1 ^ u1 ; u2 ^ u2 ÞÞ þ jððu2 u1 þ u0 ÞjCs Þ jðð^u2 ^u1 þ u0 ÞjCs ÞP0 DUðð^ u1 ; ^ for all (u1, u2) 2 C. We have DUðð^ u1 ; ^ u2 Þ; ðu1 ; u2 ÞÞ ¼ DGð^ u1 ; u1 Þ þ
Z
r^ u2 ru2 dx Xc
Z
f u1 dx X
Z
t0 u2 jC ds:
ð19Þ
ð20Þ
C
Firstly, we show that (1)–(6) is a consequence of (19). Choose ðu1 ; u2 Þ :¼ ð^u1 þ gjX ; ^u2 þ gjXc Þ 2 C for some d g 2 C1 0 ðR Þ. Then, according to (19), (20), and integration by parts,
Z Z Z o^u1 o^u2 t0 g ds: 0 6 ðf þ divð.ðjr^ u1 jÞ r^ u1 ÞÞ g dx D^u2 g dx þ .ðjr^u1 jÞ on on X Xc C 1 Varying g 2 C 1 0 ðXÞ and g 2 C 0 ðXc Þ shows that (1) and (2) hold in the sense of distributions. Varying 1 d g 2 C 0 ðR Þ shows
.ðjru1 jÞ
ou1 ou2 ¼ þ t0 : on on
d Next, let g1 ; g2 2 C 1 u1 þ g1 jX ; u^2 þ g2 jXc Þ 2 C in an analogous way to obtain 0 ðR Þ and consider ðu1 ; u2 Þ :¼ ð^ Z Z o^ u1
ðg2 g1 Þ ds þ 0 6 .ðjr^ u1 jÞ gðjg2 g1 þ ^u2 ^u1 þ u0 j j^u2 ^u1 þ u0 jÞ ds: on C Cs
u2 ^ u1 þ u0 we can rewrite Defining w :¼ g2 g1 and ^v :¼ ^ Z Z o^ u1
w ds þ u1 jÞ gðjw þ ^vj j^vjÞ ds; 0 6 .ðjr^ on C Cs Z Z o^ u1
w ds þ u1 jÞ gjwj ds: 6 .ðjr^ on C Cs
ð21Þ ð22Þ
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
Consequently, due to g P 0, there holds Z Z o^u1 o^u1 .ðjr^ 6 w ds 6g on Cs : u jÞ gjwj ds; i:e: .ðjr^ u jÞ 1 1 on on Cs
457
ð23Þ
Cs
Choosing w ¼ ^v and w ¼ ^v in (21) we derive
Z o^ u1
^v þ gj^vj ds ¼ 0: .ðjr^ u1 jÞ on Cs Next, we show that on the other hand (19) is a consequence of (1)–(6). Due to (2) and (3) we have (u1, u2) 2 C. Multiplying (1), (2) with test functions (u, v) 2 C and integrating by parts gives Z Z ou1 ujC ds ¼ DGðu1 ; uÞ .ðjru1 jÞ fu dx; ð24Þ on C X Z Z ou2 vjC ds ¼ 0: ru2 rv dx þ ð25Þ Xc C on 1 2 ou ¼ t0 we obtain Combining (24) and (25) and .ðjru1 jÞ ou on on Z ou1 ðu v þ u2 u1 ÞjC ds: DUððu1 ; u2 Þ; ðu u1 ; v u2 ÞÞ ¼ .ðjru1 jÞ on C
ð26Þ
Therefore, using (6) and (5) we obtain DUððu1 ; u2 Þ; ðu u1 ; v u2 ÞÞ þ jððu v þ u0 ÞjCs Þ jððu2 u1 þ u0 ÞjCs Þ Z Z ou1 ððv u þ u0 Þ ðu2 u1 þ u0 ÞÞ ds þ gðjv u þ u0 j ju2 u1 þ u0 jÞ ds ¼ .ðjru1 jÞ on C C Z Z ou1 ¼ .ðjru1 jÞ ðv u þ u0 Þ ds þ gjv u þ u0 j dsP0: on C C Consequently, (u1, u2) 2 C is a solution of (19). h Next, we consider the functional 1 J ðu; vÞ :¼ GðuÞ þ hSðujC þ vÞ; ujC þ vi kðu; vÞ 2
8ðu; vÞ 2 X :
ð27Þ
Then problem (J) consists in finding the minimizer ð^u; ^vÞ in D :¼ {(u, v) 2 XjhS(ujC + v u0),1i = 0 if d = 2} of J(u, v) + j(v). The following result shows the equivalence of problems (U) and (J). Theorem 2. Assuming (7) there holds: (i) Let ð^ u1 ; ^ u2 Þ 2 C solve (U). Then (u,v) solves (J) where ðu; vÞ :¼ ð^u1 c; ðu0 þ ^u2 jC ^u1 jC ÞÞ for some c 2 R (c = 0 if d P 3). (ii) Let ð^ u; ^vÞ solve (J) and take a 2 R, arbitrary, for d = 2, whereas a = 0 if d P 3. Define u2 by the representation formula u2 :¼ 12 ðKð^ ujC þ ^v u0 þ aÞ V ðSð^ujC þ ^v u0 ÞÞÞ þ a and set u1 :¼ ^u þ a. Then (u1, u2) 2 C solves (U). Proof. Let (u, v) 2 D (not necessarily a solution of (J)) and define (u1,u2) as in (ii). Then we have (u1, u2) 2 C. For d = 2 we have hSu2jC,1i = ahS1, 1i.
458
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
Due to the properties of the Steklov–Poincare´ operator there holds ou2 ¼ Sðu2 jC aÞ; on
ð28Þ
where a is the constant in (3) for d = 2, and a = 0 for d P 3. For $u2 2 L2(Xc) there holds by Greens formula Z Z ou2 2 u2 jC ds: jru2 j dx ¼ Xc C on Hence, we have Z jru2 j2 dx ¼ hSðu2 jC aÞ; u2 jC i: Xc
Therefore, using (8) we obtain 1 c0 ¼ hSu0 ; u0 i þ ht0 ; u0 i: ð29Þ 2 Conversely, let (u1, u2) 2 C with U(u1, u2) < 1; hence $u2 2 L2(Xc). Moreover, we assume Du2 = 0. Thus, we may define (u, v) as in (i). Then we obtain (u, v) 2 D. By (29) and the shown correspondence of C and D the theorem follows. h Uðu1 ; u2 Þ ¼ c0 þ J ðu; vÞ;
Finally, we prove existence and uniqueness for the solution of our friction problem (P). Theorem 3. Assuming (7) there holds: Problem (P) is equivalent to problem (J) and has a unique solution. Proof. We have J(u, v) + j(v) ! +1 when k(u, v)kX ! 1, since j(Æ) is convex, proper and lower semi-continuous. On the other hand DG is uniformly monotone in H1(X) and the bilinear term with S is positive definite in H1/2(C). This implies that the solution of (J) exists due to [8, Chapter 1, Theorem 2.1] and is unique due to [8, Chapter 1, Theorem 2.2], since J is strictly convex. Analogously to [8, Chapter 1, Section 2] problem (J) is equivalent to problem ðP~ Þ, defined as: Find ð^ u; ^vÞ 2 D such that Að^ u; ^v; u ^ u; v ^vÞ þ jðvÞ jð^vÞPkðu ^ u; v ^vÞ
ð30Þ
for all (u, v) 2 D. For d P 3, problem (P) and problem ðP~ Þ are identical due to D = X. Next, we show the equivalence of (P) and ðP~ Þ in case of d = 2. Firstly, let ð^ u; ^vÞ 2 X be a solution of problem (P). Choosing ðu; vÞ ¼ ð1; ^vÞ in (15) and using (8) we have Z hSð^ ujC þ ^vÞ; 1iP f 1 dx þ ht0 þ Su0 ; 1i; i:e:; hSð^ujC þ ^v u0 Þ; 1i ¼ 0: X
Therefore, ð^ u; ^vÞ 2 D is the unique solution of problem ðP~ Þ. Next, let ð^ u; ^vÞ 2 D be the solution of problem ðP~ Þ. Let (u, v) 2 X. Due to hS1, 1i 5 0, there exists c 2 R, such that hS((u c)jC + v u0), 1i = 0. Then we have Að^ u; ^v; u ^ u; v ^vÞ þ jðvÞ jð^vÞ kðu ^ u; v ^vÞ ¼ Að^u; ^v; u c ^u; v ^vÞ þ jðvÞ jð^vÞ kðu c u^; v ^vÞ þ Að^u; ^v; c; 0Þ kðc; 0Þ PAð^u; ^v; c; 0Þ kðc; 0Þ ¼ hSð^ujC þ ^v u0 Þ; ci ¼ 0: Therefore, problem ðP~ Þ and problem (P) are equivalent and there exists a unique solution.
h
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
459
3. Discretization, a priori error estimate and FEM/BEM coupling ~ 1=2 H h1=2 , h 2 I be a family of finite-dimensional subspaces of H 1 ðXÞ H ~ 1=2 ðCs Þ H 1=2 ðCÞ, Let H 1h H h 1=2 1=2 ~ ,! H ~ ðCÞ, k h : H h1=2 ,! H 1=2 ðCÞ dewhere I (0, 1) with 0 2 I. For h 2 I let ih : H 1h ,! H 1 ðXÞ, jh : H h note the canonical imbeddings with their duals ih , jh , and k h . In order to approximate S we define 1
S h :¼ ih c W cih þ ih c ðI K 0 Þk h ðk h Vk h Þ k h ðI KÞcih :
ð31Þ
The discretized version (Ph) of problem (P) is given by ~ 1=2 H 1 ðXÞ H ~ 1=2 ðCs Þ such that Find ð^ uh ; ^vh Þ 2 X h :¼ H 1h H h uh ; ^vh ; uh ^ uh ; vh ^vh Þ þ jðvh Þ jð^vh ÞPkh ðuh ^uh ; vh ^vh Þ Ah ð^ for all (uh, vh) 2 Xh, where Ah :
H 1h
~ 1=2 H h
!
ðH 1h
~ 1=2 Þ H h
ð32Þ
is defined via
Ah ðuh ; vh ; rh ; sh Þ :¼ DGðuh ; vh Þ þ hS h ðuh jC þ vh Þ; rh jC þ sh i and kh : H 1h
ð33Þ
~ 1=2 H h
! R is defined via Z f uh dx þ ht0 ; uh jC þ vh i þ hih c ðW þ ðI K 0 Þk h ðk h Vk h Þ1 k h ðI KÞÞu0 ; uh jC þ vh i: kh ðuh ; vh Þ :¼ X
ð34Þ
Theorem 4. There exists a h0 > 0 such that for all h < h0 problem (Ph) has a unique solution (uh, vh) 2 Xh. Proof. We recall from [1] that there exists a constant c0 > 0 and h0 2 I such that for any h 2 I with h < h0 and any uh 2 H 1h there holds 2
hS h uh ; uh iPc0 kcih uh kH 1=2 ðCÞ :
ð35Þ
Therefore the arguments in the proof of Theorem 3 can be applied here.
h
There holds the following a priori error estimate. Theorem 5. For ð^ u; ^vÞ 2 X , the solution of problem (P), and ð^uh ; ^vh Þ 2 X h , the solution of problem (Ph), there holds with a constant C > 0, independent of h, for h < h0. 2
kð^ u^ uh ; ^v ^vh ÞkH 1 ðXÞH~ 1=2 ðC Þ s ( 6C
2
2
inf k^ u uh kH 1 ðXÞ þ inf ðk^v vh kH~ 1=2 ðC Þ þ k^v vh kL2 ðCÞ Þ
uh 2H h
~ 1=2 vh 2H h
þ distH 1=2 ðCÞ ðV 1 ðI KÞð^ u þ ^v
s
) 1=2 2 u0 Þ; H h Þ
1
:
ð36Þ
Proof. Let Eh :¼ S S h ¼ ðI K 0 Þ½V 1 k h ðk h Vk h Þ k h ðI KÞ. Then the symmetric mapping Eh: H1/2(C) ! H1/2(C) is continuous, positive semidefinite and bounded independent of h. It holds 1=2 kEh ðuÞkH 1=2 ðCÞ 6c0 distH 1=2 ðCÞ ðV 1 ðI KÞu; H h Þ 8u 2 H 1=2 ðCÞ, see [13, Lemma 15], extending the results for discretized Schur complements in [1]. Let (uh, vh) 2 Xh. Using (15) of problem (P) and (32) of problem (Ph) we obtain for h < h0 and any e > 0
460
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466 2
c kvh ^vh þ uh ^ uh kH 1=2 ðCÞ 6 hS h ðvh ^vh þ uh ^ uh Þ; vh ^vh þ uh ^ uh i þ Ah ð^uh ; ^vh ; uh ^uh ; vh ^vh Þ þ jðvh Þ jð^vh Þ kh ðuh ^ uh ; vh ^vh Þ þ Að^u; ^v; ^uh ^u; ^vh ^vÞ þ jð^vh Þ jð^vÞ kð^uh ^u; ^vh ^vÞ ^h Þ DGð^ ^h Þ þ Að^u; ^v; uh u^; vh ^vÞ kðuh u^; vh ^vÞ þ jðvh Þ jð^vÞ ¼ DGð^ u h ; uh u u; uh u ^h þ vh ^vh i hEh ðvh þ uh u0 Þ; uh u^h þ vh ^vh i hSð^ u þ ^v uh vh Þ; uh u ^h Þ DGð^ 6 DGð^ uh ; uh u u; uh ^ uh Þ þ Að^u; ^v; uh ^u; vh ^vÞ kðuh ^u; vh ^vÞ þ jðvh Þ jð^vÞ c0 ð1 þ kEh kÞe c00 ð1 þ kEh kÞ 2 2 kuh ^ k^u þ ^v uh vh kH 1=2 ðCÞ uh þ vh ^vh kH 1=2 ðCÞ þ 2 2e c000 þ kEh ð^ u þ ^v u0 k2H 1=2 ðCÞ : 2e Choosing a suitable e > 0 we obtain þ
2
38; ~c kvh ^vh þ uh ^ uh kH 1=2 ðCÞ 6 DGð^ uh ; uh ^uh Þ DGð^u; uh ^uh Þ þ Að^u; ^v; uh ^u; vh ^vÞ 2
kðuh ^u; vh ^vÞ þ jðvh Þ jð^vÞ þ kEh ð^u þ ^v u0 kH 1=2 ðCÞ 2
þ k^ u þ ^v uh vh kH 1=2 ðCÞ : Due to [3, Lemma 4.1] and the triangle inequality there exists a constant ^c > 0 such that there holds 2 2 2 ^c kð^ u^ uh ; ^v ^vh ÞkX 6 k^ u uh þ ^v vh kH 1=2 ðCÞ þ kuh ^uh þ vh ^vh kH 1=2 ðCÞ þ DGð^u; ^u ^uh Þ
DGð^ uh ; ^ u^ uh Þ:
ð37Þ
Therefore we obtain k^ u uh þ ^v vh k2H 1=2 ðCÞ 6 CfDGð^ uh ; uh ^ uÞ DGð^u; uh ^uÞ þ Að^u; ^v; uh ^u; vh ^vÞ kðuh ^ u; vh ^vÞ þ jðvh Þ jð^vÞ þ kEh ð^u þ ^v u0 k2H 1=2 ðCÞ þ k^ u þ ^v uh vh k2H 1=2 ðCÞ g:
ð38Þ
Now we have to bound the individual terms in (38). u uh ; ^vÞ in (15) and we obtain First, we set ðu; vÞ ¼ ðuh ; ^vÞ and ðu; vÞ ¼ ð2^ u; 0Þ: Að^ u; ^v; uh ^u; 0Þ ¼ kðuh ^
ð39Þ
Next, using (13), (14), (28), (4), (5) and the linearity of A in its second argument, we have u; vh ^vÞ ¼ Að^u; ^v; 0; vh ^vÞ kð0; vh ^vÞ Að^ u; ^v; uh ^u; vh ^vÞ kðuh ^ ¼ ht0 Sð^u þ ^v u0 Þ; vh ^vi
o^u ¼ .ðjr^ujÞ ; vh ^v 6 kgkL2 ðCs Þ kvh ^vkL2 ðCs Þ : on Additionally, there holds Z Z jðvh Þ jð^vÞ ¼ gðjvh j j^vjÞ ds6 g jvh ^vj ds6kgkL2 ðCs Þ kvh ^vkL2 ðCs Þ : Cs
Cs
In order to obtain an upper bound for the remaining term in (38) we proceed as follows. Firstly, we define the bilinear form Aw(Æ; Æ) by Z ~ Aw ðu; vÞ :¼ .ðjrwjÞru:rv dx; X
ð40Þ
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
461
where ~ .ðxÞ 2 Rdd is the Jacobian of x ! .(jxj)x, i.e., ~ .ðxÞ ¼ .ðjxjÞI dd þ .0 ðjxjÞ
x xt jxj
ðx 2 Rd Þ:
ð41Þ
The assumptions on . and the continuity of the differential operators imply that there exists a constant m > 0 such that Aw ðu; vÞ6mkukH 1 ðXÞ kvkH 1 ðXÞ
8u; v; w 2 H 1 ðXÞ:
ð42Þ
1
Next, for all u, v, w 2 H (X) we define with Rðw; uÞ :¼ .ðjrwjÞrw .ðjrujÞru ~ .ðruÞ rðw uÞ the form Qðw u; vÞ ¼
Z
½Rðw; uÞ :rv dx:
ð43Þ
X
Note that for all ^ u; ~ u, u 2 H1(X) we have u ^u; uÞ: DGð~ u; uÞ DGð^ u; uÞ ¼ Qð~ u^ u; uÞ þ A^u ð~ Since .(jxj)x is differentiable there holds u; wÞ6dðhÞk^ uh ^ ukH 1 ðXÞ kwkH 1 ðXÞ ; Qð^ uh ^
ð44Þ
with dðhÞ :¼
kRð^ uh ; ^ uÞkL2 ðXÞ krð^ uh ^ uÞkL2 ðXÞ
;
where d(h) ! 0 for k^ uh ^ ukH 1 ðXÞ ! 0. Now using (12), (40) and (44) we can estimate the remaining term in (38) DGð^ uh ; uh ^ uÞ DGð^ u; uh ^ uÞ ¼ Qð^ uh ^ u; uh ^uÞ þ A^u ð^uh ^u; uh ^uÞ 6 ðdðhÞ þ mÞk^uh u^kH 1 ðXÞ kuh u^kH 1 ðXÞ : Combining all estimates we obtain 2
kð^ u^ uh ; ^v ^vh ÞkX 6 Cf2kgkL2 ðCs Þ kvh ^vkL2 ðCs Þ þ ðdðhÞ þ mÞkð^uh ^u; ^vh ^vÞkX kðuh ^u; vh ^vÞkX 1=2 2
þ distH 1=2 ðCÞ ðV 1 ðI KÞð^u þ ^v u0 Þ; H h
Þ þ kð^u uh ; ^v vh Þk2X g:
Using standard arguments we obtain the desired a priori error estimate (36). h Next, we apply the finite element method (FEM) and the boundary element method (BEM) to compute with (Ph) approximations ð^ uh ; ^vh Þ 2 X h to the exact solution ð^u; ^vÞ 2 X of problem (P). Here, for simplicity we assume that X is a polyhedron and there are given regular triangulations xh of X into triangles or tetrahedrals T and regular partitions ch of Cs. Writing oCs for the boundary of Cs, we choose : uh j 2 P1 ðT Þ 8T 2 xh g; H h :¼ fuh 2 C 0 ðXÞ T 1=2
H h :¼ fvh 2 C 0 ðCÞ : vh jr 2 P1 ðrÞ 1=2 Hh
~ 1=2 H h
8r 2 ch g;
2
:¼ fwh 2 L ðCÞ : wh jr 2 P0 ðrÞ 8r 2 ch g;
:¼ fvh 2 C 0 ðCs Þ : vh jr 2 P1 ðrÞ 8r 2 ch ; vh ðoCs Þ ¼ 0g:
For this choice of subspaces Theorem 5 yields the convergence of the FEM/BEM coupling.
462
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
4. Saddle point formulation and Uzawa algorithm With the subset K of L2(Cs) given by K ¼ fr 2 L2 ðCs Þ : jrðxÞj61 a:e: on Cs g; we rewrite problem (P) as the following problem (S): ~ 1=2 ðCs Þ K such that ^Þ 2 H 1 ðXÞ H Find ð^ u; ^v; r Z ~ 1=2 ðCs Þ; ^ Að^ u; v; u; vÞ þ g^ r v ds ¼ kðu; vÞ 8ðu; vÞ 2 H 1 ðXÞ H
ð45Þ
Cs
^^v ¼ j^vj a:e: on Cs : r
ð46Þ
Then we have the following characterization of the solution of (P) which can be proved analogously to [8, Chapter 4, Theorem 2.2]. Theorem 6. The solution of problem (P) is solution of problem (S) and vice versa. ^. We define the active contact zone by Cc = {x 2 This result conforms the existence of a multiplier r Csju0(x) + u2(x) u1(x) = 0}, which can also be characterized by fx 2 Cs jj^ rðxÞj < 1g Cc , i.e. the body 1 sticks. The remaining part is called the nonactive contact zone. From (6) we have ð.ðjru1 jÞ ou þ on g^ rÞðu0 þ u2 u1 Þ ¼ 0 on Cs. In the active contact zone Cc we have u0 + u2 u1 = 0, whereas in the remain1 ^ 2 K. ing part Cs n Cc we have r ¼ g1 .ðjru1 jÞ ou . From (5) we conclude that indeed this gives r on Now, in view of Theorem 6 we can give another formulation of (P) as the following saddle point problem (L): ~ 1=2 ðCs Þ K such that ^Þ 2 H 1 ðXÞ H Find ð^ u; ^v; r ~ 1=2 ðCs Þ K; ^Þ6Lðu; v; r ^Þ 8ðu; v; rÞ 2 H 1 ðXÞ H Lð^ u; ^v; rÞ6Lð^ u; ^v; r ~ where the Lagrangian L : H 1 ðXÞ H
ð47Þ
1=2
ðCs Þ L2 ðCs Þ ! R is defined by Z 1 gr v ds kðu; vÞ: Lðu; v; rÞ :¼ GðuÞ þ hSðujC þ vÞ; ujC þ vi þ 2 Cs
ð48Þ
^Þ 2 X K be a solution of problem (S). Then ð^u; ^v; r ^Þ is the unique solution of problem Theorem 7. Let ð^ u; ^v; r (L) and vice versa. ~ 1=2 ðCs Þ Þ 2 X K be a solution of (L). Due to the second inequality in (47) with any v 2 H Proof. Let ð u; v; r and 2 R there holds Z 1 Þ6Lð ^Þ ¼ Lð Þ þ hSðujC þ vÞ; vi þ Lð u; v; r u; v þ v; r u; v; r g rv ds kð0; vÞ þ 2 hSv; vi: 2 Cs Therefore we have hSð ujC þ vÞ; vi þ
Z
~ 1=2 ðCs Þ: g rv ds ¼ kð0; vÞ 8v 2 H Cs
Let u 2 H1(X) and 0 6 t 6 1, then due to the second inequality in (47) there holds
ð49Þ
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
463
Þ6Lðu þ tðu uÞ; v; r ^Þ ¼ Lðu; v; r Þ þ Gð Lðu; v; r u þ tðu uÞÞ Gð uÞ 1 ujC i tkðu u; 0Þ þ t2 hSðujC ujC Þ; ujC ujC i: ¼ þthSðujC þ vÞ; ujC 2
By the convexity of G(Æ) we have 1 0 6 ð1 tÞGð uÞ þ tGðuÞ Gð uÞ þ thSð ujC þ vÞ; ujC ujC i tkðu u; 0Þ þ t2 hSðujC ujC Þ; ujC ujC i 2 1 2 ujC i kðu u; 0Þg þ t hSðujC ujC Þ; ujC ujC i: ¼ tfGðuÞ Gð uÞ þ hSð ujC þ vÞ; ujC 2 Therefore there holds for all u 2 H1(X) uÞ þ hSðujC þ vÞ; ujC i kðu; 0Þ: GðuÞ þ hSð ujC þ vÞ; ujC i kðu; 0Þ P Gð
ð50Þ
Taking the Fre´chet derivative we get DGð u; uÞ þ hSð ujC þ vÞ; ujC i ¼ kðu; 0Þ 8u 2 H 1 ðXÞ: Now, combining (49) and (51) we obtain Eq. (45) of problem (S). On the other hand, due to the first inequality in (47) there holds Z Þ ¼ Lð Lð u; v; rÞ6Lð u; v; r u; v; rÞ þ gð r rÞv ds 8r 2 K;
ð51Þ
ð52Þ
Cs
i.e., there holds Z gð r rÞv ds 8r 2 K: 06
ð53Þ
Cs
The latter inequality immediately implies the second equation (46) of (S). ^Þ 2 X K be a solution of (S). Conversely, let ð^ u; ^v; r Therefore, in view of Theorem 3 and (45) ð^ u; ^vÞ 2 X minimizes Z 1 J~ðu; vÞ :¼ GðuÞ þ hSðujC þ vÞ; ujC þ vi kðu; vÞ þ g^ rv; 2 Cs over X. Hence ^Þ6Lðu; v; r ^Þ Lð^ u; ^v; r
8ðu; vÞ 2 X :
Therefore, the second inequality in (47) holds. ^^v ¼ j^vj. Thus, for all r 2 K, i.e. jrj 6 1, there holds Furthermore, due to (46) we have r ^Þ^v ¼ r^v r ^^v ¼ r^v j^vj60. This gives ðr r Z ^Þ þ ^Þ^v ds6Lð^u; ^v; r ^Þ; Lð^ u; ^v; rÞ ¼ Lð^ u; ^v; r gðr r Cs
and the first inequality in (47) holds.
h
From Theorem 7 it follows that problem (P) can be solved by the Uzawa algorithm. In order to introduce this algorithm let PK be the projection of L2(Cs) onto K, i.e. there holds pointwise d # PK(d) = sup{1, inf(1, d)}. Algorithm 1 (Uzawa) (1) Choose r0 2 K. ~ 1=2 ðCs Þ such that (2) For n = 0, 1, 2, . . . determine ðun ; vn Þ 2 H 1 ðXÞ H
464
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
~ 8ðu; vÞ 2 H 1 ðXÞ H
Lðun ; vn ; rn Þ6Lðu; v; rn Þ
1=2
ðCs Þ;
ð54Þ
1=2
~ ðCs Þ such that i.e. find ðun ; vn Þ 2 H 1 ðXÞ H Z ~ 1=2 ðCs Þ: Aðun ; vn ; u; vÞ þ grn v ds ¼ kðu; vÞ 8ðu; vÞ 2 H 1 ðXÞ H
ð55Þ
Cs
(3) Set rnþ1 ¼ P K ðrn þ qgvn Þ;
ð56Þ
where q > 0 is a sufficiently small parameter that will be specified later on. (4) Repeat with 2. until a convergence criterion is satisfied. There holds the following convergence result. with c from (37). The Uzawa algorithm converges for arbitrary r0 2 K, i.e. ~ 1=2 ðCs Þ. u; ^vÞ strongly in X ¼ H 1 ðXÞ H ðun ; vn Þ ! ð^
Theorem 8. Let 0 < q < kgk22c
L1 ðCs Þ
Proof. Let ð^ u; ^vÞ be the solution of problem (P). Due to Theorem 6 and (46) we have ^ ¼ P K ð^ r r þ qg^vÞ:
ð57Þ
Since PK is a contraction, is follows from (56) and (57) that kr
nþ1
^k2L2 ðCs Þ 6 krn r ^ þ qgðvn ^vÞk2L2 ðCs Þ 6 krn r ^k2L2 ðCs Þ þ 2q r 2
Z
^Þðvn ^vÞ ds gðrn r Cs
2
þ q2 kgkL1 ðCs Þ kvn ^vkL2 ðCs Þ : u, v ¼ vn ^v in (45) and (55), we obtain from (37) Choosing u ¼ un ^ n
u; v c kðu ^
n
2 ^vÞkX
n
n
n
n
n
n
6 Aðu ; v ; u ^ u; v ^vÞ Að^u; ^v; u ^u; v ^vÞ ¼
Z
^Þðvn ^vÞ ds: gðrn r Cs
2
2
Using kvn ^vkL2 ðCs Þ 6 kvn ^vkH~ 1=2 ðC Þ we get s
2
2
2
2
^kL2 ðCs Þ krnþ1 r ^kL2 ðCs Þ Pqð2c qkgkL1 ðCs Þ Þkðun ^u; vn ^vÞkX : krn r 2 2c=kgkL1 ðCs Þ
2 qkgkL1 ðCs Þ Þ
we have ð2c > 0 and the sequence krn If 0 < q < therefore converging, because it has a lower bound. Hence, we get 2
ð58Þ ^k2L2 ðCs Þ r
is decreasing and
2
^kL2 ðCs Þ krnþ1 r ^kL2 ðCs Þ Þ ¼ 0 lim ðkrn r
n!1
2 ~ 1=2 ðCs Þ. and therefore limn!1 kðun ^ u; vn ^vÞkX ¼ 0, i.e. ðun ; vn Þ ! ð^u; ^vÞ strongly in X ¼ H ðXÞ H
h
5. Computations We present a numerical example for the friction contact problem (1)–(6) with . 1 and friction coefficient g 1. For X we take the L-shaped domain with vertices (0, 0), ð0; 14Þ, ð 14 ; 14Þ, ð 14 ; 14Þ, ð14 ; 14Þ and ð14 ; 0Þ. We take vanishing body forces, i.e. f = 0, Coulomb conditions on both large edges of X, and we choose u0 ¼ r2=3 sin 23 ðu p2Þ and t0 ¼ ono u0 ; therefore ht0, 1i = 0. As solver we apply the Uzawa algorithm (q = 25) with preconditioned CG-method for the inner problem (55). As preconditioner we use V-cycle
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
465
Table 1 Error in J, number of Uzawa iterations #it and computing times scpu(s) 1=2
jH 1h j
~ j jH h
jLhj
J h ð^uh ; ^vh Þ
dJ
ah
#it
scpu(s)
21 65 225 833 3201 12,545 49,665 197,633 788,481
7 15 31 63 127 255 511 1023 2047
7 15 31 63 127 255 511 1023 2047
0.821760029 0.826351913 0.828719162 0.829941080 0.830558446 0.830864354 0.831014083 0.831086830 0.831122044
0.0093951 0.0048032 0.0024359 0.0012140 0.0005967 0.0002907 0.0001410 0.0000683 0.0000331
0.938 0.980 1.005 1.025 1.037 1.044 1.046 1.045
3 3 3 3 3 3 3 3 3
0.01 0.02 0.04 0.15 0.62 2.72 12.14 52.40 279.12
multigrid with dampened Jacobi-smoother (x = 0.5) for both variables u and v. The Lagrange multiplier r ~ 1=2 . The error dJ, depending on the degree of freedom jÆj, is given by is approximated by rh 2 Lh :¼ H h J h ð^uh ; ^vh Þ J ð^ u; ^vÞ. Jh(uh, vh) is defined by replacing S in the definition (27) of J(u, v) by the approximated Steklov–Poincare´ operator (31). J ð^ u; ^vÞ ¼ 0:8311551 was computed by extrapolating the approximate values J h ð^ uh ; ^vh Þ. ah is the numerical convergence rate of J h ð^uh ; ^vh Þ with respect to h. #it is the number of Uzawa iterations and scpu(s) the cpu-time of the solver in seconds (Table 1).
6. Summary and further remarks The coupling of finite elements and boundary elements for nonlinear interface problems in elasticity has been investigated extensively (see e.g. [5,2]). In [11] Signorini contact problems in linear elasticity have been reduced to variational inequalities involving only boundary integral operators corresponding to the Lame´ equation. Therefore the above approach of reducing the friction problem to a boundary/domain variational inequality should work for elasticity problems as well. Also, the approach in [7] of using a mixed FEM/ BEM coupling method for Signorini problems can be extended to the above friction problem and the elasticity case as well. In applications such as metal forming processes usually plasticity occurs in a small bounded region, whereas large parts of the body are deformed elastically. Therefore also in case of a finite body a boundary/domain formulation seems appropriate, where the elastic part is modeled by boundary elements on the interface and only the plastic part is modeled by finite elements.
Acknowledgment This research was partially supported by the DFG grant no. STE 573/5-1.
References [1] C. Carstensen, Interface problem in holonomic elastoplasticity, Math. Methods Appl. Sci. 16 (1993) 819–835. [2] C. Carstensen, S.A. Funken, E.P. Stephan, On the adaptive coupling of FEM and BEM in 2-d elasticity, Numer. Math. 77 (1997) 187–221. [3] C. Carstensen, J. Gwinner, FEM and BEM coupling for a nonlinear transmission problem with Signorini contact, SIAM J. Numer. Anal. 34 (1997) 1845–1864. [4] J. Ce´a, Optimisation, the´orie et algorithmes, Paris, Dunod, 1971.
466
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
[5] M. Costabel, E.P. Stephan, Coupling of finite and boundary element methods for an elastoplastic interface problem, SIAM J. Numer. Anal. 27 (1990) 1212–1226. [6] C. Eck, O. Steinbach, W. Wendland, A symmetric boundary element method for contact problems with friction, Math. Comp. Simulat. 50 (1999) 43–61. [7] G. Gatica, M. Maischak, E.P. Stephan, Mixed FEM and BEM coupling for a transmission problem with Signorini contact, ifam preprint 58, 2002. [8] R. Glowinski, J.-L. Lions, R. Tre´molie`res, Numerical analysis of variational inequalitiesStudies in Mathematics and its Applications, 8, North-Holland Publishing, Amsterdam-New York, 1981. [9] H. Guediri, On a boundary variational inequality of the second kind modelling a friction problem, Math. Methods Appl. Sci. 25 (2002) 93–114. [10] J. Gwinner, E.P. Stephan, Boundary element convergence for a variational inequality of the second kind, Parametric Optimization and Related Topics, III (Gu¨strow, 1991), 1993, pp. 227–241. [11] J. Gwinner, E.P. Stephan, A boundary element procedure for contact problems in linear elastostatics, RAIRO Math. Mod. Numer. Anal. 27 (1993) 457–480. [12] H. Han, A direct boundary element method for Signorini problems, Math. Comp. 55 (1990) 115–128. [13] M. Maischak, E.P. Stephan, Adaptive hp-versions of BEM for Signorini problems, Appl. Numer. Math., in press.
A FEM–BEM coupling method for a nonlinear transmission problem modelling Coulomb friction contact Matthias Maischak *, Ernst P. Stephan Institut fu¨r Angewandte Mathematik, Universita¨t Hannover, Welfengarten 1, D-30167 Hannover, Germany Received 12 December 2003; received in revised form 22 March 2004; accepted 22 March 2004
Abstract A nonlinear transmission problem modelling elastoplastic interface problems with contact and Coulomb friction is reduced to a boundary/domain variational inequality. We present a corresponding FEM/BEM coupling procedure and derive for its Galerkin solution an a priori error estimate. Furthermore we reformulate the problem as an equivalent saddle point problem whose discretization can be solved by the Uzawa algorithm. Convergence of the FEM/BEM coupling method is proved and numerical results are given. Ó 2004 Elsevier B.V. All rights reserved. Keywords: Transmission problem; Boundary elements; Finite elements; Steklov–Poincare´ operator; Friction
1. Introduction In this paper, we consider a nonlinear transmission problem which models unilateral contact with Coulomb friction. The problem under consideration can be seen as a simplified model of the displacement field of a nonlinear elastic body occupying a bounded domain X with smooth boundary C which is partly in unilateral frictional contact with a linear elastic body. We introduce a boundary/domain formulation (P) with a variational inequality. Using boundary integral operators the original transmission problem is reformulated as a coupling problem consisting of a weak formulation of a nonlinear differential operator in the bounded domain X and the Steklov–Poincare´ operator on the boundary of X. The use of boundary integral operators in formulations with variational inequalities goes back to the work of [3,10–12]. We present a FEM/BEM coupling method generalizing the approach in [9] which deals with a linear problem and
*
Corresponding author. E-mail addresses: [email protected] (M. Maischak), [email protected] (E.P. Stephan).
0045-7825/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.03.018
454
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
therefore uses a boundary integral operator formulation only. For another pure boundary element formulation, which uses regularization and therefore yields variational equations, see [6]. With a discrete version of the Steklov–Poincare´ operator we derive an a priori error estimate for the FEM/BEM coupling solution. Furthermore we present an equivalent saddle point formulation whose discretization can be solved with an Uzawa algorithm. Numerical experiments underline the usefulness and convergence of our approach. 2. Boundary/domain variational inequality Let X Rd be a bounded domain with Lipschitz boundary C. Let C ¼ Ct [ Cs where Ct and Cs are nonempty, disjoint and open in C. In the interior part X we consider a nonlinear partial differential equation, whereas in the exterior part Xc :¼ Rd n X we consider the Laplace equation. Then the friction contact problem under consideration reads: Find u such that divð.ðjrujÞ ruÞ ¼ f in X;
ð1Þ
Du ¼ 0 in Xc ;
ð2Þ
with the radiation condition uðxÞ ¼ a þ oð1Þ if d ¼ 2; uðxÞ ¼ Oðjxj2d Þ if dP3;
ðjxj ! 1Þ;
ð3Þ
where a 2 R is constant for any u, but varying with u. Here .: [0, 1) ! [0, 1) is a C1[0, 1) function with t Æ .(t) being monotonously increasing with t, .(t) 6 .0, (t Æ .(t)) 0 6 .1 and let .(t) + t Æ min{0, .(t)} P a > 0. 1 2 Writing u1 :¼ ujX and u2 :¼ ujXc , the tractions on C are given by .ðjru1 jÞ ou and ou with normal vector on on n pointing into Xc. We consider transmission conditions on Ct u1 jCt u2 jCt ¼ u0 jCt
and
.ðjru1 jÞ
ou1 ou2 j j ¼ t0 jCt on Ct on Ct
ð4Þ
and Coulomb friction conditions with given friction force g > 0 on Cs, ou1 ou2 ou1 6 g on Cs ; .ðjru1 jÞ j j ¼ t0 jCs ; .ðjru1 jÞ on Cs on Cs on .ðjru1 jÞ
ou1 ðu0 þ u2 u1 Þ þ gju0 þ u2 u1 j ¼ 0 on Cs : on
ð5Þ
ð6Þ
1 We remark that the above boundary conditions imply: if j.ðjru1 jÞ ou j < g then u0 + u2 u1 = 0, if on ou1 ou1 .ðjru1 jÞ on ¼ g then u0 + u2 u1 P 0, and if .ðjru1 jÞ on ¼ g then u0 + u2u1 6 0. These conditions reflect 1 the Coulomb law of friction. Here .ðjru1 jÞ ou corresponds to the tangential component of the traction in on 1 linear elasticity. The first implication means that the body sticks if the absolute value of .ðjru1 jÞ ou is less on than the friction coefficient. The remaining implications mean that if this bound is attained, the body slides. In order to derive a variational formulation of (1)–(6) we assume that the given data satisfy
f 2 L2 ðXÞ;
u0 2 H 1=2 ðCÞ;
together with Z f dx þ ht0 ; 1i ¼ 0;
t0 2 H 1=2 ðCÞ;
if d ¼ 2;
X
and we look for u 2 H 1 ðXÞ [ H 1loc ðXc Þ.
g 2 L1 ðCs Þ;
gP0;
ð7Þ
ð8Þ
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
455
To this end we introduce the following operators and forms. Let the function q be given by . through Z t s .ðsÞ ds: q : ½0; 1Þ ! ½0; 1Þ; t 7! qðtÞ ¼ 0
From the assumptions on .(t) it follows that 06qðtÞ6 12 .0 t2 and hence Z GðuÞ :¼ qðjrujÞ dx X
is finite for any u 2 H1(X) and strictly convex. The Fre´chet derivative of G is Z T DGðu; rÞ :¼ .ðjrujÞðruÞ rr dx 8u; r 2 H 1 ðXÞ
ð9Þ
X
and DG is uniformly monotone with respect to the semi-norm j jH 1 ðXÞ ¼ kr kL2 ðXÞ in H1(X), i.e., there exists a constant c > 0 such that cju rj2H 1 ðXÞ 6 DGðu; u rÞ DGðr; u rÞ
8u; r 2 H 1 ðXÞ
(see [3, Proposition 2.1]). The Coulomb friction is described by the following functional Z jðvÞ :¼ g jvj ds:
ð10Þ
ð11Þ
Cs
Next we define the boundary integral operators Z Z o V /ðzÞ :¼ 2 /ðxÞ Gðz; xÞ dsx ; K/ðzÞ :¼ 2 /ðxÞ
Gðz; xÞ dsx ; on x C C Z o o K 0 /ðzÞ :¼ 2 /ðxÞ
Gðz; xÞ dsx ; W /ðzÞ :¼ K/ðzÞ; onz onz C with the fundamental solution 8 1 > > log jx yj if d ¼ 2; > < 2p Gðx; yÞ ¼ d 2 > C > > 2 :
jx yj2d if dP3: 4pd=2 The Steklov–Poincare´ operator for the exterior problem S ¼ 12ðW þ ðK 0 IÞV 1 ðK IÞÞ;
ð12Þ
maps H1/2(C) continuously into its dual space H1/2(C) and the bilinear form hSv, vi is positive definite for ~ 1=2 ðCs Þ :¼ all v 2 H1/2(C), where hÆ, Æi is understood as L2-duality on C. Additionally, we define H s g. fv 2 H 1=2 ðCÞ : supp v C ~ 1=2 ðCs Þ we set Finally, for all ðu; vÞ; ðr; sÞ 2 X :¼ H 1 ðXÞ H Aðu; v; r; sÞ :¼ DGðu; rÞ þ hSðujC þ vÞ; rjC þ si;
ð13Þ
and we introduce for all (u, v) 2 X the linear functional k Z kðu; vÞ :¼ f u dx þ ht0 þ Su0 ; ujC þ vi:
ð14Þ
X
456
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
Then the variational formulation (P) of (1)–(6) reads as boundary/domain variational inequality: Find ð^ u; ^vÞ 2 X such that Að^ u; ^v; u ^ u; v ^vÞ þ jðvÞ jð^vÞPkðu ^ u; v ^vÞ 8ðu; vÞ 2 X :
ð15Þ
In order to describe the connection between (P) and the original friction problem we define the auxiliary problem (U). We consider the functional Z Z Z 1 2 Uðu1 ; u2 Þ :¼ Gðu1 Þ þ jru2 j dx f u1 t0 u2 jC ds: ð16Þ 2 Xc X C Then problem (U) consists in finding the minimizer ð^u1 ; ^u2 Þ in C of Uðu1 ; u2 Þ þ jððu2 u1 þ u0 ÞjCs Þ, where the convex set C is given by C :¼ fðu1 ; u2 Þ 2 H 1 ðXÞ H 1loc ðXc Þju1 jCt u2 jCt ¼ u0 jCt and u2 2 L2 g;
ð17Þ
L2 :¼ fv 2 H 1loc ðXc ÞjDv ¼ 0 in H 1 ðXc Þ and ðfor d ¼ 2 there exists a 2 R such thatÞ v satisfies (3)g: ð18Þ Theorem 1. Problem (U) is equivalent to (1)–(6) in the sense of distributions. Proof. C is convex and closed in H 1 ðXÞ H 1loc ðXc Þ (c.f. [3, Remark 4]). Due to standard arguments (see [3,4]) (U) is equivalent to: Find ð^ u1 ; ^ u2 Þ 2 C such that u2 Þ; ðu1 ^ u1 ; u2 ^ u2 ÞÞ þ jððu2 u1 þ u0 ÞjCs Þ jðð^u2 ^u1 þ u0 ÞjCs ÞP0 DUðð^ u1 ; ^ for all (u1, u2) 2 C. We have DUðð^ u1 ; ^ u2 Þ; ðu1 ; u2 ÞÞ ¼ DGð^ u1 ; u1 Þ þ
Z
r^ u2 ru2 dx Xc
Z
f u1 dx X
Z
t0 u2 jC ds:
ð19Þ
ð20Þ
C
Firstly, we show that (1)–(6) is a consequence of (19). Choose ðu1 ; u2 Þ :¼ ð^u1 þ gjX ; ^u2 þ gjXc Þ 2 C for some d g 2 C1 0 ðR Þ. Then, according to (19), (20), and integration by parts,
Z Z Z o^u1 o^u2 t0 g ds: 0 6 ðf þ divð.ðjr^ u1 jÞ r^ u1 ÞÞ g dx D^u2 g dx þ .ðjr^u1 jÞ on on X Xc C 1 Varying g 2 C 1 0 ðXÞ and g 2 C 0 ðXc Þ shows that (1) and (2) hold in the sense of distributions. Varying 1 d g 2 C 0 ðR Þ shows
.ðjru1 jÞ
ou1 ou2 ¼ þ t0 : on on
d Next, let g1 ; g2 2 C 1 u1 þ g1 jX ; u^2 þ g2 jXc Þ 2 C in an analogous way to obtain 0 ðR Þ and consider ðu1 ; u2 Þ :¼ ð^ Z Z o^ u1
ðg2 g1 Þ ds þ 0 6 .ðjr^ u1 jÞ gðjg2 g1 þ ^u2 ^u1 þ u0 j j^u2 ^u1 þ u0 jÞ ds: on C Cs
u2 ^ u1 þ u0 we can rewrite Defining w :¼ g2 g1 and ^v :¼ ^ Z Z o^ u1
w ds þ u1 jÞ gðjw þ ^vj j^vjÞ ds; 0 6 .ðjr^ on C Cs Z Z o^ u1
w ds þ u1 jÞ gjwj ds: 6 .ðjr^ on C Cs
ð21Þ ð22Þ
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
Consequently, due to g P 0, there holds Z Z o^u1 o^u1 .ðjr^ 6 w ds 6g on Cs : u jÞ gjwj ds; i:e: .ðjr^ u jÞ 1 1 on on Cs
457
ð23Þ
Cs
Choosing w ¼ ^v and w ¼ ^v in (21) we derive
Z o^ u1
^v þ gj^vj ds ¼ 0: .ðjr^ u1 jÞ on Cs Next, we show that on the other hand (19) is a consequence of (1)–(6). Due to (2) and (3) we have (u1, u2) 2 C. Multiplying (1), (2) with test functions (u, v) 2 C and integrating by parts gives Z Z ou1 ujC ds ¼ DGðu1 ; uÞ .ðjru1 jÞ fu dx; ð24Þ on C X Z Z ou2 vjC ds ¼ 0: ru2 rv dx þ ð25Þ Xc C on 1 2 ou ¼ t0 we obtain Combining (24) and (25) and .ðjru1 jÞ ou on on Z ou1 ðu v þ u2 u1 ÞjC ds: DUððu1 ; u2 Þ; ðu u1 ; v u2 ÞÞ ¼ .ðjru1 jÞ on C
ð26Þ
Therefore, using (6) and (5) we obtain DUððu1 ; u2 Þ; ðu u1 ; v u2 ÞÞ þ jððu v þ u0 ÞjCs Þ jððu2 u1 þ u0 ÞjCs Þ Z Z ou1 ððv u þ u0 Þ ðu2 u1 þ u0 ÞÞ ds þ gðjv u þ u0 j ju2 u1 þ u0 jÞ ds ¼ .ðjru1 jÞ on C C Z Z ou1 ¼ .ðjru1 jÞ ðv u þ u0 Þ ds þ gjv u þ u0 j dsP0: on C C Consequently, (u1, u2) 2 C is a solution of (19). h Next, we consider the functional 1 J ðu; vÞ :¼ GðuÞ þ hSðujC þ vÞ; ujC þ vi kðu; vÞ 2
8ðu; vÞ 2 X :
ð27Þ
Then problem (J) consists in finding the minimizer ð^u; ^vÞ in D :¼ {(u, v) 2 XjhS(ujC + v u0),1i = 0 if d = 2} of J(u, v) + j(v). The following result shows the equivalence of problems (U) and (J). Theorem 2. Assuming (7) there holds: (i) Let ð^ u1 ; ^ u2 Þ 2 C solve (U). Then (u,v) solves (J) where ðu; vÞ :¼ ð^u1 c; ðu0 þ ^u2 jC ^u1 jC ÞÞ for some c 2 R (c = 0 if d P 3). (ii) Let ð^ u; ^vÞ solve (J) and take a 2 R, arbitrary, for d = 2, whereas a = 0 if d P 3. Define u2 by the representation formula u2 :¼ 12 ðKð^ ujC þ ^v u0 þ aÞ V ðSð^ujC þ ^v u0 ÞÞÞ þ a and set u1 :¼ ^u þ a. Then (u1, u2) 2 C solves (U). Proof. Let (u, v) 2 D (not necessarily a solution of (J)) and define (u1,u2) as in (ii). Then we have (u1, u2) 2 C. For d = 2 we have hSu2jC,1i = ahS1, 1i.
458
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
Due to the properties of the Steklov–Poincare´ operator there holds ou2 ¼ Sðu2 jC aÞ; on
ð28Þ
where a is the constant in (3) for d = 2, and a = 0 for d P 3. For $u2 2 L2(Xc) there holds by Greens formula Z Z ou2 2 u2 jC ds: jru2 j dx ¼ Xc C on Hence, we have Z jru2 j2 dx ¼ hSðu2 jC aÞ; u2 jC i: Xc
Therefore, using (8) we obtain 1 c0 ¼ hSu0 ; u0 i þ ht0 ; u0 i: ð29Þ 2 Conversely, let (u1, u2) 2 C with U(u1, u2) < 1; hence $u2 2 L2(Xc). Moreover, we assume Du2 = 0. Thus, we may define (u, v) as in (i). Then we obtain (u, v) 2 D. By (29) and the shown correspondence of C and D the theorem follows. h Uðu1 ; u2 Þ ¼ c0 þ J ðu; vÞ;
Finally, we prove existence and uniqueness for the solution of our friction problem (P). Theorem 3. Assuming (7) there holds: Problem (P) is equivalent to problem (J) and has a unique solution. Proof. We have J(u, v) + j(v) ! +1 when k(u, v)kX ! 1, since j(Æ) is convex, proper and lower semi-continuous. On the other hand DG is uniformly monotone in H1(X) and the bilinear term with S is positive definite in H1/2(C). This implies that the solution of (J) exists due to [8, Chapter 1, Theorem 2.1] and is unique due to [8, Chapter 1, Theorem 2.2], since J is strictly convex. Analogously to [8, Chapter 1, Section 2] problem (J) is equivalent to problem ðP~ Þ, defined as: Find ð^ u; ^vÞ 2 D such that Að^ u; ^v; u ^ u; v ^vÞ þ jðvÞ jð^vÞPkðu ^ u; v ^vÞ
ð30Þ
for all (u, v) 2 D. For d P 3, problem (P) and problem ðP~ Þ are identical due to D = X. Next, we show the equivalence of (P) and ðP~ Þ in case of d = 2. Firstly, let ð^ u; ^vÞ 2 X be a solution of problem (P). Choosing ðu; vÞ ¼ ð1; ^vÞ in (15) and using (8) we have Z hSð^ ujC þ ^vÞ; 1iP f 1 dx þ ht0 þ Su0 ; 1i; i:e:; hSð^ujC þ ^v u0 Þ; 1i ¼ 0: X
Therefore, ð^ u; ^vÞ 2 D is the unique solution of problem ðP~ Þ. Next, let ð^ u; ^vÞ 2 D be the solution of problem ðP~ Þ. Let (u, v) 2 X. Due to hS1, 1i 5 0, there exists c 2 R, such that hS((u c)jC + v u0), 1i = 0. Then we have Að^ u; ^v; u ^ u; v ^vÞ þ jðvÞ jð^vÞ kðu ^ u; v ^vÞ ¼ Að^u; ^v; u c ^u; v ^vÞ þ jðvÞ jð^vÞ kðu c u^; v ^vÞ þ Að^u; ^v; c; 0Þ kðc; 0Þ PAð^u; ^v; c; 0Þ kðc; 0Þ ¼ hSð^ujC þ ^v u0 Þ; ci ¼ 0: Therefore, problem ðP~ Þ and problem (P) are equivalent and there exists a unique solution.
h
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
459
3. Discretization, a priori error estimate and FEM/BEM coupling ~ 1=2 H h1=2 , h 2 I be a family of finite-dimensional subspaces of H 1 ðXÞ H ~ 1=2 ðCs Þ H 1=2 ðCÞ, Let H 1h H h 1=2 1=2 ~ ,! H ~ ðCÞ, k h : H h1=2 ,! H 1=2 ðCÞ dewhere I (0, 1) with 0 2 I. For h 2 I let ih : H 1h ,! H 1 ðXÞ, jh : H h note the canonical imbeddings with their duals ih , jh , and k h . In order to approximate S we define 1
S h :¼ ih c W cih þ ih c ðI K 0 Þk h ðk h Vk h Þ k h ðI KÞcih :
ð31Þ
The discretized version (Ph) of problem (P) is given by ~ 1=2 H 1 ðXÞ H ~ 1=2 ðCs Þ such that Find ð^ uh ; ^vh Þ 2 X h :¼ H 1h H h uh ; ^vh ; uh ^ uh ; vh ^vh Þ þ jðvh Þ jð^vh ÞPkh ðuh ^uh ; vh ^vh Þ Ah ð^ for all (uh, vh) 2 Xh, where Ah :
H 1h
~ 1=2 H h
!
ðH 1h
~ 1=2 Þ H h
ð32Þ
is defined via
Ah ðuh ; vh ; rh ; sh Þ :¼ DGðuh ; vh Þ þ hS h ðuh jC þ vh Þ; rh jC þ sh i and kh : H 1h
ð33Þ
~ 1=2 H h
! R is defined via Z f uh dx þ ht0 ; uh jC þ vh i þ hih c ðW þ ðI K 0 Þk h ðk h Vk h Þ1 k h ðI KÞÞu0 ; uh jC þ vh i: kh ðuh ; vh Þ :¼ X
ð34Þ
Theorem 4. There exists a h0 > 0 such that for all h < h0 problem (Ph) has a unique solution (uh, vh) 2 Xh. Proof. We recall from [1] that there exists a constant c0 > 0 and h0 2 I such that for any h 2 I with h < h0 and any uh 2 H 1h there holds 2
hS h uh ; uh iPc0 kcih uh kH 1=2 ðCÞ :
ð35Þ
Therefore the arguments in the proof of Theorem 3 can be applied here.
h
There holds the following a priori error estimate. Theorem 5. For ð^ u; ^vÞ 2 X , the solution of problem (P), and ð^uh ; ^vh Þ 2 X h , the solution of problem (Ph), there holds with a constant C > 0, independent of h, for h < h0. 2
kð^ u^ uh ; ^v ^vh ÞkH 1 ðXÞH~ 1=2 ðC Þ s ( 6C
2
2
inf k^ u uh kH 1 ðXÞ þ inf ðk^v vh kH~ 1=2 ðC Þ þ k^v vh kL2 ðCÞ Þ
uh 2H h
~ 1=2 vh 2H h
þ distH 1=2 ðCÞ ðV 1 ðI KÞð^ u þ ^v
s
) 1=2 2 u0 Þ; H h Þ
1
:
ð36Þ
Proof. Let Eh :¼ S S h ¼ ðI K 0 Þ½V 1 k h ðk h Vk h Þ k h ðI KÞ. Then the symmetric mapping Eh: H1/2(C) ! H1/2(C) is continuous, positive semidefinite and bounded independent of h. It holds 1=2 kEh ðuÞkH 1=2 ðCÞ 6c0 distH 1=2 ðCÞ ðV 1 ðI KÞu; H h Þ 8u 2 H 1=2 ðCÞ, see [13, Lemma 15], extending the results for discretized Schur complements in [1]. Let (uh, vh) 2 Xh. Using (15) of problem (P) and (32) of problem (Ph) we obtain for h < h0 and any e > 0
460
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466 2
c kvh ^vh þ uh ^ uh kH 1=2 ðCÞ 6 hS h ðvh ^vh þ uh ^ uh Þ; vh ^vh þ uh ^ uh i þ Ah ð^uh ; ^vh ; uh ^uh ; vh ^vh Þ þ jðvh Þ jð^vh Þ kh ðuh ^ uh ; vh ^vh Þ þ Að^u; ^v; ^uh ^u; ^vh ^vÞ þ jð^vh Þ jð^vÞ kð^uh ^u; ^vh ^vÞ ^h Þ DGð^ ^h Þ þ Að^u; ^v; uh u^; vh ^vÞ kðuh u^; vh ^vÞ þ jðvh Þ jð^vÞ ¼ DGð^ u h ; uh u u; uh u ^h þ vh ^vh i hEh ðvh þ uh u0 Þ; uh u^h þ vh ^vh i hSð^ u þ ^v uh vh Þ; uh u ^h Þ DGð^ 6 DGð^ uh ; uh u u; uh ^ uh Þ þ Að^u; ^v; uh ^u; vh ^vÞ kðuh ^u; vh ^vÞ þ jðvh Þ jð^vÞ c0 ð1 þ kEh kÞe c00 ð1 þ kEh kÞ 2 2 kuh ^ k^u þ ^v uh vh kH 1=2 ðCÞ uh þ vh ^vh kH 1=2 ðCÞ þ 2 2e c000 þ kEh ð^ u þ ^v u0 k2H 1=2 ðCÞ : 2e Choosing a suitable e > 0 we obtain þ
2
38; ~c kvh ^vh þ uh ^ uh kH 1=2 ðCÞ 6 DGð^ uh ; uh ^uh Þ DGð^u; uh ^uh Þ þ Að^u; ^v; uh ^u; vh ^vÞ 2
kðuh ^u; vh ^vÞ þ jðvh Þ jð^vÞ þ kEh ð^u þ ^v u0 kH 1=2 ðCÞ 2
þ k^ u þ ^v uh vh kH 1=2 ðCÞ : Due to [3, Lemma 4.1] and the triangle inequality there exists a constant ^c > 0 such that there holds 2 2 2 ^c kð^ u^ uh ; ^v ^vh ÞkX 6 k^ u uh þ ^v vh kH 1=2 ðCÞ þ kuh ^uh þ vh ^vh kH 1=2 ðCÞ þ DGð^u; ^u ^uh Þ
DGð^ uh ; ^ u^ uh Þ:
ð37Þ
Therefore we obtain k^ u uh þ ^v vh k2H 1=2 ðCÞ 6 CfDGð^ uh ; uh ^ uÞ DGð^u; uh ^uÞ þ Að^u; ^v; uh ^u; vh ^vÞ kðuh ^ u; vh ^vÞ þ jðvh Þ jð^vÞ þ kEh ð^u þ ^v u0 k2H 1=2 ðCÞ þ k^ u þ ^v uh vh k2H 1=2 ðCÞ g:
ð38Þ
Now we have to bound the individual terms in (38). u uh ; ^vÞ in (15) and we obtain First, we set ðu; vÞ ¼ ðuh ; ^vÞ and ðu; vÞ ¼ ð2^ u; 0Þ: Að^ u; ^v; uh ^u; 0Þ ¼ kðuh ^
ð39Þ
Next, using (13), (14), (28), (4), (5) and the linearity of A in its second argument, we have u; vh ^vÞ ¼ Að^u; ^v; 0; vh ^vÞ kð0; vh ^vÞ Að^ u; ^v; uh ^u; vh ^vÞ kðuh ^ ¼ ht0 Sð^u þ ^v u0 Þ; vh ^vi
o^u ¼ .ðjr^ujÞ ; vh ^v 6 kgkL2 ðCs Þ kvh ^vkL2 ðCs Þ : on Additionally, there holds Z Z jðvh Þ jð^vÞ ¼ gðjvh j j^vjÞ ds6 g jvh ^vj ds6kgkL2 ðCs Þ kvh ^vkL2 ðCs Þ : Cs
Cs
In order to obtain an upper bound for the remaining term in (38) we proceed as follows. Firstly, we define the bilinear form Aw(Æ; Æ) by Z ~ Aw ðu; vÞ :¼ .ðjrwjÞru:rv dx; X
ð40Þ
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
461
where ~ .ðxÞ 2 Rdd is the Jacobian of x ! .(jxj)x, i.e., ~ .ðxÞ ¼ .ðjxjÞI dd þ .0 ðjxjÞ
x xt jxj
ðx 2 Rd Þ:
ð41Þ
The assumptions on . and the continuity of the differential operators imply that there exists a constant m > 0 such that Aw ðu; vÞ6mkukH 1 ðXÞ kvkH 1 ðXÞ
8u; v; w 2 H 1 ðXÞ:
ð42Þ
1
Next, for all u, v, w 2 H (X) we define with Rðw; uÞ :¼ .ðjrwjÞrw .ðjrujÞru ~ .ðruÞ rðw uÞ the form Qðw u; vÞ ¼
Z
½Rðw; uÞ :rv dx:
ð43Þ
X
Note that for all ^ u; ~ u, u 2 H1(X) we have u ^u; uÞ: DGð~ u; uÞ DGð^ u; uÞ ¼ Qð~ u^ u; uÞ þ A^u ð~ Since .(jxj)x is differentiable there holds u; wÞ6dðhÞk^ uh ^ ukH 1 ðXÞ kwkH 1 ðXÞ ; Qð^ uh ^
ð44Þ
with dðhÞ :¼
kRð^ uh ; ^ uÞkL2 ðXÞ krð^ uh ^ uÞkL2 ðXÞ
;
where d(h) ! 0 for k^ uh ^ ukH 1 ðXÞ ! 0. Now using (12), (40) and (44) we can estimate the remaining term in (38) DGð^ uh ; uh ^ uÞ DGð^ u; uh ^ uÞ ¼ Qð^ uh ^ u; uh ^uÞ þ A^u ð^uh ^u; uh ^uÞ 6 ðdðhÞ þ mÞk^uh u^kH 1 ðXÞ kuh u^kH 1 ðXÞ : Combining all estimates we obtain 2
kð^ u^ uh ; ^v ^vh ÞkX 6 Cf2kgkL2 ðCs Þ kvh ^vkL2 ðCs Þ þ ðdðhÞ þ mÞkð^uh ^u; ^vh ^vÞkX kðuh ^u; vh ^vÞkX 1=2 2
þ distH 1=2 ðCÞ ðV 1 ðI KÞð^u þ ^v u0 Þ; H h
Þ þ kð^u uh ; ^v vh Þk2X g:
Using standard arguments we obtain the desired a priori error estimate (36). h Next, we apply the finite element method (FEM) and the boundary element method (BEM) to compute with (Ph) approximations ð^ uh ; ^vh Þ 2 X h to the exact solution ð^u; ^vÞ 2 X of problem (P). Here, for simplicity we assume that X is a polyhedron and there are given regular triangulations xh of X into triangles or tetrahedrals T and regular partitions ch of Cs. Writing oCs for the boundary of Cs, we choose : uh j 2 P1 ðT Þ 8T 2 xh g; H h :¼ fuh 2 C 0 ðXÞ T 1=2
H h :¼ fvh 2 C 0 ðCÞ : vh jr 2 P1 ðrÞ 1=2 Hh
~ 1=2 H h
8r 2 ch g;
2
:¼ fwh 2 L ðCÞ : wh jr 2 P0 ðrÞ 8r 2 ch g;
:¼ fvh 2 C 0 ðCs Þ : vh jr 2 P1 ðrÞ 8r 2 ch ; vh ðoCs Þ ¼ 0g:
For this choice of subspaces Theorem 5 yields the convergence of the FEM/BEM coupling.
462
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
4. Saddle point formulation and Uzawa algorithm With the subset K of L2(Cs) given by K ¼ fr 2 L2 ðCs Þ : jrðxÞj61 a:e: on Cs g; we rewrite problem (P) as the following problem (S): ~ 1=2 ðCs Þ K such that ^Þ 2 H 1 ðXÞ H Find ð^ u; ^v; r Z ~ 1=2 ðCs Þ; ^ Að^ u; v; u; vÞ þ g^ r v ds ¼ kðu; vÞ 8ðu; vÞ 2 H 1 ðXÞ H
ð45Þ
Cs
^^v ¼ j^vj a:e: on Cs : r
ð46Þ
Then we have the following characterization of the solution of (P) which can be proved analogously to [8, Chapter 4, Theorem 2.2]. Theorem 6. The solution of problem (P) is solution of problem (S) and vice versa. ^. We define the active contact zone by Cc = {x 2 This result conforms the existence of a multiplier r Csju0(x) + u2(x) u1(x) = 0}, which can also be characterized by fx 2 Cs jj^ rðxÞj < 1g Cc , i.e. the body 1 sticks. The remaining part is called the nonactive contact zone. From (6) we have ð.ðjru1 jÞ ou þ on g^ rÞðu0 þ u2 u1 Þ ¼ 0 on Cs. In the active contact zone Cc we have u0 + u2 u1 = 0, whereas in the remain1 ^ 2 K. ing part Cs n Cc we have r ¼ g1 .ðjru1 jÞ ou . From (5) we conclude that indeed this gives r on Now, in view of Theorem 6 we can give another formulation of (P) as the following saddle point problem (L): ~ 1=2 ðCs Þ K such that ^Þ 2 H 1 ðXÞ H Find ð^ u; ^v; r ~ 1=2 ðCs Þ K; ^Þ6Lðu; v; r ^Þ 8ðu; v; rÞ 2 H 1 ðXÞ H Lð^ u; ^v; rÞ6Lð^ u; ^v; r ~ where the Lagrangian L : H 1 ðXÞ H
ð47Þ
1=2
ðCs Þ L2 ðCs Þ ! R is defined by Z 1 gr v ds kðu; vÞ: Lðu; v; rÞ :¼ GðuÞ þ hSðujC þ vÞ; ujC þ vi þ 2 Cs
ð48Þ
^Þ 2 X K be a solution of problem (S). Then ð^u; ^v; r ^Þ is the unique solution of problem Theorem 7. Let ð^ u; ^v; r (L) and vice versa. ~ 1=2 ðCs Þ Þ 2 X K be a solution of (L). Due to the second inequality in (47) with any v 2 H Proof. Let ð u; v; r and 2 R there holds Z 1 Þ6Lð ^Þ ¼ Lð Þ þ hSðujC þ vÞ; vi þ Lð u; v; r u; v þ v; r u; v; r g rv ds kð0; vÞ þ 2 hSv; vi: 2 Cs Therefore we have hSð ujC þ vÞ; vi þ
Z
~ 1=2 ðCs Þ: g rv ds ¼ kð0; vÞ 8v 2 H Cs
Let u 2 H1(X) and 0 6 t 6 1, then due to the second inequality in (47) there holds
ð49Þ
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
463
Þ6Lðu þ tðu uÞ; v; r ^Þ ¼ Lðu; v; r Þ þ Gð Lðu; v; r u þ tðu uÞÞ Gð uÞ 1 ujC i tkðu u; 0Þ þ t2 hSðujC ujC Þ; ujC ujC i: ¼ þthSðujC þ vÞ; ujC 2
By the convexity of G(Æ) we have 1 0 6 ð1 tÞGð uÞ þ tGðuÞ Gð uÞ þ thSð ujC þ vÞ; ujC ujC i tkðu u; 0Þ þ t2 hSðujC ujC Þ; ujC ujC i 2 1 2 ujC i kðu u; 0Þg þ t hSðujC ujC Þ; ujC ujC i: ¼ tfGðuÞ Gð uÞ þ hSð ujC þ vÞ; ujC 2 Therefore there holds for all u 2 H1(X) uÞ þ hSðujC þ vÞ; ujC i kðu; 0Þ: GðuÞ þ hSð ujC þ vÞ; ujC i kðu; 0Þ P Gð
ð50Þ
Taking the Fre´chet derivative we get DGð u; uÞ þ hSð ujC þ vÞ; ujC i ¼ kðu; 0Þ 8u 2 H 1 ðXÞ: Now, combining (49) and (51) we obtain Eq. (45) of problem (S). On the other hand, due to the first inequality in (47) there holds Z Þ ¼ Lð Lð u; v; rÞ6Lð u; v; r u; v; rÞ þ gð r rÞv ds 8r 2 K;
ð51Þ
ð52Þ
Cs
i.e., there holds Z gð r rÞv ds 8r 2 K: 06
ð53Þ
Cs
The latter inequality immediately implies the second equation (46) of (S). ^Þ 2 X K be a solution of (S). Conversely, let ð^ u; ^v; r Therefore, in view of Theorem 3 and (45) ð^ u; ^vÞ 2 X minimizes Z 1 J~ðu; vÞ :¼ GðuÞ þ hSðujC þ vÞ; ujC þ vi kðu; vÞ þ g^ rv; 2 Cs over X. Hence ^Þ6Lðu; v; r ^Þ Lð^ u; ^v; r
8ðu; vÞ 2 X :
Therefore, the second inequality in (47) holds. ^^v ¼ j^vj. Thus, for all r 2 K, i.e. jrj 6 1, there holds Furthermore, due to (46) we have r ^Þ^v ¼ r^v r ^^v ¼ r^v j^vj60. This gives ðr r Z ^Þ þ ^Þ^v ds6Lð^u; ^v; r ^Þ; Lð^ u; ^v; rÞ ¼ Lð^ u; ^v; r gðr r Cs
and the first inequality in (47) holds.
h
From Theorem 7 it follows that problem (P) can be solved by the Uzawa algorithm. In order to introduce this algorithm let PK be the projection of L2(Cs) onto K, i.e. there holds pointwise d # PK(d) = sup{1, inf(1, d)}. Algorithm 1 (Uzawa) (1) Choose r0 2 K. ~ 1=2 ðCs Þ such that (2) For n = 0, 1, 2, . . . determine ðun ; vn Þ 2 H 1 ðXÞ H
464
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
~ 8ðu; vÞ 2 H 1 ðXÞ H
Lðun ; vn ; rn Þ6Lðu; v; rn Þ
1=2
ðCs Þ;
ð54Þ
1=2
~ ðCs Þ such that i.e. find ðun ; vn Þ 2 H 1 ðXÞ H Z ~ 1=2 ðCs Þ: Aðun ; vn ; u; vÞ þ grn v ds ¼ kðu; vÞ 8ðu; vÞ 2 H 1 ðXÞ H
ð55Þ
Cs
(3) Set rnþ1 ¼ P K ðrn þ qgvn Þ;
ð56Þ
where q > 0 is a sufficiently small parameter that will be specified later on. (4) Repeat with 2. until a convergence criterion is satisfied. There holds the following convergence result. with c from (37). The Uzawa algorithm converges for arbitrary r0 2 K, i.e. ~ 1=2 ðCs Þ. u; ^vÞ strongly in X ¼ H 1 ðXÞ H ðun ; vn Þ ! ð^
Theorem 8. Let 0 < q < kgk22c
L1 ðCs Þ
Proof. Let ð^ u; ^vÞ be the solution of problem (P). Due to Theorem 6 and (46) we have ^ ¼ P K ð^ r r þ qg^vÞ:
ð57Þ
Since PK is a contraction, is follows from (56) and (57) that kr
nþ1
^k2L2 ðCs Þ 6 krn r ^ þ qgðvn ^vÞk2L2 ðCs Þ 6 krn r ^k2L2 ðCs Þ þ 2q r 2
Z
^Þðvn ^vÞ ds gðrn r Cs
2
þ q2 kgkL1 ðCs Þ kvn ^vkL2 ðCs Þ : u, v ¼ vn ^v in (45) and (55), we obtain from (37) Choosing u ¼ un ^ n
u; v c kðu ^
n
2 ^vÞkX
n
n
n
n
n
n
6 Aðu ; v ; u ^ u; v ^vÞ Að^u; ^v; u ^u; v ^vÞ ¼
Z
^Þðvn ^vÞ ds: gðrn r Cs
2
2
Using kvn ^vkL2 ðCs Þ 6 kvn ^vkH~ 1=2 ðC Þ we get s
2
2
2
2
^kL2 ðCs Þ krnþ1 r ^kL2 ðCs Þ Pqð2c qkgkL1 ðCs Þ Þkðun ^u; vn ^vÞkX : krn r 2 2c=kgkL1 ðCs Þ
2 qkgkL1 ðCs Þ Þ
we have ð2c > 0 and the sequence krn If 0 < q < therefore converging, because it has a lower bound. Hence, we get 2
ð58Þ ^k2L2 ðCs Þ r
is decreasing and
2
^kL2 ðCs Þ krnþ1 r ^kL2 ðCs Þ Þ ¼ 0 lim ðkrn r
n!1
2 ~ 1=2 ðCs Þ. and therefore limn!1 kðun ^ u; vn ^vÞkX ¼ 0, i.e. ðun ; vn Þ ! ð^u; ^vÞ strongly in X ¼ H ðXÞ H
h
5. Computations We present a numerical example for the friction contact problem (1)–(6) with . 1 and friction coefficient g 1. For X we take the L-shaped domain with vertices (0, 0), ð0; 14Þ, ð 14 ; 14Þ, ð 14 ; 14Þ, ð14 ; 14Þ and ð14 ; 0Þ. We take vanishing body forces, i.e. f = 0, Coulomb conditions on both large edges of X, and we choose u0 ¼ r2=3 sin 23 ðu p2Þ and t0 ¼ ono u0 ; therefore ht0, 1i = 0. As solver we apply the Uzawa algorithm (q = 25) with preconditioned CG-method for the inner problem (55). As preconditioner we use V-cycle
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
465
Table 1 Error in J, number of Uzawa iterations #it and computing times scpu(s) 1=2
jH 1h j
~ j jH h
jLhj
J h ð^uh ; ^vh Þ
dJ
ah
#it
scpu(s)
21 65 225 833 3201 12,545 49,665 197,633 788,481
7 15 31 63 127 255 511 1023 2047
7 15 31 63 127 255 511 1023 2047
0.821760029 0.826351913 0.828719162 0.829941080 0.830558446 0.830864354 0.831014083 0.831086830 0.831122044
0.0093951 0.0048032 0.0024359 0.0012140 0.0005967 0.0002907 0.0001410 0.0000683 0.0000331
0.938 0.980 1.005 1.025 1.037 1.044 1.046 1.045
3 3 3 3 3 3 3 3 3
0.01 0.02 0.04 0.15 0.62 2.72 12.14 52.40 279.12
multigrid with dampened Jacobi-smoother (x = 0.5) for both variables u and v. The Lagrange multiplier r ~ 1=2 . The error dJ, depending on the degree of freedom jÆj, is given by is approximated by rh 2 Lh :¼ H h J h ð^uh ; ^vh Þ J ð^ u; ^vÞ. Jh(uh, vh) is defined by replacing S in the definition (27) of J(u, v) by the approximated Steklov–Poincare´ operator (31). J ð^ u; ^vÞ ¼ 0:8311551 was computed by extrapolating the approximate values J h ð^ uh ; ^vh Þ. ah is the numerical convergence rate of J h ð^uh ; ^vh Þ with respect to h. #it is the number of Uzawa iterations and scpu(s) the cpu-time of the solver in seconds (Table 1).
6. Summary and further remarks The coupling of finite elements and boundary elements for nonlinear interface problems in elasticity has been investigated extensively (see e.g. [5,2]). In [11] Signorini contact problems in linear elasticity have been reduced to variational inequalities involving only boundary integral operators corresponding to the Lame´ equation. Therefore the above approach of reducing the friction problem to a boundary/domain variational inequality should work for elasticity problems as well. Also, the approach in [7] of using a mixed FEM/ BEM coupling method for Signorini problems can be extended to the above friction problem and the elasticity case as well. In applications such as metal forming processes usually plasticity occurs in a small bounded region, whereas large parts of the body are deformed elastically. Therefore also in case of a finite body a boundary/domain formulation seems appropriate, where the elastic part is modeled by boundary elements on the interface and only the plastic part is modeled by finite elements.
Acknowledgment This research was partially supported by the DFG grant no. STE 573/5-1.
References [1] C. Carstensen, Interface problem in holonomic elastoplasticity, Math. Methods Appl. Sci. 16 (1993) 819–835. [2] C. Carstensen, S.A. Funken, E.P. Stephan, On the adaptive coupling of FEM and BEM in 2-d elasticity, Numer. Math. 77 (1997) 187–221. [3] C. Carstensen, J. Gwinner, FEM and BEM coupling for a nonlinear transmission problem with Signorini contact, SIAM J. Numer. Anal. 34 (1997) 1845–1864. [4] J. Ce´a, Optimisation, the´orie et algorithmes, Paris, Dunod, 1971.
466
M. Maischak, E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 194 (2005) 453–466
[5] M. Costabel, E.P. Stephan, Coupling of finite and boundary element methods for an elastoplastic interface problem, SIAM J. Numer. Anal. 27 (1990) 1212–1226. [6] C. Eck, O. Steinbach, W. Wendland, A symmetric boundary element method for contact problems with friction, Math. Comp. Simulat. 50 (1999) 43–61. [7] G. Gatica, M. Maischak, E.P. Stephan, Mixed FEM and BEM coupling for a transmission problem with Signorini contact, ifam preprint 58, 2002. [8] R. Glowinski, J.-L. Lions, R. Tre´molie`res, Numerical analysis of variational inequalitiesStudies in Mathematics and its Applications, 8, North-Holland Publishing, Amsterdam-New York, 1981. [9] H. Guediri, On a boundary variational inequality of the second kind modelling a friction problem, Math. Methods Appl. Sci. 25 (2002) 93–114. [10] J. Gwinner, E.P. Stephan, Boundary element convergence for a variational inequality of the second kind, Parametric Optimization and Related Topics, III (Gu¨strow, 1991), 1993, pp. 227–241. [11] J. Gwinner, E.P. Stephan, A boundary element procedure for contact problems in linear elastostatics, RAIRO Math. Mod. Numer. Anal. 27 (1993) 457–480. [12] H. Han, A direct boundary element method for Signorini problems, Math. Comp. 55 (1990) 115–128. [13] M. Maischak, E.P. Stephan, Adaptive hp-versions of BEM for Signorini problems, Appl. Numer. Math., in press.