A Kind of FM-BEM Penalty Function Method for a 3D Elastic Frictional Contact Nonlinear System

In this paper, a kind of node_face frictional contact FM-BEM penalty function method is presented for 3D elastic frictional contact nonlinear problems. According to the principle of minimum potential energy, nonpenetrating constraints are introduced into the elastic frictional contact system as a penalty term. By using the least square method and penalty function method, an optimization mathematical model and a mathematical programming model with a penalty factor are established for the node_face frictional contact nonlinear system. For the two models, a penalty optimization IGMRES (m) algorithm is proposed, and the inﬂuences of diﬀerent penalty factors on the solution of the whole system are analyzed. Finally, a numerical simulation is carried out for two elastic frictional contact objects, and some important results including displacements, pressures, friction forces, and friction slips in the contact area are presented. Theoretical analysis and numerical experiment show that the newly presented FM-BEM penalty function method not only is eﬃcient and practical but also has much superiority. It is easy to implement, and it is fast convergent with good stability.


Introduction
Elastic frictional contact is a multiple nonlinear problem [1,2], and it is necessary to accurately track the motion of the object before contact and the interaction between objects after contact, which includes the correct simulation of friction and deformation behavior between contact surfaces and the analysis of the possible energy conversion problem. For the contact problems, only very few of them can be solved by analytical methods, and most of them need to be simulated by numerical methods such as the Finite Element Method (FEM) [3,4] and the Boundary Element Method (BEM) [5,6]. e FEM is relatively mature and widely used [7][8][9][10]. However, the BEM has the advantages of dimension reduction, singularity adaptation, high precision, and so on [11][12][13][14].
e penalty function method [15,16] is a common method to solve optimization problems, and it is also one of the effective methods to solve an elastic contact problem [17][18][19]. Without increasing the system's Degree of Freedom (DOF), this method can be used to directly apply constraints to the two contact objects. Many scholars have used it to solve the frictional contact problems in different fields [20][21][22][23]. In engineering, gradient-based optimization algorithms, for example, the existing FEM such as the Lagrange multiplier method and penalty function method, are often used to solve the contact problems. For the case of nonfrictional contact, sufficiently stable results can be obtained. For the case of frictional contact, severe numerical oscillation may occur with the change of loads or meshes, and it will be very difficult to obtain a stable result unless special treatments are made. In addition, the procedures of existing numerical algorithms are usually complicated and much memory space and computing time are required, so repeated checking and revision are needed to obtain suitable results. At present, various kinds of commercial computing software often fail to give accurate and reliable results for the analysis of frictional contact. erefore, it is very urgent to develop some stable and efficient numerical algorithms [24][25][26][27].
In recent years, the Fast Multipole Boundary Element Method (FM-BEM) [28,29] has attracted much attention as a kind of new and efficient numerical method [30][31][32][33][34]. Our research group studied the mathematical and mechanical theories of the FM-BEM from the perspective of fundamental solution. By using the superiorities of FM-BEM such as high precision, high computational efficiency and being suitable for large-scale computing, we have successfully applied it to the numerical analysis of elastic frictional contact problems and have completed some simulations [35][36][37][38], for example, the interference fit between taper sleeve and roll neck of an oil film bearing and a surface force field of screw pair in a rolling mill.
For the study of elastic frictional contact problems, the penalty function method in the existing literature was used to solve some optimization problems with a node_node contact mode. e BEM and FM-BEM focused on the modeling and numerical analysis for the nonpenetrating contact mode and often failed to give numerical results for the penetrating contact mode. According to the abovementioned analysis, we will present a kind of FM-BEM penalty function method to solve the elastic node_face frictional contact problems. As the same time, we will establish a mathematical programming model with a penalty factor and propose a penalty optimization algorithm. In this method, some important factors will be synthetically considered, which include the deformation and stress condition in a contact process, the nonlinearity of boundary condition for the contact surface, the size and mutual position of the contact area, the change of contact state, and so on. e research work will involve some mathematical, mechanical, and physical problems that are closely related to the frictional contact. e purpose is to provide new ideas and numerical methods for the solution of elastic frictional contact problems. is paper is organized as follows. In Section 1, basic thought of the Penalty Function Method is introduced. In Section 2, fundamental formulas and frictional contact condition for the 3D elastic frictional contact FM-BEM are presented. In Section 3, interpolation constraints are analyzed for the node_face frictional contact nonlinear system. en, an optimization mathematical model and a mathematical programming model with a penalty factor are established by using the least square method and penalty function method. In Section 4, a penalty optimization IGMRES (m) algorithm is proposed. In Section 5, a simulation of two elastic objects' frictional contact process is provided and numerical analysis is completed. At last, the concluding remarks are presented.

Basic Idea of the Penalty Function Method
For the optimization problem we introduce a parameter λ and define an augmented objective function as follows: where F(x, λ) is called a penalty function and the parameter λ is called a penalty factor that is a very large positive number. When h i (x) � 0, (i � 1, 2, . . . l), the penalty function F(x, λ) is just equal to the objective function f(x) in equation (1); otherwise, its value will be very large and equation (1) will be transformed into the following unconstrained problem:

