Variational Method for Solving Contact Problem of Elasticity Robert Namma, Georgiy Tsoya, Ellina Vikhtenkob and Gyungsoo Wooc a Computing Center of Far Eastern Branch Russian Academy of Sciences, 65 Kim Yu Chen Street, Khabarovsk, 680000, Russia b Pacific National University, 136 Tikhookeanskaya Street, Khabarovsk, 680035, Russia c Changwon National University, Changwon, 51140, Republic of Korea Abstract Variational inequalities corresponding to nonlinear contact problems in mechanics often arise in engineering practice. To solve them, duality methods are widely used. As a rule, they are based on the classical methods of constructing Lagrange functionals with linear dependence in dual variables. This approach is typical for determining the saddle point – the displacement vector and normal stress in the contact area. The linear dependence in the dual variables does not allow to prove the theoretical convergence to the saddle point of the well-known iterative methods. It is possible to justify the convergence only in the primal variable under the condition that the shift in the dual variable is sufficiently small. Keywords 1 Variational inequality, contact problem, modified Lagrange functional, saddle point, Uzawa algorithm. 1. Introduction The contact problem for an elastic body with a rigid support is represented in the form of a variational inequality or an equivalent constrained minimization problem of a convex energy functional. Using the classical duality scheme, this problem can be reduced to the problem of finding a saddle point for the Lagrange functional [1-2]. Saddle point search methods for classical Lagrange functionals have been deeply and in detail investigated in many works [3-5], but as a rule, they guarantee convergence only in a primal variable. The question of convergence in a dual variable remains open. To overcome this drawback, a modified duality scheme is considered in the paper, on the basis of which a saddle point search method is constructed, which guarantees convergence in the dual variable as well. 2. Two-dimensional contact problem of elasticity Let Ξ© βŠ‚ 𝑅 2 be a bounded domain with Lipschitz boundary Ξ“. Figure 1: Elastic body with the contact zone VI International Conference Information Technologies and High-Performance Computing (ITHPC-2021), September 14–16, 2021, Khabarovsk, Russia EMAIL: rnamm@yandex.ru (A. 1); tsoy.dv@mail.ru (A. 2); vikht.el@gmail.com (A. 3); gswoo@changwon.ac.kr (A. 4) ORCID: 0000-0001-7226-3166 (A. 1); 0000-0002-4209-1284 (A. 2); 0000-0002-7152-2311 (A. 3) ©️ 2020 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) 98 We assume that on part of the boundary Ξ“D the body is rigidly fixed, on part Γ𝑁 the surface forces are given. The contact zone of an elastic body with a rigid foundation will be denoted by Γ𝐾 . 2 For the displacement vector 𝑣 = (𝑣1 , 𝑣2 ), define the strain tensor πœ€ = {πœ€π‘–π‘— }𝑖,𝑗=1 1 πœ•π‘£π‘– πœ•π‘£π‘— πœ€π‘–π‘— (𝑣) = ( + ), 𝑖, 𝑗 = 1,2, 2 πœ•π‘₯𝑗 πœ•π‘₯𝑖 2 and the stress tensor 𝜎 = {πœŽπ‘–π‘— }𝑖,𝑗=1 πœŽπ‘–π‘— (𝑣) = π‘π‘–π‘—π‘˜π‘š πœ€π‘˜π‘š (𝑣), where 𝐢 = {π‘π‘–π‘—π‘˜π‘š} is a given elasticity tensor with the usual properties of positive definiteness and symmetry π‘π‘–π‘—π‘˜π‘š = π‘π‘—π‘–π‘˜π‘š = π‘π‘˜π‘šπ‘–π‘— , 𝑖, 𝑗, π‘˜, π‘š=1,2; π‘π‘–π‘—π‘˜π‘š 𝛼𝑖𝑗 π›Όπ‘˜π‘š β‰₯ 𝑐0 𝛼𝑖𝑗 𝛼𝑖𝑗 , 𝑐0 > 0 – const. Summation over repeated indices is assumed. Let us specify vector-functions of the body and surface forces 𝑓 = (𝑓1 , 𝑓2 ) and 𝑝 = (𝑝1 , 𝑝2 ), respectively. The boundary value problem is formulated as follows βˆ’π‘‘π‘–π‘£ 𝜎(𝑒) = 𝑓 in Ξ©, (1) 𝑒 = 0 on Γ𝐷 , (2) Οƒ(𝑒)𝑛 = 𝑝 on Γ𝑁 , (3) π‘’πœˆ ≀ 0, 𝜎𝜈 (𝑒) ≀ 0, 𝜎𝜈 (𝑒)π‘’πœˆ = 0 on Γ𝐾 , (4) 𝜎𝜏 = 0 π‘œπ‘› Γ𝐾 , (5) where 𝑛 = (𝑛1 , 𝑛2 ) is the unit outward normal vector to Γ𝑁 , 𝜈 = (𝜈1 , 𝜈2 ) is the unit outward normal vector to Γ𝐾 , 𝜎𝜈 (𝑒) = πœŽπ‘–π‘— (𝑒)πœˆπ‘– πœˆπ‘— , 𝜎𝜏 (𝑒) = 𝜎(𝑒)𝜈 βˆ’ 𝜎𝜈 (𝑒) βˆ™ 𝜈, where 𝜎𝜈 (𝑒) and 𝜎𝜏 (𝑒) are normal and tangential components of the surface traction on Γ𝐾 , respectively. The boundary value problem (1)-(5) belongs to the class of problems with a free boundary since the adhesion region (π‘’πœˆ = 0) on Γ𝐾 is not known in advance and is found simultaneously with the desired solution of the problem. Condition (4) means non-penetration of an elastic body into a rigid foundation. The problem (1)-(5) has a variational formulation. Let 𝑓 ∈ [𝐿2 (Ξ©)]2 . Define the set of admissible displacements 2 𝐾 = {𝑣 ∈ [𝐻Γ1𝐷 (Ξ©)] : π‘£πœˆ ≀ 0 on Γ𝐾 }, where 𝐻Γ1𝐷 (Ξ©) = {𝑣 ∈ 𝐻1 (Ξ©): 𝑣 = 0 on Γ𝐷 }. The problem (1)-(5) corresponds to the variational inequality [6] π‘Ž(𝑒, 𝑣 βˆ’ 𝑒) β‰₯ ∫ 𝑓 β‹… (𝑣 βˆ’ 𝑒) 𝑑Ω + ∫ 𝑝 β‹… (𝑣 βˆ’ 𝑒) 𝑑Γ βˆ€ 𝑣 ∈ 𝐾, (6) Ξ© Γ𝑁 here π‘Ž(𝑒, 𝑣) = ∫Ω πœŽπ‘–π‘— (𝑣)πœ€π‘–π‘— (𝑣) 𝑑Ω = ∫Ω π‘π‘–π‘—π‘˜π‘š πœ€π‘˜π‘š (𝑣)πœ€π‘–π‘— (𝑣) 𝑑Ω. Variational inequality (6) is equivalent to the minimization problem 𝐽(𝑣) β†’ min, { (7) 𝑣 ∈ 𝐾, 1 where 𝐽(𝑣) = 2 π‘Ž(𝑒, 𝑣) βˆ’ ∫Ω 𝑓 β‹… 𝑣 𝑑Ω βˆ’ βˆ«Ξ“ 𝑝 β‹… 𝑣 𝑑Γ. 𝑁 It is known that the solution 𝑒 ∈ 𝐾 to the problem (7) exists and is unique, and it satisfies the equilibrium equation (1) and boundary conditions (2)-(5) in the generalized sense [6]. Let us formulate the dual problem for the problem (7) using the classical duality scheme. For this, we define the Lagrange functional 2 𝐿(𝑣, 𝑙) = 𝐽(𝑣) + ∫ 𝑙 π‘£πœˆ 𝑑Γ βˆ€ (𝑣, π‘˜) ∈ [𝐻Γ1𝐷 (Ξ©)] Γ— 𝐿2 (Γ𝐾 ) (8) Γ𝐾 and the corresponding dual functional πœ‘(𝑙) = inf 𝐿(𝑣, 𝑙). π‘£βˆˆ[𝐻Γ1 (Ξ©)] 2 (9) 𝐷 Problem πœ‘(𝑙) β†’ sup, { 𝑙 ∈ 𝐿+2 (Ξ“π‘˜ ), + (Ξ“ ) where 𝐿2 π‘˜ = {𝑙 ∈ 𝐿2 (Γ𝐾 ): 𝑙 β‰₯ 0 on Γ𝐾 }, is called dual to the problem (7). 99 2 A pair (𝑣 βˆ— , 𝑙 βˆ— ) ∈ [𝐻Γ1𝐷 (Ξ©)] Γ— 𝐿+2 (Γ𝐾 ) is called a saddle point of the Lagrange functional 𝐿(𝑣, 𝑙) if the following two-sided inequality takes place 2 𝐿(𝑣 βˆ— , 𝑙) ≀ 𝐿(𝑣 βˆ— , 𝑙 βˆ— ) ≀ 𝐿(𝑣, 𝑙 βˆ— ) βˆ€(𝑣, 𝑙) ∈ [𝐻Γ1𝐷 (Ξ©)] Γ— 𝐿+2 (Γ𝐾 ) In this case, 𝑣 βˆ— is the desired solution 𝑒 to the problem (7), and 𝑙 βˆ— is the solution of the dual problem (9) and coincides on Γ𝐾 with the normal stress 𝜎𝜈 (𝑒). 3. Modified Lagrange functional As already noted, the solution of the contact problem of the theory of elasticity is closely related to the search for the saddle point of the classical Lagrange functional. The well-known saddle point search algorithms for classical Lagrange functionals do not guarantee convergence in dual variables. This situation occurs, for example, in the well-known Uzawa method [3-4]. To overcome this serious drawback, let us consider a modified duality scheme that allows one to construct algorithms for finding saddle points that provide convergence in both primal and dual variables. Consider the modified Lagrange functional [7-8] 1 2 𝑀(𝑣, 𝑙) = 𝐽(𝑣) + ∫ (((𝑙 + π‘Ÿπ‘£πœˆ )+ )2 βˆ’ 𝑙 2 )𝑑Γ βˆ€ (𝑣, π‘˜) ∈ [𝐻Γ1𝐷 (Ξ©)] Γ— 𝐿2 (Γ𝐾 ). (10) 2π‘Ÿ Γ𝐾 Here (𝑙 + π‘Ÿπ‘£πœˆ )+ = π‘šπ‘Žπ‘₯{0; 𝑙 + π‘Ÿπ‘£πœˆ }, π‘Ÿ > 0 is arbitrary positive constant. For the modified functional 𝑀(𝑣, 𝑙), we define a saddle point as follows. 2 Definition. A pair (𝑣 βˆ— , 𝑙 βˆ— ) ∈ [𝐻Γ1𝐷 (Ξ©)] Γ— 𝐿2 (Γ𝐾 ) is called a saddle point of the modified Lagrange functional 𝑀(𝑣, 𝑙) if the following two-sided inequality takes place 2 𝑀(𝑣 βˆ— , 𝑙) ≀ 𝑀(𝑣 βˆ— , 𝑙 βˆ— ) ≀ 𝑀(𝑣, 𝑙 βˆ— ) βˆ€(𝑣, 𝑙) ∈ [𝐻Γ1𝐷 (Ξ©)] Γ— 𝐿2 (Γ𝐾 ). (11) The definition of the saddle point for the modified Lagrange functional 𝑀(𝑣, 𝑙) differs from the definition of the saddle point for the classical one in that in the two-sided inequality (11) the domain of variation of the dual variable 𝑙 coincides with the entire functional space 𝐿2 (Γ𝐾 ), in contrast to the corresponding inequality for the classical analogue, whereas the domain variation of the dual variable 𝑙 is taken by 𝐿+2 (Γ𝐾 ). Despite this, the sets of saddle points for 𝐿(𝑣, 𝑙) and 𝑀(𝑣, 𝑙) coincide. This important property of the modified Lagrange functionals is provided by the nonlinear dependence of the 𝑀(𝑣, 𝑙) on the dual variable in formula (10). Introduce the dual functional 𝑀(𝑙) = inf 2 𝑀(𝑣, 𝑙) π‘£βˆˆ[𝐻Γ1 (Ξ©)] 𝐷 and the corresponding dual problem 𝑀(𝑙) β†’ sup, { (12) 𝑙 ∈ 𝐿2 (Γ𝐾 ). The following statement holds [7, 8]. Theorem. The dual functional 𝑀(𝑙) is Gateaux differentiable in 𝐿2 (Γ𝐾 ) and its derivative βˆ‡π‘€(𝑙) satisfies the Lipschitz condition with the constant 1/π‘Ÿ, that is 1 β€–βˆ‡π‘€(𝑙1 ) βˆ’ βˆ‡π‘€(𝑙2 ) ‖𝐿 (Ξ“ ) ≀ ‖𝑙1 βˆ’ 𝑙2 ‖𝐿2 (Γ𝐾 ) βˆ€ 𝑙1 , 𝑙2 ∈ 𝐿2 (Ξ“π‘˜ ). 2 𝐾 π‘Ÿ It can be shown that βˆ‡π‘€(𝑙) = max{π‘’πœˆ , βˆ’π‘™/π‘Ÿ} βˆ€ 𝑙 ∈ 𝐿2 (Ξ“π‘˜ ), where 1 𝑒 = arg min 2 {𝐽(𝑣) + ∫ (((𝑙 + π‘Ÿπ‘£πœˆ )+ )2 βˆ’ 𝑙 2 )𝑑Γ} . π‘£βˆˆ[𝐻 1 (Ξ©)] 2π‘Ÿ Γ𝐾 Γ𝐷 To solve the dual problem (12), taking into account the above theorem, we can consider the gradient method 𝑙 π‘˜+1 = 𝑙 π‘˜ + π‘Ÿ βˆ‡π‘€(𝑙 π‘˜ ), π‘˜ = 1,2, … with any initial 𝑙 0 ∈ 𝐿2 (Γ𝐾 ). 100 Since βˆ‡π‘€(𝑙 π‘˜ ) = max{π‘’πœˆπ‘˜+1 , βˆ’π‘™ π‘˜ /π‘Ÿ}, the gradient method is transformed into the following Uzawa algorithm 1 + 2 2 (𝑖) π‘’π‘˜+1 =arg min 2 {𝐽(𝑣) + ∫ (((𝑙 π‘˜ + π‘Ÿπ‘£πœˆ ) ) βˆ’ (𝑙 π‘˜ ) ) 𝑑Γ} , π‘£βˆˆ[𝐻 1 (Ξ©)] 2π‘Ÿ Γ𝐾 (13) Γ𝐷 + (𝑖𝑖) 𝑙 π‘˜+1 = (𝑙 π‘˜ + π‘Ÿ π‘’πœˆπ‘˜+1 ) . Under the condition of solvability of the dual problem (12), it is possible to prove the weak convergence of the sequence {𝑙 π‘˜ } generated by the Uzawa algorithm to the solution 𝑙 βˆ— of the dual problem. In this case, the sequence {π‘’π‘˜ } converges to the desired solution π‘’βˆ— with respect to the minimized functional, that is lim 𝐽(π‘’π‘˜ ) = 𝐽(π‘’βˆ— ) [7]. It can be proved, if the sequence {π‘’π‘˜ } belongs to π‘˜β†’βˆž [𝐻 2 (Ξ©)]2 and is bounded, then the sequence {(π‘’π‘˜ , 𝑙 π‘˜ )} converges to (π‘’βˆ— , 𝑙 βˆ— ) at the norm of 2 [𝐻Γ1𝐷 (Ξ©)] Γ— 𝐿2 (Γ𝐾 ) and, at the same time, 𝑙 βˆ— = βˆ’πœŽπœˆ (π‘’βˆ— ). 4. Numerical solution of the contact problem of elasticity Research on the numerical analysis of variational inequalities is carried out, as a rule, on the basis of the finite element method [2, 3]. Let Ξ© βŠ‚ 𝑅 2 be a bounded polygonal domain. Let us carry out a finite element approximation of problem (7) using piecewise bilinear basis functions [9]. By standard transformations, problem (7) is transformed into a finite-dimensional quadratic programming problem 1 π’₯(π‘₯) = < 𝐴π‘₯, π‘₯ > βˆ’ βˆ‘ (𝑓𝑗1 π‘₯𝑗1 + 𝑓𝑗2 π‘₯𝑗2 ) βˆ’ βˆ‘ (𝑝𝑗1 π‘₯𝑗1 + 𝑝𝑗2 π‘₯𝑗2 ) β†’ min , 2 (14) π‘—βˆˆπ’© π‘—βˆˆβ„³ { 𝑙(π‘₯ β‹… 𝜈)𝑗 ≑ π‘₯𝑗1 𝜈1 + π‘₯𝑗2 𝜈2 ≀ 0 𝑗 ∈ 𝒫. where 𝒩 is the set of indices of quadrilateral mesh nodes, |𝒩| is the cardinality of the set 𝒩, β„³ is the set of indices of mesh nodes on Γ𝑁 , 𝒫 is the set of mesh nodes indices on Γ𝐾 , 𝐴 = (π‘Žπ‘–π‘— ), 𝑖, 𝑗 = 1, … ,2|𝒩| is the stiffness matrix, (𝑓𝑗1 , 𝑓𝑗2 ), (𝑝𝑗1 , 𝑝𝑗2 ) are the coordinates of the expansion of the vectors of volume and surface forces in the finite element basis of each node 𝑗, respectively. Under the natural condition that π‘šπ‘’π‘Žπ‘  Γ𝐷 > 0, the matrix 𝐴 is symmetric and positive definite. The function to be minimized in problem (14) corresponds to a finite-dimensional approximation of the functional 𝐽(𝑣) in problem (7). To solve the problem (7), we use the Uzawa algorithm with a modified Lagrange functional in a finite-dimensional version. Let us apply one of the quadrature formulas for the finite-dimensional approximation of the expression 1 ∫ (((𝑙 + π‘Ÿπ‘£πœˆ )+ )2 βˆ’ 𝑙 2 )𝑑Γ. 2π‘Ÿ Γ𝐾 As a result, we obtain a continuously differentiable piecewise-quadratic function in the variables π‘₯𝑗1 , π‘₯𝑗2 , 𝑗 ∈ 𝒫, of the form 1 + 2 2 βˆ‘ (((𝑙𝑗 + π‘Ÿ(π‘₯ β‹… 𝜈)𝑗 ) ) βˆ’ (𝑙𝑗 ) ) β„Žπ‘— , 2π‘Ÿ π‘—βˆˆπ’« where 𝑙𝑗 , β„Žπ‘— are known quantities. Let us set an arbitrary 𝑙 0 ∈ 𝑅 |𝒫| , where |𝒫| is the cardinality of the set 𝒫. Uzawa algorithm (13) in the finite-dimensional case has the form 1 + 2 2 (𝑖)β€² π‘₯ π‘˜+1 =arg min {π’₯(π‘₯) + βˆ‘ (((π‘™π‘—π‘˜ + π‘Ÿ(π‘₯ β‹… 𝜈)𝑗 ) ) βˆ’ (π‘™π‘—π‘˜ ) ) β„Žπ‘— } , 2|𝒩| π‘₯βˆˆπ‘… 2π‘Ÿ (15) π‘—βˆˆπ’« + (𝑖𝑖)β€²β€² π‘™π‘—π‘˜+1 = (π‘™π‘—π‘˜ + π‘Ÿ(π‘₯ β‹… 𝜈)𝑗 ) , 𝑗 ∈ 𝒫. Let us consider the step (𝑖)β€². It is the problem of minimizing a continuously differentiable piecewise quadratic function. A feature of the function to be minimized is that its Hessian has discontinuities on 101 some linear manifolds and, at the same time, the gradient of the function is continuous. To minimize function (16) we apply a natural generalization of Newton’s method [10], [12] 1 + 2 2 π’₯(π‘₯) + βˆ‘ (((π‘™π‘—π‘˜ + π‘Ÿ(π‘₯ β‹… 𝜈)𝑗 ) ) βˆ’ (π‘™π‘—π‘˜ ) ) β„Žπ‘— (16) 2π‘Ÿ π‘—βˆˆπ’« under a fixed 𝑙 π‘˜ . Let us assume that at some step π‘š the monotony of generalized Newton method is broken, i.e. 1 + 2 2 π’₯(π‘₯ π‘šβˆ’1 ) + βˆ‘ (((π‘™π‘—π‘˜ + π‘Ÿ(π‘₯ π‘šβˆ’1 β‹… 𝜈)𝑗 ) ) βˆ’ (π‘™π‘—π‘˜ ) ) < 2π‘Ÿ π‘—βˆˆπ’« 1 + 2 2 < π’₯(π‘₯ π‘š ) + βˆ‘ (((π‘™π‘—π‘˜ + π‘Ÿ(π‘₯ π‘š β‹… 𝜈)𝑗 ) ) βˆ’ (π‘™π‘—π‘˜ ) ). 2π‘Ÿ π‘—βˆˆπ’« + Let 𝒫̃ βŠ‚ 𝒫 and be such that (π‘™π‘—π‘˜ + π‘Ÿ(π‘₯ π‘šβˆ’1 β‹… 𝜈)𝑗 ) = π‘™π‘—π‘˜ + π‘Ÿ(π‘₯ π‘šβˆ’1 β‹… 𝜈)𝑗 for all 𝑗 ∈ 𝒫̃ . Consider the convex quadratic function 1 2 2 π’₯(π‘₯) + βˆ‘ ((π‘™π‘—π‘˜ + π‘Ÿ(π‘₯ β‹… 𝜈)𝑗 ) βˆ’ (π‘™π‘—π‘˜ ) ). (17) 2π‘Ÿ π‘—βˆˆπ’«Μƒ Its Hessian and gradient will be denoted as π΄π‘šβˆ’1 and π΄π‘šβˆ’1 π‘₯ βˆ’ π΅π‘šβˆ’1 , respectively, where π΅π‘šβˆ’1 is some vector. Taking into account that according to Newton’s method π‘₯ π‘š = π‘₯ π‘šβˆ’1 βˆ’ π΄βˆ’1 π‘šβˆ’1 (π΄π‘šβˆ’1 π‘₯ π‘šβˆ’1 βˆ’ π΅π‘šβˆ’1 ), then βˆ’βŒ©(π΄π‘šβˆ’1 π‘₯ π‘šβˆ’1 βˆ’ π΅π‘šβˆ’1 ), π‘₯ π‘š βˆ’ π‘₯ π‘šβˆ’1 βŒͺ = βŒ©π΄π‘šβˆ’1 (π‘₯ π‘š βˆ’ π‘₯ π‘šβˆ’1 ), π‘₯ π‘š βˆ’ π‘₯ π‘šβˆ’1 βŒͺ > 0 or βŒ©π΄π‘šβˆ’1 π‘₯ π‘šβˆ’1 βˆ’ π΅π‘šβˆ’1 , π‘₯ π‘š βˆ’ π‘₯ π‘šβˆ’1 βŒͺ < 0. Thus, in the vicinity of the point π‘₯ π‘šβˆ’1 , the minimized function (16) decreases locally in the direction π‘₯ π‘š βˆ’ π‘₯ π‘šβˆ’1. The point π‘₯ π‘š is the minimum point of the quadratic function (17), but is not the desired minimum point of the convex function (16). Using the well-known rule of Armijo [11] we find a number π›Όπ‘š > 0 such that the value of the minimized function (17) at the element π‘₯ π‘šβˆ’1 + π›Όπ‘š (π‘₯ π‘š βˆ’ π‘₯ π‘šβˆ’1 ) will be less that at the element π‘₯ π‘šβˆ’1. Next, we take π‘₯ π‘šβˆ’1 + π›Όπ‘š (π‘₯ π‘š βˆ’ π‘₯ π‘šβˆ’1 ) as a new π‘₯Μƒ π‘š and return to the usual Newtonian method. We ensured monotonicity in the process of minimizing the piecewise quadratic continuously differentiable function (16). This ensures the convergence in the minimized function of the generalized Newton’s methods, and from strong convexity, the convergence in the arguments π‘₯ in a finite number of steps (due to the finiteness of the set 𝒫). The criterion that the minimum point of function (16) is found is the repeatability of the set 𝒫̃ at two successive steps of the generalized Newton method. The use of the generalized Newton method at step (𝑖)β€² of the Uzawa algorithm with a modified Lagrange function significantly accelerates the search for a solution to the problem (14). As a rule, the monotonicity of the minimization process at the (𝑖)β€² step is achieved automatically after the first or second application of the Armijo rule. Overall, the Uzawa algorithm (15) rapidly converges to a saddle point due to the fast stabilization of the sequence {𝑙 π‘˜ }. Note that, within the framework of solving the linear programming problem, several authors previously investigated similar algorithms for minimizing piecewise quadratic functions [12]. 5. Numerical experiments We shall consider the following model example: the body Ω = (0, 3) Γ— (0, 1) (in m) is made of an elastic isotropic, homogeneous material characterized by Young’s modulus E = 21.19x104 MPa and Poisson’s ratio ΞΌ = 0.277. It is fixed along Γ𝐷 = {0} Γ— (0, 1) and linearly distributed surface tractions of density 𝑝 = (𝑝1 , 𝑝2 ) are applied of Γ𝑁 = Γ𝑁1 βˆͺ Γ𝑁2 , where Γ𝑁1 = {3} Γ— (0, 1) and Γ𝑁2 = (0, 3) Γ— {1}. We consider the following traction forces: 𝑝1 (x) = 0, 𝑝2 (x) = 1 MPa on Γ𝑁1 , 𝑝1 (x) = 0, 𝑝2 (x) = βˆ’ 2⁄3 (3 βˆ’ x1 ) MPa on Γ𝑁2 , parameter π‘Ÿ = 108 . 102 The body is discretized into 𝑁π‘₯ Γ— 𝑁𝑦 4-node quadrilateral finite elements, where 𝑁π‘₯ (= 3𝑁𝑦 ) is varied to generate problem instances with different sizes. The number of degrees of freedom of displacements is 𝑛𝑝 = 2(𝑁π‘₯ + 1)(𝑁𝑦 + 1) and the number of contact candidate nodes is 𝑛𝑑 = 𝑁π‘₯ + 1. We assume the small deformation and solve the finite-dimensional optimization problem using the generalized Newton method (GNM). Finding an inverse matrix is a very computationally expensive operation. Therefore, rather than computing the inverse of the generalized Hessian matrix, one may save time and increase numerical stability by solving the system of linear equations. For this purpose, we use the conjugate gradient method implemented in the SciPy [13] and CuPy [14] packages and compare the computation time. In all numerical experiments considered below, we choose the following condition: β€–π‘₯ π‘š βˆ’ π‘₯ π‘šβˆ’1 β€–2 max (‖𝐺(π‘₯ π‘š )β€–2 , ) < 10βˆ’10 β€–π‘₯ π‘š β€–2 as a stopping criterion for the generalized Newton method. Uzawa algorithm terminates if ‖𝑙 π‘˜ βˆ’ 𝑙 π‘˜βˆ’1 β€–2 < 10βˆ’8 . ‖𝑙 π‘˜ β€–2 All experiments are implemented in Python, using the scikit-fem library [15] for performing finite element assembly and CuPy library for GPU-accelerated computing. Computation was carried out on IBM Power Systems S822LC 8335-GTB server, which is based on two 10-core IBM POWER8 processors with a maximum operating frequency of 4.023 GHz and two NVIDIA Tesla P100 GPU accelerators. Table 1 shows how the total number of Uzawa method iterations and the number of generalized Newton method iterations depend on 𝑛𝑝 ΠΈ 𝑛𝑑 . The number of GNM iterations on the first step of the Uzawa method is represented by the first integer in the respective column. The second integer represents the number of iterations in subsequent steps of the Uzawa method. We can see that the number of iterations slightly increases with the increasing number of primal and dual variables. Calculations show that the use of CuPy library can speed up the execution time for problems of large size up to 12.6 times. Table 1 Computational results 𝑁π‘₯ Γ— 𝑁𝑦 𝑛𝑝 𝑛𝑑 GNM it Uzawa it Time CPU, s Time GPU, s M(uβˆ— , lβˆ— ) 60x20 2562 61 7/2 6 1.89 5.06 -6.700472e-05 120x40 9922 121 8/2 7 7.70 8.76 -6.711822e-05 240x80 39042 241 9/2 8 56.35 17.24 -6.714898e-05 480x160 154882 481 10/2 11 480.60 38.06 -6.715744e-05 Figure 2 shows normal and tangential displacements of the body on the contact zone for 𝑛𝑑 = 481. Figure 2: Normal and tangential contact displacements 103 The value of the dual variable (normal contact stress) is depicted in Figure 3. We see that the penetration of the elastic body into a rigid foundation does not occur and at the points where the body is in contact, the dual variable is positive. The resulting domain deformation in Lagrange coordinates x+1000u(x) with an amplification factor 1000 and Von Mises stresses are presented in Figure 4. Figure 3: Normal contact stress Figure 4: Distribution of the von Mises stresses 6. Conclusion In the paper, the numerical algorithm of solving the contact problem was proposed. The algorithm based on the modified Lagrange functionals and Uzawa method. The algorithm was implemented by using the finite element method. The numerical experiments illustrating the fast convergence of the algorithm by primal and dual variables were presented. This circumstance can be explained by the good differential properties of the modified dual functional, which makes it possible to implement the gradient method for solving the dual problem. After a discretization, we used the generalized Newton method with Armijo line search to solve the minimization problem of piecewise quadratic functional. 7. Acknowledgements This study was supported by the Russian Foundation for Basic Research (Project 20-01-00450 A), by the Ministry of Science and Higher Education of the Russian Federation, Supplementary Agreement No. 075-02-2020-1529/1 dated April 21, 2020. 104 The studies were carried out using the resources of the Center for Shared Use of Scientific Equipment "Center for Processing and Storage of Scientific Data of the Far Eastern Branch of the Russian Academy of Sciences", funded by the Russian Federation represented by the Ministry of Science and Higher Education of the Russian Federation under project No. 075-15-2021-663. 8. References [1] I. Hlavačhek, Y. Haslinger, I. Nechas, Y. LoviΕ‘hek, Solution of Variational Inequalities in Mechanics, Springer-Verlag, New York, 1986. [2] N. Kikuchi, T. Oden, Contact Problem in Elasticity: a Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, 1988. [3] R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer, New York, 1984. [4] I. Ekland, R. Temam, Convex Analysis and Variational Problems, North-Holland, Amsterdam, 1976. [5] R. Glowinski R., J. L. Lions, R. Tremolier, Numerical Analysis of Variational Inequalities, North- Holland, Amsterdam, 1981. [6] A. M. Khludnev, Elasticity Problems in Nonsmooth Domains, Fizmatlit, Moscow, 2010. [7] R. V. Namm, G. I. Tsoy, A modified dual scheme for solving an elastic crack problem, Numer. Anal. Appl. 10(1) (2017) 37–46. doi.org/10.1134/S1995423917010050. [8] R. V. Namm, E. M. Vikhtenko, N. N. Maksimova, Sensitivity functionals in variational inequalities of mechanics and their applications to duality schemes, Numer. Anal. Appl. 7(1) (2014) 36–44. doi.org/10.1134/S1995423914010042. [9] P. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, Amsterdam, 1978. [10] O. L. Mangasarian, A generalized Newton method for absolute value equations, Optim. Lett. 3(1) (2009) 101-108. doi.org/10.1007/s11590-008-0094-5. [11] L. Armijo, Minimization of functions having Lipschitz continuous first partial derivatives, Pacific J. Math. 16(1) (1966) 1-3. doi.org/10.2140/pjm.1966.16.1. [12] A. I. Golikov, Yu. G. Evtushenko, Generalized Newton’s method for linear optimization problems with inequality constraints, Proc. Steklov Inst. Math. 284 (2014) 96–107. doi.org/10.1134/S0081543814020096. [13] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, et al., SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Methods 17(3) (2020) 261– 272. doi.org/10.1038/s41592-019-0686-2. [14] R. Okuta, Y. Unno, D. Nishino, S. Hido, and C. Loomis, CuPy: A NumPy-Compatible Library for NVIDIA GPU Calculations, in: Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems(NIPS), ), Long Beach, CA, USA, 2017. [15] T. Gustafsson, G. D. McBain, scikit-fem: A Python package for finite element assembly, J. Open Source Softw. 5(52) (2020) 2369. doi.org/10.21105/joss.02369. 105