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

Authors: Chunxiao Yu; Dinghui Jing; Chang Fu; Yanfang Yang
Journal: Journal of Mathematics
 (2021)

---

## Abstract

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 influences of different 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 efficient and practical but also has much superiority. It is easy to implement, and it is fast convergent with good stability.


---

## Body

## 1. 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]. The FEM is relatively mature and widely used [7–10]. However, the BEM has the advantages of dimension reduction, singularity adaptation, high precision, and so on [11–14].The 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–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–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. Therefore, it is very urgent to develop some stable and efficient numerical algorithms [24–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–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–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. The 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. The research work will involve some mathematical, mechanical, and physical problems that are closely related to the frictional contact. The purpose is to provide new ideas and numerical methods for the solution of elastic frictional contact problems.This paper is organized as follows. In Section1, 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. Then, 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.


## 2. Basic Idea of the Penalty Function Method

For the optimization problem(1)minfxs.t.hix=0,i=1,2,…,l,we introduce a parameter λ and define an augmented objective function as follows:(2)Fx,λ=fx+λ∑i=1lhix2,where Fx,λ is called a penalty function and the parameter λ is called a penalty factor that is a very large positive number. When hix=0,i=1,2,…l, the penalty function Fx,λ is just equal to the objective function fx in equation (1); otherwise, its value will be very large and equation (1) will be transformed into the following unconstrained problem:(3)minFx,λ=minfx+λ∑i=1lhix2.


## 3. 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]:(4)cijujx+∫ΓTijx,yujydΓ=∫ΓUijx,ytjydΓ,where x indicates a source node, y indicates an arbitrary node in boundary Γ, cij indicates a boundary shape coefficient, and Uijx,y and Tijx,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]:(5)cijujxq+Rijmxq∑k,l,s1xq−ycujklϕlξsnmyξsJyξsωs+Simxq∑k,l,s1xq−ycujklϕlξsnmyξsyjξsJyξsωs−Pijxq∑k,l,s1xq−yctjklϕlξsJyξsωs−Qixq∑k,l,s1xq−yctjklϕlξsyjξsJyξsωs=0,where xq indicates a source node, yc indicates a multiple central node, s indicates an element integral node, ωs indicates the integral weight function of ξs, and Jyξs indicates a Jacobian determinant.With the given boundary conditions, equation (5) can be transformed into the following system of equation:(6)Ax=b,where x indicates an unknown column vector for displacements and surface force.


### 3.2. 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 Figure1.Figure 1
A sketch of nonpenetrating constraints.Namely, the following expression is satisfied:(7)ΔuA⋅n≤τ,where ΔuA 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.


## 3.1. Fundamental Formulas

For 3D elastic frictional contact problems, the boundary integral equation without consideration of body force is expressed as follows [6]:(4)cijujx+∫ΓTijx,yujydΓ=∫ΓUijx,ytjydΓ,where x indicates a source node, y indicates an arbitrary node in boundary Γ, cij indicates a boundary shape coefficient, and Uijx,y and Tijx,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]:(5)cijujxq+Rijmxq∑k,l,s1xq−ycujklϕlξsnmyξsJyξsωs+Simxq∑k,l,s1xq−ycujklϕlξsnmyξsyjξsJyξsωs−Pijxq∑k,l,s1xq−yctjklϕlξsJyξsωs−Qixq∑k,l,s1xq−yctjklϕlξsyjξsJyξsωs=0,where xq indicates a source node, yc indicates a multiple central node, s indicates an element integral node, ωs indicates the integral weight function of ξs, and Jyξs indicates a Jacobian determinant.With the given boundary conditions, equation (5) can be transformed into the following system of equation:(6)Ax=b,where x indicates an unknown column vector for displacements and surface force.


## 3.2. 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 Figure1.Figure 1
A sketch of nonpenetrating constraints.Namely, the following expression is satisfied:(7)ΔuA⋅n≤τ,where ΔuA 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.


## 4. Modeling and Optimization for the Node_Face Frictional Contact System Using the FM-BEM Penalty Function Method

### 4.1. Analysis of Node_Face Frictional Contact

