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][4][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.

Two-dimensional contact problem of elasticity

Let Ω ⊂ 𝑅 2 be a bounded domain with Lipschitz boundary Γ. 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 Γ 𝐾 .

For the displacement vector 𝑣 = (𝑣 1 , 𝑣 2 ), define the strain tensor 𝜀 = {𝜀 𝑖𝑗 } 𝑖,𝑗=1 𝜎 𝑖𝑗 (𝑣) = 𝑐 𝑖𝑗𝑘𝑚 𝜀 𝑘𝑚 (𝑣), where 𝐶 = {𝑐 𝑖𝑗𝑘𝑚} is a given elasticity tensor with the usual properties of positive definiteness and symmetry 𝑐 𝑖𝑗𝑘𝑚 = 𝑐 𝑗𝑖𝑘𝑚 = 𝑐 𝑘𝑚𝑖𝑗 , 𝑖, 𝑗, 𝑘, 𝑚=1,2; 𝑐 𝑖𝑗𝑘𝑚 𝛼 𝑖𝑗 𝛼 𝑘𝑚 ≥ 𝑐 0 𝛼 𝑖𝑗 𝛼 𝑖𝑗 , 𝑐 0 > 0const.

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

where

𝐽(𝑣) = 1 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

𝐿(𝑣, 𝑙) = 𝐽(𝑣) + ∫ 𝑙 𝑣 𝜈 𝑑Γ Γ 𝐾 ∀ (𝑣, 𝑘) ∈ [𝐻 Γ 𝐷 1 (Ω)] 2 × 𝐿 2 (Γ 𝐾 )(8)

and the corresponding dual functional 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 𝜎 𝜈 (𝑢).

𝜑(𝑙) = inf 𝑣∈[𝐻 Γ 𝐷 1 (Ω)] 2 𝐿(𝑣, 𝑙).(9)

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 × 𝐿 2 (Γ 𝐾 ).(10)

Here (𝑙 + 𝑟𝑣 𝜈 ) + = 𝑚𝑎𝑥{0; 𝑙 + 𝑟𝑣 𝜈 }, 𝑟 > 0 is arbitrary positive constant.

For the modified functional 𝑀(𝑣, 𝑙), we define a saddle point as follows.

Definition. A pair (𝑣 * , 𝑙 * ) ∈ [𝐻 Γ 𝐷 1 (Ω)] 2 × 𝐿 2 (Γ 𝐾

) is called a saddle point of the modified Lagrange functional 𝑀(𝑣, 𝑙) if the following two-sided inequality takes place

𝑀(𝑣 * , 𝑙) ≤ 𝑀(𝑣 * , 𝑙 * ) ≤ 𝑀(𝑣, 𝑙 * ) ∀(𝑣, 𝑙) ∈ [𝐻 Γ 𝐷 1 (Ω)] 2 × 𝐿 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 𝑣∈[𝐻 Γ 𝐷 1 (Ω)] 2 𝑀(𝑣, 𝑙)

and the corresponding dual problem

{ 𝑀(𝑙) → sup, 𝑙 ∈ 𝐿 2 (Γ 𝐾 ).(12)

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 ) − ∇𝑀(𝑙 2 ) ‖ 𝐿 2 (Γ 𝐾 ) ≤ 1 𝑟 ‖𝑙 1 − 𝑙 2 ‖ 𝐿 2 (Γ 𝐾 ) ∀ 𝑙 1 , 𝑙 2 ∈ 𝐿 2 (Γ 𝑘 ).

It can be shown that ∇𝑀(𝑙) = max{𝑢 𝜈 , −𝑙/𝑟} ∀ 𝑙 ∈ 𝐿 2 (Γ 𝑘 ), where

𝑢 = arg min 𝑣∈[𝐻 Γ 𝐷 1 (Ω)] 2 {𝐽(𝑣) + 1 2𝑟 ∫ (((𝑙 + 𝑟𝑣 𝜈 ) + ) 2 − 𝑙 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 (Γ 𝐾 ).

Since ∇𝑀(𝑙 𝑘 ) = max{𝑢 𝜈 𝑘+1 , −𝑙 𝑘 /𝑟}, the gradient method is transformed into the following Uzawa algorithm

(𝑖) 𝑢 𝑘+1 =arg min 𝑣∈[𝐻 Γ 𝐷 1 (Ω)] 2 {𝐽(𝑣) + 1 2𝑟 ∫ (((𝑙 𝑘 + 𝑟𝑣 𝜈 ) + ) 2 − (𝑙 𝑘 ) 2 ) 𝑑Γ Γ 𝐾 } , (𝑖𝑖) 𝑙 𝑘+1 = (𝑙 𝑘 + 𝑟 𝑢 𝜈 𝑘+1 ) + .(13)

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 and is bounded, then the sequence {(𝑢 𝑘 , 𝑙 𝑘 )} converges to (𝑢 * , 𝑙 * ) at the norm of

[𝐻 2 (Ω)][𝐻 Γ 𝐷 1 (Ω)] 2 × 𝐿 2 (Γ 𝐾

) and, at the same time, 𝑙 * = −𝜎 𝜈 (𝑢 * ).

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 2 < 𝐴𝑥, 𝑥 > − ∑ (𝑓 𝑗1 𝑥 𝑗1 + 𝑓 𝑗2 𝑥 𝑗2 ) − 𝑗∈𝒩 ∑ (𝑝 𝑗1 𝑥 𝑗1 + 𝑝 𝑗2 𝑥 𝑗2 ) → min 𝑗∈ℳ , 𝑙(𝑥 ⋅ 𝜈) 𝑗 ≡ 𝑥 𝑗1 𝜈 1 + 𝑥 𝑗2 𝜈 2 ≤ 0 𝑗 ∈ 𝒫. (14)

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)

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 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] (

)17

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].

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.19x10 4 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 − x 1 ) MPa on Γ 𝑁 2 , parameter 𝑟 = 10 8 .

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:

max (‖𝐺(𝑥 𝑚 )‖ 2 , ‖𝑥 𝑚 − 𝑥 𝑚−1 ‖ 2 ‖𝑥 𝑚 ‖ 2 ) < 10 −10

as a stopping criterion for the generalized Newton method. Uzawa algorithm terminates if

‖𝑙 𝑘 − 𝑙 𝑘−1 ‖ 2 ‖𝑙 𝑘 ‖ 2 < 10 −8 .

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. 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.

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.