Solution of the Contact Elasticity Problem Based on an Iterative Proximal Regularization Method for the Modified Lagrangian Functional Robert Namm, George Tsoy Computing Center of Far Eastern Branch Russian Academy of Sciences, Kim Yu Chen 65, 680000 Khabarovsk, Russia rnamm@yandex.ru,tsoy.dv@mail.ru http://www.ccfebras.ru/ Abstract. The method of successive approximation is considered for solving the contact elasticity problem which corresponds to the quasivariational Sig- norini problem. Auxiliary problems with given frictions arising from each exter- nal step of this method are solved by the Uzawa method with iterative prox- imal regularization of the modified Lagrangian functional. Stabilization of the sequence of auxiliary finite-element solutions of external steps of successive ap- proximation is investigated. Numerical results are considered. Keywords: proximal regularization, Lagrangian functional, saddle point, Uzawa method, Delaunay triangulation, finite element method 1 Introduction In this paper, we consider the contact problem of elasticity with friction between an elastic body Ω and absolutely rigid foundation (see [1,2]). Such problem corresponds to quasivariational Signorini inequality. The method of successive approximation is used for solving such inequality, on each external step of which the auxiliary contact problem with given friction arises. For solving semicoercive auxiliary problem we consider Uzawa method, based on modified Lagrangian functional and investigated in [3]. To overcome the problem of singularity (semicoercivity) we use the iterative proxi- mal regularization method of modified Lagrangian functional, considered in [4]. In this paper, we consider the finite element solution of quasivariational inequality. Domain Ω is decomposed by Delaunay triangulation with the help of mesh concentra- tion near the contact zone between elastic body and absolutely rigid foundation. 2 Formulation of the Problem Let Ω ⊂ R2 be a domain with a sufficiently regular boundary Γ , Γ = Γ0 Γ1 Γ2 , S S where Γ0 , Γ1 , Γ2 are open, mutually disjoint subsets of Γ ; moreover, Γ0 and Γ2 are Copyright c by the paper’s authors. Copying permitted for private and academic purposes. In: A. Kononov et al. (eds.): DOOR 2016, Vladivostok, Russia, published at http://ceur-ws.org Solution of the Contact Elasticity Problem 243 nonempty (that is, mesΓ0 > 0 and mesΓ2 > 0) (see fig. 1). We consider the con- tact problem of elasticity with friction between elastic body Ω and absolutely rigid foundation. Fig. 1: Contact between elastic body and absolutely rigid foundation For displacement vector v = (v1 , v2 ), we define the strain tensor   1 ∂vi ∂vj εij (v) = + , i, j = 1, 2, 2 ∂xj ∂xi and the stress tensor σij (v) = cijkm εkm (v). Here, i, j, k, m = 1, 2; cijkm = cjimk = ckmij , and repeated indices indicate summation. For given functions f = (f1 , f2 ), p = (p1 , p2 ) and F consider the following boundary value problem (see [1]): ∂σij − = fi in Ω, i = 1, 2 ∂xj (1) un = 0, στ = 0 on Γ0 σij nj = pi on Γ1 , i = 1, 2 On the contact surface Γ2 of the elastic body and the absolutely rigid support, we impose the conditions: un ≤ 0, σn ≤ 0, un σn = 0 on Γ2 (2) |στ | ≤ F |σn |, (F |σn | − |στ |)uτ = 0, uτ · στ ≤ 0 on Γ2 244 Robert Namm, George Tsoy Here, n = (n1 , n2 ) is the unit outward-pointing normal to Γ , un = u · n, uτ = u − un n; σi = σij nj , i = 1, 2; σ = (σ1 , σ2 ), σn = σij ni nj , στ = σ − σn n; and frictional coefficient F ≥ 0 on Γ2 . The main difficulty in the study and construction of numerical methods for solving this nonlinear boundary value problem is that the frictional force F |σn (u)| is a function of desired solution u. Define the sets V = v ∈ [W21 (Ω)]2 : vn ≡ v2 = 0 on Γ0 ,  K = {v ∈ V : vn ≤ 0 on Γ2 } . Assume that functions cijkm ∈ L∞ (Ω) (i, j, k, m = 1, 2), f ∈ [L2 (Ω)]2 , p ∈ [L2 (Γ1 )]2 and F ∈ L2 (Γ2 ). Suppose that the solution u to the boundary value problem (1), (2) exists and belongs to the space [W22 (Ω)]2 . Then u satisfies the quasivariational Signorini inequality for ∀ v ∈ K (see [1,2]) Z Z Z a(u, v − u) + F |σn (u)|(|vτ | − |uτ |) dΓ ≥ f · (v − u) dΩ + p · (v − u) dΓ , (3) Γ2 Ω Γ1 where Z Z a(u, v) = σij (u)εij (v) dΩ = cijkm εkm (u)εij (v) dΩ . (4) Ω Ω We use the successive approximation method for solving the quasivariational in- equality (see [1,2]): 1/2 1. Set starting friction force g 0 ∈ W2 (Γ2 ), g 0 ≥ 0 . 2. Find uk as a solution of auxiliary variational inequality Z Z Z a(u , v − u ) + g (|vτ | − |uτ |) dΓ ≥ f · (v − u ) dΩ + p · (v − uk ) dΓ . (5) k k k k k Γ2 Ω Γ1 3. Correct friction force g k+1 = F |σn (uk )|. The convergence of the successive approximation method to the solution to the quasivariational Signorini inequality (3) is still an open question. The existence of the solution is proved in the coercive case for a sufficiently small coefficient F only (see [1]). Variational inequality (5) is called the problem with given friction g k . It is equivalent to the following constrained non-differentiable minimization problem: J(v) = 1 a(v, v) + g k |vτ | dΓ − f · v dΩ − p · v dΓ → min ,  R R R 2 Γ2 Ω Γ1 (6) v∈K .  Under geometric form of domain Ω shown in fig. 1, functional J(v) is not strongly convex on all space V (see [1,2]). The kernel of bilinear form a(u, v) is not trivial and consists of vector-functions ρ = (a, 0), where a is an arbitrary number. Solution of the Contact Elasticity Problem 245 However, if Z Z f1 dΩ + p1 dΓ > 0 Ω Γ1 then J(v) → +∞ under kvk[W21 (Ω)]2 → ∞, v ∈ K, that is, J is a coercive functional on the set K and therefore problem (6) is solvable. 3 Application of Dual Schemes for Solving the Auxiliary Problem with Given Friction For problem (6), on set V × L2 (Γ2 ) we define the classical Lagrangian functional Z L(v, l) = J(v) + lvn dΓ . Γ2 Denote by (L2 (Γ2 ))+ the set of nonnegative square integrable functions on Γ2 . Definition 1. A pair (v ∗ , l∗ ) ∈ V ×(L2 (Γ2 ))+ is called a saddle point of the Lagrangian functional L(v, l) if it satisfies the two sided inequalities L(v ∗ , l) ≤ L(v ∗ , l∗ ) ≤ L(v, l∗ ) ∀ (v, l) ∈ V × (L2 (Γ2 ))+ . It was shown in [3] that, if a solution u of auxiliary problem (6) belongs to space [W22 (Ω)]2 and mes{x ∈ Γ2 : σn (u) < 0} > 0, then u is a unique solution to the problem (6) and the pair (u, −σn (u)) is a unique saddle point of L(v, l). Because of the second part of saddle point is equal to −σn (u), we can find accurately the frictional force on the next step of successive approximation method. However, application of the dual Uzawa method with the classical Lagrangian func- tional L(v, l) does not guarantee convergence to a saddle point in the semicoercive case (see [6,7]). To overcome this difficulty the modified Lagrangian functional M (v, l) was consid- ered on space V × L2 (ΓK ) (see [8],[9]): Z n 1 2 o (l + rvn )+ − l2 dΓ ,  M (v, l) = J(v) + 2r Γ2 where (l + rvn )+ ≡ max{0, l + rvn }, r > 0 – const. Definition 2. A pair (v, l) ∈ V × (L2 (Γ2 )) is called a saddle point of the modified Lagrangian functional M (v, l) if it satisfies the two sided inequalities M (v, l) ≤ M (v, l) ≤ M (v, l) ∀ (v, l) ∈ V × (L2 (Γ2 )) . 246 Robert Namm, George Tsoy Definition 2 is different from Definition 1. Here we use all space L2 (Γ2 ) instead of L2 (Γ2 )+ in the first definition. However, it is known that the functionals L(v, l) and M (v, l) have the same set of saddle points (see [3], [10]). For finding a saddle point of modified Lagrangian functional M (v, l) we use the method based on the combination of Uzawa method with the proximal regularization (see [4]). We note that similar method in finite-dimensional case was investigated in [5]. According to this method, the sequence {(um , lm )} is generated as follows. 1 1. Assign the initial approximation (u0 , l0 ) ∈ V × W22 (Γ2 ). 2. Find um+1 such that kum+1 − um+1 k[W21 (Ω)]2 ≤ δm , (7) where   m+1 1 u = arg min M (v, l ) + kv − um k2[L2 (Ω)]2 m , (8) v∈V 2 ∞ X δm > 0, δm < ∞, (9) m=1 3. Calculate the next value of the dual variable by the formula lm+1 = (lm + rum + n) . (10) Criterion (7) implies that the exact solution um+1 is replaced by its approximation m+1 u , obtained by solving problem (8) numerically using the finite element method. In this case parameter δm can be interpreted as the error of the numerical solution. The regularizing term 12 kv−um k2[L2 (Ω)]2 ensures that the functional to be minimized in (8) is strongly convex on V . This guarantees that the auxiliary problem (5) is uniquely solvable. Application of the modified Lagrangian fucntional allows as to find saddle point efficiently in comparison with duality methods based on classical Lagrangian functional. It should be noted that in article [11] it was built and investigated a wide class of iterative solution methods of the semicoercive variational inequalities including the iterative proximal regularization method. Discuss discretization of the variational problem (8). 4 Finite Element Discretization √ Domain Ω is taken in the form of a trapezoid (fig. 1) with sides 2, 1, 1.5, 5/2. Under the proposed algorithm (7)–(10) on each step of the iteration process, we regard the minimization problem of the strongly convex functional 1 Mm (v) = M (v, lm ) + kv − um k2[L2 (Ω)]2 → min, v∈V (11) 2 Solution of the Contact Elasticity Problem 247 For the solving of the problem (11) finite element method have been applied. Using the Delaunay triangulation, we triangulate Ω into a set of triangles Tk (fig. 2), so SN that Ω = 1 t Tk , where Nt is the number of triangles. Thus we have a finite element consisting of triangles, which have one point in common at the triangulation nodes. Condensation of the mesh occurs near the contact zone Γ2 . Fig. 2: Delaunay triangulation of Ω Enumerate triangulation nodes from the top down, from 1 to N . For each node i it is defined basis function ϕi (x, y), for which ϕi (xi , y i ) = 1 and for all neighboring nodes j: ϕi (xj , y j ) = 0. For basis functions ϕi we take piecewise linear functions (see [12]). Let us introduce the following notation: h - maximum edge length of the triangula- tion Tk , Ph = {D1 , · · · , DN } - triangulation node set , Ih = {M1 , · · · , MR } - boundary node set on Γ2 , Vh - linear span of the basis functions ϕi (x, y), uh = (uh1 , uh2 ) - piecewise interpolation of the exact solution u: N (i) (i) X uhi (x, y) = tj ϕj (x, y), for i = 1, 2 and tj ∈ R . j=1 Note that since Ω is a polygon, then embedding Vh ⊂ V is provided. Thus we can substitute the problem (8) by finite-element problem: um+1 = arg min {Mm (v)} . (12) v∈Vh There is shown in [13] that error estimate is true for the exact solution sequence {um }: kum − um k[W21 (Ω)]2 ≤ Ch1/2 , C > 0 − const, m = 1, 2, . . . , (13) and finite-element solution sequence {um } converges on V to solution of the auxiliary problem (6) under h → 0. 248 Robert Namm, George Tsoy Let us introduce the vector t = (t1 , t2 , . . . , t2N ), where the first N its components correspond to uh1 (Di (xi , y i )) and last N components correspond to uh2 (Di (xi , y i )). Then the minimization problem (11) reduces to finding optimal values ti . For this purpose, we use coordinate descent method. Let stiffness matrix A and vector F be defined as follows    Z  A = a(ψi , ψj ) + ψi ψj dΩ ,   Ω i,j=1,2N   Z Z  F = (f + um )ψi dΩ + pψi dΓ ,   Ω Γ1 i=1,2N where ψi are vector-functions such that ( (ϕi , 0) , i≤N , ψi = (0, ϕ(i−N ) ) , i>N . Thus we have the computational formula for all nodes Di ∈ Ph \ Ih   1  X X ts+1 i =− Aij ts+1 j + Aij tsj − Fi  . (14) Aii ji Consider next a problem for nodes Di ∈ Ih . For this purpose, we pass from variables u1 , u2 to un , uτ , where α is an acute angle of the domain Ω: ( u1 = un cos α + uτ sin α , u2 = −un sin α + uτ cos α . For evaluation the boundary integral on Γ2 Z n 1  m 2 o (l + rvn )+ − (lm )2 dΓ 2r Γ2 m we use the trapezoidal rule for vn > lr Z n Z  1  m + 2 m 2 o r  lm vn + vn2 dΓ ≈  (l + rvn ) − (l ) dΓ = 2r 2 Γ2 Γ2 R−1  |M , M | X r j j+1 lm (Mj )vn (Mj ) + lm (Mj+1 )vn (Mj+1 ) + (vn2 (Mj ) + vn2 (Mj+1 )) . j=1 2 2 Let |Mj , Mj+1 | = H for j = 1, R − 1, then introduce some notations An = Aii cos2 α + Ai+N i+N sin2 α − 2Aii+N cos α sin α , Solution of the Contact Elasticity Problem 249 Aτ = Aii cos α sin α + Aii+N cos2 α − Ai+N i sin2 α − Ai+N i+N cos α sin α ,    ∂  m Z lm (Di )H , Di ∈ {M2 , · · · , MR−1 } , L(Di ) = l vn dΓ ≈ m  H ∂tn l (Di ) , Di ∈ {M1 , MR } , Γ2 2    ∂2  r 2 Z rH , Di ∈ {M2 , · · · , MR−1 } , B(Di ) = 2 vn dΓ  ≈ H ∂tn 2 r , Di ∈ {M1 , MR } , Γ2 2 where tn = vn (Di ). To find the optimal tn compute  1 X X φi = − tsτ Aτ + Aij cos α ts+1 j + Aij cos α tsj − An ji j6=i+N X X  − Ai+N j sin α ts+1 j − Ai+N j sin α tsj − Fi cos α + Fi+N sin α . ji j6=i+N Let us denote an expression within the parentheses by ωi . Then m φi ≤ − l (Di)   φi , r ,  ts+1 n = (15) − ωi + L(Di ) , m φi > − l (Di)   r , An + B(Di ) A similar approach can be used for tτ on (i+N )th step of the Gauss-Seidel method. Approximate R−1  |Mj , Mj+1 | Z X g k |vt | dΓ ≈ g k (Mj )|vτ (Mj )| + g k (Mj+1 )|vτ (Mj+1 )| . j=1 2 Γ2 Let A∗n = Aii cos α sin α − Aii+N sin2 α + Ai+N i cos2 α − Ai+N i+N cos α sin α , A∗τ = Aii sin2 α + Ai+N i+N cos2 α + 2Aii+N cos α sin α ,  g(Di )H , Di ∈ {M2 , · · · , MR−1 } , G(Di ) = H g(Di ) , Di ∈ {M1 , MR } , 2 X X X s ∗ s+1 ξ i = t n An + Aij sin α tj + Aij sin α tsj + Ai+N j cos α ts+1 j + ji+N ji+N 250 Robert Namm, George Tsoy (a) un on Γ2 (b) |σn | on Γ2 Fig. 3: Computational result #1 Then the computational formula for tτ can be written as follows  − ξi + G(Di ) ,  ξi < −G(Di ) , A∗τ         tτ = − ξi − G(Di ) , s+1 ξi > G(Di ) , (16)     A∗τ      0 , −G(Di ) ≤ ξi ≤ G(Di ) . Iterations over (14) – (16) stop when min |ts+1 i − tsi | < 10−10 , i = 1, 2N . i Cessation condition of the iteration process for the Uzawa method and successive approximation method can be defined as follows min |lim+1 − lim | < 10−8 , min |gik+1 − gik | < 10−8 . i=1,R i=1,R 4.1 Numerical Solution of the Quasivariational Inequality We present the results of numerical computations for solving problem (3). We assume a volume load f = (f1 , f2 ) = (0, 0), boundary loading from the left p1 |Γ1 = 27 mPa and from the right p1 |Γ1 = −27 mPa sides, p2 |Γ1 = 0, frictional coefficient F = 0.3, Young’s modulus E = 73000 mPa, Poisson’s ratio µ = 0.34, constant r = 108 . Numerical solution of the considering problem is shown in fig. 3. The graphs of un and |σn | show that body Ω is detached from absolutely rigid foundation at the vertex of the obtuse angle. It follows from the fact that un < 0 and |σn | = 0 at this vertex. Let us introduce the example, when body Ω sticks together with absolutely rigid foundation. We change boundary loading from the right side p1 |Γ1 = −21.6 mPa. Solution of the Contact Elasticity Problem 251 Figure 4 shows that un takes on a value close to zero on Γ2 (less than calculation accuracy ε = 10−10 ), |σn | is everywhere positive, whence it follows that elastic body comes in contact with a rigid foundation at all points. The graph of |σn | states that body is compressed and the maximum stress is achieved at the vertex of the obtuse angle. (a) un on Γ2 (b) |σn | on Γ2 Fig. 4: Computational result #2 Numerical results confirmed that modified Lagrangian functionals effectively remove conditions like un ≤ 0 on Γ2 when passing to the unconstrained minimization problem. Besides, it was revealed that successive approximation method is more effective at larger r (r = 106 , 108 ). 252 Robert Namm, George Tsoy References 1. Hlaváček, I., Haslinger, J., Nečas, J., Lovı́šek, J.: Solution of Variational Inequalities in Mechanics. Springer, New York (1988) 2. Kikuchi, N., Oden, T.: Contact problem in elasticity: a study of variational inequalities and finite element methods. SIAM, Philadelphia (1988) 3. Vikhtenko, E.M., Namm, R.V.: Duality scheme for solving the semicoercive signorini problem with friction. Comput. Math. Math. Phys. 47(12), 2023–2036 (2007) 4. Vikhtenko, E.M., Namm, R.V.: Iterative proximal regularization of the modified La- grangian functional for solving the quasi-variational Signorini inequality. Comput. Math. Math. Phys. 48(9), 1–9 (2008) 5. Rockafellar, R.T.: Augmented lagrangians and applications of the proximal point algo- rithm in convex programming, Math. Oper. Res. 1, 97–116 (1976) 6. Glowinski, R., Lions, J.L., Tremolieres, R.: Numerical analysis of variational inequalities. North-Holland, Amsterdam (1981) 7. Glowinski, R.: Numerical methods for nonlinear variational problems. Springer-Verlag, New York (1984) 8. Woo, G., Namm, R.V., Sachkoff, S.A.: An iterative method based on a modified La- grangian functional for finding a saddle point in the semicoercive Signorini problem. Comput. Math. Math. Phys. 46(1), 26–36 (2006) 9. Woo, G., Kim, S., Namm, R.V., Sachkoff, S.A.: Iterative proximal regularization method for finding a saddle point in the semicoercive Signorini problem. Comput. Math. Math. Phys. 46(11), 2024–2031 (2006) 10. Vikhtenko, E.M., Woo, G., Namm, R.V.: Solution methods for semicoerive variational inequalities of mechanics based on modified Lagrangian functionals. Dal’nevost. Mat. Zh. 14(1), 6–17 (2014) 11. Konnov, I., Gwinner, J.: A strongly convergent combined relaxation in Hilbert Spaces. Numer. Func. Anal. Opt. 35, 1066–1077 (2014) 12. Marchuk, G.I., Agoshkov, V.I.: An Introduction to Projective Grid Methods [in Russian]. Fizmatlit, Moscow (1981) 13. Namm, R.V., Sachkoff, S.A.: Solving the quasi-variational Signorini inequality by the method of successive approximations. Comput. Math. Math. Phys. 49(5), 805–814 (2009)