We consider two objectsA and B in contact with each other. We suppose that object A (with fixed displacement constraints) is passive and object B is active. The numbers of discrete nodes are represented as NA and NB, respectively. Also, the numbers of contact nodes are represented as NAc and NBc, respectively. For the traditional BEM, the DOF of the final system of equations is 3NA+NB 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 objectB, it contacts with some element of object A, and its displacement can be obtained by the interpolation of its contact element nodes’ displacements. Then, 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. The node_face frictional contact constraints are shown in the following expressions:Stick state:(8)UkB=∑l=1Mφlξ1,ξ2UklA,k=1,2,3.Slip state:(9)T1B=−Ttcosθ,T2B=−Ttsinθ,U3B=∑l=1Mφlξ1,ξ2U3lA.For each contact node of objectA, it contacts with some element of object B, and its surface force can be obtained by the interpolation of its contact element nodes’ force. Then, 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. The node_face frictional contact constraints are shown in the following expressions:Stick state:(10)TkA=−∑l=1Mφlξ1,ξ2TklB,k=1,2,3.Slip state:(11)T1A=Ttcosθ,T2A=Ttsinθ,T3A=−∑l=1Mφlξ1,ξ2T3lB.In equations (8)–(11), UkB indicates k-direction displacement of each node in object B, UklA indicates k-direction displacement of node l in object A, TkA 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, Tt 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. Then, the total DOF of the contact system becomes 3NA+NB+NAc+NBc. For convenience, it can be written as NF.


### 4.2. 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. Then, mathematical programming is conducted on the frictional contact system and an optimization mathematical model is established. The detailed process is as follows:For the objectB, the system of equations formed by the traditional BEM can be expressed as A′x=b′, where A′=aij′,b′=bi′. Let(12)fiBxj=bi′−∑j=1NFaij′xj,i=1,2,…,3NB;j=1,2,…,NF.After equations (8) and (9) have been linearized, according to equation (12), they can be rewritten as(13)fikBxj=UikB−∑l=1Mφlξ1,ξ2UiklA,i=1,2,…,3NBc;j=1,2,…,NF;k=1,2,3,(14)fi1Bxj=Ti1B+Titcosθ,fi2Bxj=Ti2B+Titsinθ,fi3Bxj=Ui3B−∑l=1Mφlξ1,ξ2Ui3lA,i=1,2,…,3NBc;j=1,2,…,NF.For the objectA, the system of equations formed by the traditional BEM can be expressed as A〞x=b〞, where A〞=aij〞,b〞=bi〞. Let(15)fiAxj=bi〞−∑j=1NFaij〞xj,i=1,2,…,3NA;j=1,2,…,NF.After equations (10) and (11) have been linearized, according to equation (15), they can be rewritten as(16)fikAxj=TikA+∑l=1Mφlξ1,ξ2TiklB,i=1,2,…,3NAc;j=1,2,…,NF;k=1,2,3,(17)fi1Axj=Ti1A−Titcosθ,fi2Axj=Ti2A−Titsinθ,fi3Axj=Ti3A+∑l=1Mφlξ1,ξ2Ti3lB,i=1,2,…,3NAc;j=1,2,…,NF.According to equations (12)–(17), let x=xjj=1,2,…,NF, and we define an objective function for the nonlinear analysis of node_face frictional contact as follows:(18)fx=∑i=13NBfiBx2+∑i=13NBc∑k=13fikBx2+∑i=13NAfiAx2+∑i=13NAc∑k=13fikAx2.According to equations (12)–(18), an optimization mathematical model for the node_face frictional contact system can be established as follows:(19)minfx,fiBx=0,i=1,2,…,3NB,fiBx=0,i=1,2,…,3NBc,fiAx=0,i=1,2,…,3NA,fiAx=0,i=1,2,…,3NAc.


### 4.3. 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. The 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 E0). Otherwise, the spring will prevent the objects from contacting and do work, so the potential energy (written as Ec) will sharply increase. For the node_face frictional contact system, let(20)fx=b−Ax,x=xj,j=1,2,…,NF.We construct an energy objective function as follows:(21)E=E0+Ec,where(22)E0=fx22=∑i=13NBfiBx2+∑i=13NBcfiBx2+∑i=13NAfiAx2+∑i=13NAcfiAx2,Ec=12α∑i=1NFδidi−π2,δi=1,di−τ>0,0,di−τ≤0,where di 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:(23)minE=minE0+Ec.So, the solution of node_face frictional contact is transformed into an unconstrained optimization problem.


