In this paper, boundary element and augmented Lagrangian methods for Coulomb friction contact problems are presented. Based on the projection technique, both unilateral contact and Coulomb friction conditions are reformulated as fixed point problems. The original problem is deduced to a variational formulation with boundary integral operators. Then, we propose a new augmented Lagrangian method which can be dealt with the semismooth Newton method. Short theoretical results and the algorithm description are given. Numerical simulations show the performance of the method proposed.

1. Introduction

In many industrial applications or engineering problems, the contact between a deformable elastic body and a rigid obstacle plays an important role. Due to the nonlinear boundary condition, these problems are rather complicated from both the theoretical and numerical points of view. Therefore, reliable and efficient methods for the numerical simulation of friction contact problems are quite necessary for many areas of solid mechanics. For the existence and uniqueness of a solution, the weak formulation of the considered problem leads to a variational inequality which is a traditional method for the problem [14]. During the past decades, a vast majority of methods developed have been based on the variational inequalities. We note that there exist two types of approaches for the numerical solution of the contact problem. An option is to discrete the problem by the finite element method (FEM) [25] or the boundary element method (BEM) [610] and obtain a convex optimization problem in the finite dimensional space. Another option is to use the Lagrange multiplier which replaces the nonlinear problem with a sequence of linear problems in function spaces, and this idea has been introduced in [1118]. Recently, some new methods for the numerical simulation of the friction contact problem have also been developed, and we mention the penalty method and Nitsche’s method [1927].

Fixed point methods such as projection techniques are also widely applied to complementary problems including contact problems in linear elasticity [1214, 22, 25, 26]. The basic idea of these methods is to rewrite the inequality system as an equivalent equality system by using the projection. Although it is still nonlinear and not differentiable in the classical sense, semismooth Newton methods can be applied and easily implemented in terms of a primal-dual active set strategy. The main advantage of these methods is that the generalized derivatives of the nonlinear projection function are extremely easy to compute. This is quite important to deal with the nonlinearity of the problem. We refer to the series of recent papers by Stadler et al. [12, 13] and Chouly and Hild [22, 23], where the augmented Lagrangian method (ALM) and the semismooth Newton (SSN) methods have been applied to different contact problems.

The BEM has been also employed to solve the contact problem in elasticity. The main benefit of BEM is the significant reduction of expense mesh generation because the formulation of the problem is reduced to the boundary of the domain. Consequently, the boundary variational formulation expresses in terms of boundary unknowns only. For contact problems, the key unknowns are displacement and stress on the contact boundary, which are considered as primary variables and can be obtained directly in the BEM [2830]. Therefore, the BEM seems to be the natural way for these problems [610, 3137].

In this paper, we focus on the combination of augmented Lagrangian methods and boundary element methods for the Coulomb friction contact problem. To be precise, we use the projection technique to deal with the nonlinear boundary conditions and obtain the corresponding equivalent formulations. Then, the weak formulation of the original problem is deduced by using boundary integral operators. When compared with other methods, there is no inequality constraint, and a boundary integral equation is introduced, which is useful from both theoretical and numerical points of view. At each iteration step, our method only needs updating boundary values and solving a boundary variational problem, and we can apply the semismooth Newton method for the solution [11, 12, 15]. Using properties of projection and the Steklov–Poincaré operator (also known as the Dirichlet-to-Neumann mapping), we obtain the convergence in function spaces.

The structure of this paper is as follows. In Section 2, we introduce two equivalent formulations to replace the nonlinear condition, for the classical Coulomb frictional contact problems. A weak formulation with the boundary integral operators is presented in Section 3. Section 4 is devoted to an ALM with unconditional monotone convergence for all penalty parameters. In Section 5, we give two numerical examples to illustrate the behavior of the method, and some conclusions and perspectives are drawn in the last section.

2. Setting of the Problem