3D Elastic Frictional Contact FM-BEM
3.1. Fundamental Formulas. For 3D elastic frictional contact problems, the boundary integral equation without consideration of body force is expressed as follows [6]: where x indicates a source node, y indicates an arbitrary node in boundary Γ, c ij indicates a boundary shape coefficient, and U ij (x, y) and T ij (x, y) indicate the kernel functions of displacement and surface force fundamental solutions, respectively. By the FM-BEM, equation (4) can be discretized as follows [29]: where x q indicates a source node, y c indicates a multiple central node, s indicates an element integral node, ω s indicates the integral weight function of ξ s , and J [y(ξ s )] indicates a Jacobian determinant. With the given boundary conditions, equation (5) can be transformed into the following system of equation: where x indicates an unknown column vector for displacements and surface force.

Frictional Contact Condition.
When two objects contact each other, in order to ensure the balance and stability, the contact system must be satisfied with nonpenetrating constraints, as is shown in Figure 1. Namely, the following expression is satisfied: where Δ u A indicates the displacement vector increment of the node A, n indicates the unit normal vector, and τ indicates the tolerance of contact distance. Otherwise, once penetration occurs in the contact area, the system solution will not be carried out normally.

Analysis of Node_Face Frictional Contact.
We consider two objects A and B in contact with each other. We suppose that object A (with fixed displacement constraints) is passive and object B is active. e numbers of discrete nodes are represented as N A and N B , respectively. Also, the numbers of contact nodes are represented as N c A and N c B , respectively. For the traditional BEM, the DOF of the final system of equations is 3 (N A + N B ) and the displacements and surface force for each contact node are unknown. As a result, for each contact node, three supplement equations must be established.
For each contact node of object B, it contacts with some element of object A, and its displacement can be obtained by the interpolation of its contact element nodes' displacements.
en, displacement constraints are established. According to Coulomb's Law of Friction, if relative slip occurs between the contact node and its contact surface, tangential displacement constraints can be replaced by tangential friction ones. e node_face frictional contact constraints are shown in the following expressions: Stick state: Slip state: For each contact node of object A, it contacts with some element of object B, and its surface force can be obtained by the interpolation of its contact element nodes' force. en, surface force constraints are established. Similarly, if relative slip occurs between the contact node and its contact surface, tangential surface force constraints can be replaced by tangential friction ones. e node_face frictional contact constraints are shown in the following expressions: Stick state: Slip state: In equations (8)- (11), U B k indicates k-direction displacement of each node in object B, U A k l indicates k-direction displacement of node l in object A, T A k indicates k-direction surface force of each node in object A, M indicates the node number of a contact element, φ l indicates the interpolation function, T t indicates the friction at t moment, (ξ 1 , ξ 2 ) indicates the local coordinate, and θ indicates a slip angle. Here, ξ 1 , ξ 2 , θ can be predetermined by the least square method. According to equations (8)- (11), three supplement equations can be established for each contact node. en, the total DOF of the contact system becomes . For convenience, it can be written as NF.

Optimization Mathematical Model for Node_Face Frictional Contact.
Node_face frictional contact constraints show high nonlinearity, which results in a very difficult and time-consuming solution procedure. To accelerate the iterative convergence, nonlinear contact constraints will be linearized. At first, the least square method is applied to equations (8)-(11) to obtain ξ 1 , ξ 2 , θ while the contact behavior is precisely simulated. en, mathematical programming is conducted on the frictional contact system and an optimization mathematical model is established. e detailed process is as follows: For the object B, the system of equations formed by the traditional BEM can be expressed as After equations (8) and (9) have been linearized, according to equation (12), they can be rewritten as For the object A, the system of equations formed by the traditional BEM can be expressed as After equations (10) and (11) have been linearized, according to equation (15), they can be rewritten as According to equations (12)- (17), let x � (x j ) (j � 1, 2, . . . , NF), and we define an objective function for the nonlinear analysis of node_face frictional contact as follows: ‖f According to equations (12)- (18), an optimization mathematical model for the node_face frictional contact system can be established as follows:

Penalty Factor Programming Model for Node_Face
Frictional Contact. From the abovementioned analysis, when the contact system is stable, the involved objects satisfy nonpenetrating constraints shown in equation (7); otherwise, the penalty function method can be used to apply contact constraints. We suppose that there is a "spring" between a possible contact node and its contact surface and its compressive stiffness is very large while the tensile stiffness is zero. e stiffness is taken as a penalty factor and written as α. According to the principle of minimum potential energy, when two objects contact each other, if equation (7) is satisfied, the work performed by the spring will be zero, namely, the penalty factor α � 0, and the contact system will be stable with minimum potential energy (written as E 0 ). Otherwise, the spring will prevent the objects from contacting and do work, so the potential energy (written as E c ) will sharply increase. For the node_face frictional contact system, let We construct an energy objective function as follows: where where d i indicates the distance between a contact node and its contact surface. For the node_face frictional contact system, we suppose the contact surface is smooth and the deformation is very small. According to equations (3) and (19)-(21), a penalty factor programming model can be established as follows: So, the solution of node_face frictional contact is transformed into an unconstrained optimization problem. (3) and (23), we know the optimization of penalty factor α is very important. For each factor α, a corresponding objective function value can be obtained, and it will increase with the increase of α. When α � 0, equation (23) has the same solution as equation (6). While α ⟶ ∞, the solution of equation (23) will converge to the analytical solution, and abnormalities may occur for a too large factor α. So, the penalty factor α should not be taken as a too large value. From energetics point of view, the penalty factor α is equivalent to spring stiffness. When an object is subjected to fixed loads, the factor α will be inversely proportional to the deformation increments within the elastic range. For example, if two elastic objects contact each other, the relationship of the factor α and the objective function value can be as shown in Figure 2.

