Variational Method for Solving Contact Problem of Elasticity

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. Variational inequality, contact problem, modified Lagrange functional, saddle point, Uzawa 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.

1. Introduction

2020 Copyright for this paper by its authors. are given. The contact zone of an elastic body with a rigid foundation will be denoted by Γ .

We assume that on part of the boundary ΓD the body is rigidly fixed, on part Γ the surface forces For the displacement vector = ( 1, 2), define the strain tensor = { 2 } , =1 ( ) = 2 } , =1 1 +

( 2 ( ) = ) ,

, = 1,2, ( ), and the stress tensor = { where = { symmetry

} is a given elasticity tensor with the usual properties of positive definiteness and = = , , , , =1,2; ≥ 0 , 0 > 0 – const. respectively. The boundary value problem is formulated as follows

Summation over repeated indices is assumed.

Let us specify vector-functions of the body and surface forces = ( 1, 2) and = ( 1, 2), −

( ) = in Ω, = 0 on Γ , σ( ) = on Γ , = 0

Γ , ≤ 0,

( ) ≤ 0, ( ) = 0 on Γ , displacements 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 = { ∈ [ Γ1 (Ω)] : ≤ 0 on Γ },

2 where Γ1 (Ω) = { ∈ 1(Ω): = 0 on Γ }.

The problem ( 1 )-( 5 ) corresponds to the variational inequality [6] ( , − ) ≥ ∫ ⋅ ( − ) Ω + ∫ ⋅ ( − ) Γ ∀ ∈ , here ( , ) = ∫Ω ( ) ( ) Ω = ∫Ω

Variational inequality ( 6 ) is equivalent to the minimization problem Ω Γ

( ) ( ) Ω. { ( ) → min, ∈ ,

inf

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, ( , ) = ( )+ ∫ Γ ∀ ( , ) ∈ [ Γ1 (Ω)] × 2(Γ ) 2 we define the Lagrange functional and the corresponding dual functional where +2(Γ ) = { ∈ 2(Γ ): ≥ 0 on Γ }, is called dual to the problem ( 7 ). the following two-sided inequality takes place A pair ( ∗, ∗) ∈ [ Γ1 (Ω)] × +2(Γ ) is called a saddle point of the Lagrange functional ( , ) if 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 ( ).

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

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 ( ∗, ) ≤ ( ∗, ∗) ≤ ( , ∗) ∀( , ) ∈ [ Γ1 (Ω)] × 2(Γ ).

2

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 and the corresponding dual problem ( ) =

inf ( 11 ) ( 12 ) 2(Γ ) ≤ 1 ‖∇ ( 1)− ∇ ( 2)‖

‖ 1 − 2 ‖ 2(Γ ) ∀ 1, 2 ∈ 2(Γ ).

It can be shown that ∇ ( ) = max{ , − / } ∀ ∈ 2(Γ ), where = arg

min with any initial 0 ∈ 2(Γ ).

+1 = + ∇ ( ), = 1,2, …

Since ∇ ( ) = max{ +1, − / }, the gradient method is transformed into the following Uzawa algorithm ( ) +1=arg ( )

min ∫ ((( + ) ) − ( )2) Γ}, +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 ( )′ ( )′′ +1=arg ∈ 2| | min { ( )+ 1 2 ∑ ((( + ( ⋅ ) ) ) − (

)2)ℎ } , + 2 +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

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 ∈ . 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 As a result, we obtain a continuously differentiable piecewise-quadratic function in the variables 1, 2, ∈ , of the form 1 ∫ ((( + )+)2 − 2) Γ. 1 2 ∑ ((( + ( ⋅ ) ) ) − ( ) )ℎ , where , ℎ are known quantities. the finite-dimensional case has the form

Let us set an arbitrary 0 ∈ | |, where | | is the cardinality of the set . Uzawa algorithm ( 13 ) in ( 13 ) ( 14 ) ( 15 ) 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) under a fixed .

Let us assume that at some step the monotony of generalized Newton method is broken, i.e. ( −1)+ Let ̃ ⊂ and be such that ( + ( −1 ⋅ ) )+ = + ( −1 ⋅ ) for all ∈ ̃ . Consider the convex quadratic function

∈ ̃ 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), ( )+ ∑ (( + ( ⋅ ) )2 − ( )2). (17) then or −〈( −1 −1 − −1), − −1〉 = 〈 −1( − −1), − −1〉 > 0

〈 −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]. 2 1 2

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.

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) < 10−10

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

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.

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.

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

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