### 4.4. Selection of Penalty Factor

From equations (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.Figure 2
Influence of a penalty factorα on objective function values.According to Figure2, when the penalty factor α varies within the range between 10 and 108, 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 108. While penalty factor α is larger than 108, the objective function value will increase sharply, namely, equation (7) cannot be satisfied, and the “spring” will carry out punishments on the contact. Then, the system cannot be solved properly.


## 4.1. Analysis of Node_Face Frictional Contact

We consider two objectsA and B in contact with each other. We suppose that object A (with fixed displacement constraints) is passive and object B is active. The numbers of discrete nodes are represented as NA and NB, respectively. Also, the numbers of contact nodes are represented as NAc and NBc, respectively. For the traditional BEM, the DOF of the final system of equations is 3NA+NB 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 objectB, it contacts with some element of object A, and its displacement can be obtained by the interpolation of its contact element nodes’ displacements. Then, 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. The node_face frictional contact constraints are shown in the following expressions:Stick state:(8)UkB=∑l=1Mφlξ1,ξ2UklA,k=1,2,3.Slip state:(9)T1B=−Ttcosθ,T2B=−Ttsinθ,U3B=∑l=1Mφlξ1,ξ2U3lA.For each contact node of objectA, it contacts with some element of object B, and its surface force can be obtained by the interpolation of its contact element nodes’ force. Then, 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. The node_face frictional contact constraints are shown in the following expressions:Stick state:(10)TkA=−∑l=1Mφlξ1,ξ2TklB,k=1,2,3.Slip state:(11)T1A=Ttcosθ,T2A=Ttsinθ,T3A=−∑l=1Mφlξ1,ξ2T3lB.In equations (8)–(11), UkB indicates k-direction displacement of each node in object B, UklA indicates k-direction displacement of node l in object A, TkA 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, Tt 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. Then, the total DOF of the contact system becomes 3NA+NB+NAc+NBc. For convenience, it can be written as NF.


## 4.2. 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. Then, mathematical programming is conducted on the frictional contact system and an optimization mathematical model is established. The detailed process is as follows:For the objectB, the system of equations formed by the traditional BEM can be expressed as A′x=b′, where A′=aij′,b′=bi′. Let(12)fiBxj=bi′−∑j=1NFaij′xj,i=1,2,…,3NB;j=1,2,…,NF.After equations (8) and (9) have been linearized, according to equation (12), they can be rewritten as(13)fikBxj=UikB−∑l=1Mφlξ1,ξ2UiklA,i=1,2,…,3NBc;j=1,2,…,NF;k=1,2,3,(14)fi1Bxj=Ti1B+Titcosθ,fi2Bxj=Ti2B+Titsinθ,fi3Bxj=Ui3B−∑l=1Mφlξ1,ξ2Ui3lA,i=1,2,…,3NBc;j=1,2,…,NF.For the objectA, the system of equations formed by the traditional BEM can be expressed as A〞x=b〞, where A〞=aij〞,b〞=bi〞. Let(15)fiAxj=bi〞−∑j=1NFaij〞xj,i=1,2,…,3NA;j=1,2,…,NF.After equations (10) and (11) have been linearized, according to equation (15), they can be rewritten as(16)fikAxj=TikA+∑l=1Mφlξ1,ξ2TiklB,i=1,2,…,3NAc;j=1,2,…,NF;k=1,2,3,(17)fi1Axj=Ti1A−Titcosθ,fi2Axj=Ti2A−Titsinθ,fi3Axj=Ti3A+∑l=1Mφlξ1,ξ2Ti3lB,i=1,2,…,3NAc;j=1,2,…,NF.According to equations (12)–(17), let x=xjj=1,2,…,NF, and we define an objective function for the nonlinear analysis of node_face frictional contact as follows:(18)fx=∑i=13NBfiBx2+∑i=13NBc∑k=13fikBx2+∑i=13NAfiAx2+∑i=13NAc∑k=13fikAx2.According to equations (12)–(18), an optimization mathematical model for the node_face frictional contact system can be established as follows:(19)minfx,fiBx=0,i=1,2,…,3NB,fiBx=0,i=1,2,…,3NBc,fiAx=0,i=1,2,…,3NA,fiAx=0,i=1,2,…,3NAc.


## 4.3. 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. The 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 E0). Otherwise, the spring will prevent the objects from contacting and do work, so the potential energy (written as Ec) will sharply increase. For the node_face frictional contact system, let(20)fx=b−Ax,x=xj,j=1,2,…,NF.We construct an energy objective function as follows:(21)E=E0+Ec,where(22)E0=fx22=∑i=13NBfiBx2+∑i=13NBcfiBx2+∑i=13NAfiAx2+∑i=13NAcfiAx2,Ec=12α∑i=1NFδidi−π2,δi=1,di−τ>0,0,di−τ≤0,where di 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:(23)minE=minE0+Ec.So, the solution of node_face frictional contact is transformed into an unconstrained optimization problem.


## 4.4. Selection of Penalty Factor

From equations (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.Figure 2
Influence of a penalty factorα on objective function values.According to Figure2, when the penalty factor α varies within the range between 10 and 108, 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 108. While penalty factor α is larger than 108, the objective function value will increase sharply, namely, equation (7) cannot be satisfied, and the “spring” will carry out punishments on the contact. Then, the system cannot be solved properly.


## 5. Penalty Optimization IGMRES (m) Algorithm

To solve equation (23), a penalty optimization IGMRES (m) algorithm is presented. The detailed process is as follows:(1)
Initialization: for a fixed parameterm, we set an appropriate precision ε and a parameter q (2≤q≤m). We take an initial value x0 and compute(24)r0=b−Ax0,β=r0,v1=r0β,V1=v1,E0=Ex0.(2)
Iteration: forj=1,2,…,m, we have(1)
Incomplete orthogonalization:(25)hij=Avj,vii=j0,…,j,v^j+1=Avj−∑i=j0jhijvi,j0=max1,j−q−1.(2)
Standardization:(26)hj+1,j=v^j+1,vj+1=v^j+1hj+1,j.(3)
Updation ofVj+1 and H¯j:(27)Vj+1=Vj,vj+1,H¯j=H¯j−1hij0hj+1,jj+1×j,whereH¯j indicates an upper Hessenberg matrix. Then, we have(28)AVm=Vm+1H¯m.Whenj=1, the first column will be omitted.(3)
We solve the following least square problem to obtainym:(29)rm=minym∈Cmβe1−H¯mym.(4)
We construct the approximate solution:(30)xm=x0+Vmym.(5)
The modules of residual vectors and the value of energy objective function are computed.(31)rm=f−Axm,Em=Exm.(6)
Restart judgment: ifrm≤ε and Em<E0, then let x=xm and stop. Otherwise, reset x0=xm and return to the initialization step.


## 6. Numerical Example

We consider two elastic objectsA 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. The two objects are isotropic with Young’s modulus E = 210 Gpa, Poisson’s ratio ν = 0.3, and a frictional coefficient f = 0.1. The object B is subjected to a uniform load. The total load is divided into six steps, and the contact distance tolerance is τ = 0.0001 mm. The computation model and discrete mesh are shown in Figure 3, and the discrete data are shown in Table 1.Figure 3
Computation model and discrete mesh.Table 1
Discrete data.
ObjectsTotalityContact areaDOFNode numberElement numberNode numberElement numberA165016485675206651B119811964594684971Total28482844102698811622When the objectB 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 108, and the obtained results agree well with the theoretical analysis. The 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. These 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.Figure 4
Displacement distributions in the contact area.Figure 5
Pressure distributions in the contact area.Figure 6
Friction distributions in the contact area.Figure 7
Friction slip field in the contact area.Figure 8
Comparison of contact displacements under different loads.Figure 9
Comparison of contact pressures under different loads.WhenP = 1 GPa and the penalty factor α is taken as 109, 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 1010, 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×1011 to 2.1×1011, the solution will be impossible. The 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.Figure 10
Friction slip field in the contact area. (a) Friction slip (P = 1 GPa). (b) Friction slip (P = 10 GPa).
(a)(b)Figure 11
Influence of a penalty factor on the computational efficiency.


## 7. 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. The influence of penalty factor on the solution process was analyzed, and a penalty optimization IGMRES (m) algorithm was presented. The 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. Theoretical analysis and numerical experiment showed that the new method had much superiority in efficiency, applicability, easy numerical implementation, fast convergence, stability, etc. The 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.



---

*Source: 6626647-2021-01-13.xml*