Selection of Penalty Factor. From equations
According to Figure 2, when the penalty factor α varies within the range between 10 and 10 8 , the objective function value will be close to zero, namely, the nonpenetrating constraints expressed in equation (7) can be satisfied and the system will be stable. So, the penalty factor α can be taken as 10 8 . While penalty factor α is larger than 10 8 , the objective function value will increase sharply, namely, equation (7) cannot be satisfied, and the "spring" will carry out punishments on the contact. en, the system cannot be solved properly.
Journal of Mathematics 5

Penalty Optimization IGMRES (m) Algorithm
To solve equation (23), a penalty optimization IGMRES (m) algorithm is presented. e detailed process is as follows: (1) Initialization: for a fixed parameter m, we set an appropriate precision ε and a parameter q (2 ≤ q ≤ m). We take an initial value x (0) and compute (2) Iteration: for j � 1, 2, . . . , m, we have (1) Incomplete orthogonalization: (2) Standardization: (3) Updation of V j+1 and H j : where H j indicates an upper Hessenberg matrix. en, we have When j � 1, the first column will be omitted. (3) We solve the following least square problem to obtain y m : (4) We construct the approximate solution: (5) e modules of residual vectors and the value of energy objective function are computed.

Numerical Example
We consider two elastic objects A and B in contact with each other. A is a support with a width W � 50 mm, a height H � 30 mm, and a length L1 � 50 mm. B is a half cylinder with a radius R � 15 mm and a length L2 � 60 mm. e two objects are isotropic with Young's modulus E � 210 Gpa, Poisson's ratio ] � 0.3, and a frictional coefficient f � 0.1. e object B is subjected to a uniform load. e total load is divided into six steps, and the contact distance tolerance is τ � 0.0001 mm. e computation model and discrete mesh are shown in Figure 3, and the discrete data are shown in Table 1.
When the object B is subjected to a uniform load P that is not more than 1 GPa, the penalty factor is taken as a value that ranges from 10 to 10 8 , and the obtained results agree well with the theoretical analysis. e solution process is very stable. When P � 100 MPa, the distributions of contact displacement, pressure and friction force, and the friction slip field are presented, as is shown in Figures 4-7. ese results agree well with the experimental analysis. In addition, the contact displacements and pressures under different loads are compared, as is shown in Figures 8, and 9, which shows that the edge effect is becoming more and more obvious with the load increase.
When P � 1 GPa and the penalty factor α is taken as 10 9 , the solution procedure is abnormal and the friction directions of some contact nodes change, as is shown in Figure 10(a). When P � 10 GPa and the penalty factor α is taken as 10 10 , the solution procedure is more abnormal, as is shown in Figure 10(b). When P ≥ 100 GPa and the penalty factor α is taken as a value ranging from 1.0 × 10 11 to 2.1 × 10 11 , the solution will be impossible. e reason why the solution is impossible or abnormal is that penetration occurs when two objects contact each other, namely, equation (7) is not satisfied.          In addition, numerical experiments show that when the loads increase gradually, if the penalty factor α is taken as an inappropriate value, the computation time will sharply increase, as is shown in Figure 11.

Conclusions
Based on the FM-BEM, a node_face frictional contact mode was analyzed and nonpenetrating constraints were presented for 3D elastic frictional contact problems. For the case of frictional contact without penetration, nonlinear contact constraints were linearized by use of the least square method and an optimization mathematical model with node_face frictional contact mode was established. For the case of frictional contact with penetration, according to the principle of minimum potential energy, a penalty function method was used to apply the contact constraints and an energy objective function was constructed; then, a node_face frictional contact analysis was transformed into an unconstrained optimization problem. For the elastic frictional contact FM-BEM problem, nonpenetrating constraints were introduced into the system as a penalty term. Without increasing the system variables, a penalty factor mathematical programming model was established by the penalty function method. e influence of penalty factor on the solution process was analyzed, and a penalty optimization IGMRES (m) algorithm was presented. e frictional contact of two elastic objects under different loads was numerically simulated, and the results of displacements, pressures, friction forces and friction slips in the contact area were obtained. eoretical analysis and numerical experiment showed that the new method had much superiority in efficiency, applicability, easy numerical implementation, fast convergence, stability, etc. e proposed FM-BEM penalty functional method could provide new ideas and methods for the solution of frictional contact problems and related mathematical, mechanical, and physical problems.

Data Availability
e data used to support the results of this study are included within the article.