We consider the classical Coulomb frictional contact problem with a rigid foundation. Let be an open and bounded domain with a Lipschitz boundary . The boundary consists of three mutually disjoint parts , , and , where Dirichlet, Neumann, and frictional contact conditions are prescribed. Suppose that there are no volume forces acting on the body. The case of nonvanishing volume forces can be treated using Newton potentials. Let and be the outward normal and tangential vectors of , respectively. For a given boundary traction and obstacle , the problem consists in finding the displacement satisfyingwhere is the stress tensor, is the linearized strain tensor, and is the fourth-order elasticity tensor which satisfies usual conditions of symmetry, coercivity, and boundedness.

On , we use decomposition of the displacement and the stress vector fields in normal and tangential as follows:

For a given obstacle , the unilateral contact condition is given byand the Coulomb friction condition reads aswhere is the friction coefficient.

Let us introduce the following Hilbert spaces:and notations:

In the above definition, for any , the product is defined as .

Problems (1)–(8) are then equivalent to the following variational inequality:or a convex minimization problem:

It has been proved that problem (11) or equivalently (12) has a unique solution in the theory of variational inequalities [1, 13, 23].

For this problem, the major difficulty is due to the coupling between the contact condition and the friction condition. Moreover, the transition between contact and noncontact is characterized by a change in the type of the boundary condition with constraint. In this paper, we rewrite the friction contact conditions (6)–(8) as equivalent projection formulas. Let us introduce the projection notation for :and the projection notation for , :

As a result, we obtain the following results for the contact and friction conditions on .

Lemma 1. For all , the contact conditions (6) on are equivalent to

The proof can be found in [22, 3537]. For the friction condition, we also have the following result according to [23, 37].

Lemma 2. For all and , the contact conditions (7) and (8) on are equivalent to

Proof. Firstly, let us assume that (7) and (8) hold. Let , and then from (7), we haveFor the case , from (8), we haveIt results that , which means thatNote that and have the same orientation. Then, using (8) yieldsOn the other hand, let and in (16) hold. We first assume that , that is, . Thus, (16) isso that and . Then, (7) is satisfied. Finally, we consider the case , which means that and there exists such that . Then, from (16), we haveso that in this case,This implies, due to and relationship (16), that , , and . Hence, the friction conditions (7) and (8) hold.

3. Boundary Weak Formulation of the Frictional Contact Problem

In order to give a boundary variational formulation to the frictional contact problem, we apply Green’s formula to (1)–(8) and obtain the following variational formulation:

Let be the fundamental solution of the two-dimensional Lamé equation:where Lamé constants and are given by Young’s modulus and Poisson ratio asand is the unit matrix. We then introduce boundary integral operators: the single layer potential , the double layer potential , the adjoint double layer potential , and the hypersingular integral operator :where is the boundary traction operator defined by . Now, we can define the Steklov–Poincaré operator by

Properties of the boundary integral operators yield that the Steklov–Poincaré operator is linear, symmetric, and positive definite on [3032], i.e., there exists a positive constant such that for any ,

Furthermore, we note the Dirichlet-to-Neumann mapping on :and , and then we get

Let us define

By substituting (31) into (24), we use and obtain the boundary weak formulation of (1)–(4):

Consequently, we obtain boundary weak formulation (33) and projection fixed point problems (15) and (16), which are equivalent to the original problems (1)–(8) and useful for the numerical and theoretical analysis.

4. Boundary Element and Augmented Lagrangian Methods for the Contact Problem

With the above results, we present our boundary element and augmented Lagrangian methods for the Coulomb frictional contact problem as follows.

We will detail the above algorithm in Section 5 using the semismooth Newton method.

Let and be the unique solution of the problems (1)–(8) and corresponding stress tensor on the boundary , respectively. From Lemma 1 and Lemma 2, and satisfy (15) and (16). In order to analyze the convergence of Algorithm 1, we introduce the following notation:and will use the following projection properties on [12, 13, 35].

Step 0: choose initial function and large enough , and set .
Step 1: find and such that with
Step 2: check some stopping criterion, and if it is satisfied then stop, else update and go to Step 1.

Lemma 3. For all generated by the first condition in Step 1 in Algorithm 1, it holds that

Lemma 4. For all generated by the second condition in Step 1 in Algorithm 1, it holds that

Proof. Let us separate into four subparts , , , and , whereFrom (16) and the second condition in Step 1 in Algorithm 1, we haveSince and have the same orientation, it follows that

Theorem 1. Let be the sequence generated by Algorithm 1; then, converges to in as .

Proof. Let , , and , then , , and . We take , , and in (33). This yieldsLet in the equation in Step 1 in Algorithm 1, then we haveSubtracting (40) from (41) results inUsing Lemmas 3 and 4 and Young’s inequality, we obtainFrom (31) and (42), it follows thatCombining (44) and (45) yieldsSumming up (46) with respect to , we obtainthis shows thatConsequently, converges to in .
Note that in (46), we can observe that the larger value of parameters and results in faster convergence to the algorithm. Therefore, this method is reliable and efficient for friction contact problems.

5. Numerical Examples

In this section, we present two test examples for the algorithm proposed in Section 4 for the solution of friction contact problems. In order to simplify the numerical process, we use constant boundary elements to approximate and and solve the corresponding linear systems. Here, we choose as the stopping criterion, and all computations are performed in Matlab codes. To deal with the nonlinearity of Algorithm 1, we apply the semismooth Newton method to Step 1 [12, 13] as given in Algorithm 2.

Step 0: choose and large enough , and set .
Step 1: determine
Step 2: solve with the boundary conditions and obtain and on .
Step 3: update and go to Step 1.
5.1. First Example

We consider an elastic square body with a rigid foundation and . The Neumann boundary condition is given by on , on , and on the rest of . And, the symmetry conditions (i.e., ) are applied on . Young’s modulus and Poisson’s ratio are  = 10000 and , respectively. This problem is not -elliptic but only -elliptic [38].

Let denote the number of boundary elements on . First, we apply our method to this problem with  = 320,  = 10000, and  = 10000, and the initial and deformed configurations of the body are given in Figure 1. In addition, the surface displacement and traction are depicted in Figures 2 and 3, respectively. We observe that is divided into two parts: the right part where the body remains in contact with the axis y = 0 and the left part of where it separates from this axis with a separation point near (0.26, 0), and the contact part is divided into a slip part and a stick part with a transition point from the slip to the stick near (0.47, 0). It can be seen that our results are in a good agreement with those in ref. [38].

Next, we solve the problem by choosing different parameters and to investigate the convergence behavior of our method. Table 1 gives the number of iterations for  = 320. We note that numerical results converge quickly when the parameters and are sufficiently large.

5.2. Second Example

In this example, we consider another model problem of the elastic square body in contact with a rigid foundation on the bottom side. In addition to the Neumann data applied on the top of the body, we set symmetric Neumann data on the left side and on the right side, respectively. Similarly, the symmetry conditions (i.e., ) are applied on for this example [38]. We use Young’s modulus  = 10000 and Poisson’s ratio for the test.

We choose  = 320,  = 10000, and  = 10000 again and apply our method to this problem. Figure 4 depicts the initial and deformed configurations of the half of the body (because of the symmetry). Figures 5 and 6 give the numerical results of the surface displacement and traction. One can observe that there is a transition point between contact and separation near (0.08, 0) on , which matches well with the numerical results in ref. [38]. In Table 2, we also list the number of iterations for the convergence with respect to different parameters and . As it can be seen from our tests, our method also converges quickly when and are sufficiently large.

6. Conclusion

In this paper, we have studied a new method for the solution of contact and friction problems and its convergence analysis. The method is developed by the ALM and the BEM. In particular, we have shown that our method shows unconditional convergence and converges quickly when the parameters and are sufficiently large. Although this method requires solving a nonlinear elasticity problem at each iteration, it can be easily dealt with by using the SSN. The numerical examples demonstrate the remarkable efficiency and reliability of the method.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was funded by the National Natural Science Foundation Project of CQ CSTC of China (Grant no. cstc2017jcyjAX0316) and the School Foundation Project of Chongqing University of Science and Technology of China (Grant no. CK2016Z07).