Abstract
Frictionless contact problems are the simplest and classical contact problems, and the contact conditions of sticking, slipping, and separation mode all can be ascribed to complementary problems. Consequently, a smooth system of equations approach for the design and analysis of complementarity problems for frictionless contacts is presented. A compute program based on boundary element technique is given and applied to two practical contact examples. The validity and accuracy of the proposed method are demonstrated.
1. Introduction
Contact problems are of particular importance in many engineering applications [1, 2] such as gears, piles, retaining walls, and tunnel lining. Since the hertz theory was developed in 1881, much research has been developed in this area including both theoretical and experimental work [3–7]. Hence it establishes the foundation for modern contact mechanics. In the analysis of contact problems, special attention must be required because the actual contact area between the contacting bodies is usually not known in advance, and the character of interface between contact bodies largely determines the deformation, movement, and stress distribution.
There is a large literature on numerical methods for contact problems. Roughly speaking, numerical algorithms can be classified into three categories. The first class is known as direct iterative algorithm [8–11], which assumes the contact region and contact status firstly and then solves the problem and checks whether the assumption is correct or not. It solves the problem iteratively by trial-and-error, thus requiring much computational effort; the second class is contact constraint algorithm [12–14], which deals with the contact boundary properly, and transforms the constrained optimization problem into nonconstrained optimization problem. According to the different methods for unconstrained optimization, the penalty method [12], Lagrange method [13], and augmented Lagrange method [14] have been introduced. The calculation accuracy based on penalty method strongly depends on the penalty coefficient, and the coefficient is decided by experience; the number of unknown variables, computation time, and memory requirement using Lagrange method are troublesome; to foster strengths and circumvent weaknesses, the augmented Lagrange method is proposed for contact problems, which decreases ill-conditioning of governing equations, and satisfies exactly constraints with finite penalties [15]. However, these methods can be generalized into the iterative algorithms. The third class is known as the mathematical programming method [16, 17], and the solutions are obtained by using either linear programming or quadratic programming techniques [18–20]. The advantage is that the original problem can be converted to programming problem by expressing the normal and tangential conditions into the complementary expression, and mature mathematical methods can be applied. Furthermore, it is known that frictionless/frictional contact problems can be formulated as complementarity problems. Hence, a numerical method for complementarity problems can be applicable to contact problems [21]. Here, the system of nonlinear equations method [22] has been widely used, which transforms the complementarity problem into a system of nonlinear equations by employing the nonlinear complementary functions.
For numerical discretization, the boundary element method (BEM) [23] is particularly well suited to solve contact problems [24]. Except for the reduced dimensionality by one, the most striking feature of the BEM is that the tractions and displacements are computed to the same degree of accuracy, which is an important feature if reliable solutions are to be obtained [25]. The BEM was first applied to contact problems by Andersson et al. [26] in two-dimensional frictionless problems and later extended the applications to frictional problems [27]. Afterwards, other contact problems have been studied such as elastoplastic contact problems [28] and 3D frictional problems [29]. However, these problems had been solved by the iterative algorithms, and the trial-and-error, calculation accuracy and computation time, and so forth had frequently appeared. So far, the system of nonlinear equations method has not been used to solve contact problems by BEM.
In this paper, the smooth system of equations is employed to solve the two-dimensional elastic frictionless contact problems. It is comprised of the nonlinear complementary functions describing the relationship between the gap and contact pressure for any node pair and boundary integral equations. The presented approach does not need repeatedly trial calculation to decide the contact mode for any node pair. Meanwhile, the proposed method is easy to be accepted and used. According to the results of given load, the contact state can be observed obviously, and it does not need your judgment. This algorithm is implemented in a 2D BEM code and verified using two numerical examples. The results by the proposed algorithm match well with the analytical solutions and the FEM results and clearly demonstrate the feasibility and flexibility of the proposed approach for 2D contact analysis.
2. Contact Problems and Its Complementarity Description
2.1. General Description of Contact Problems
In this section, two contact bodies and are considered as shown in Figure 1. The boundary of any body is composed of three disjoint parts: displacement boundary , traction boundary , and potential contact region . Nevertheless, the region has been taken sufficiently large to contain actual contact regions.
Considering a pair of points and on the bodies and , respectively, on the contact boundary , the following contact modes shown in Figure 2 may happen: separation, stick and slip mode. For the convenient description of contact conditions, it is necessary to define a local coordinate system for the pair of points and . If body is selected as the target body, the coordinate system should be established in point as show in Figure 3, and and denote unit normal and tangential directions, respectively. Selecting (the subscript refers to direction of force, and superscript refers to the contact point) and (the subscript refers to the pair of points) to represent the traction and gap between the point pair, respectively, the contact conditions are listed as follows.
(a) Separation mode
(b) Stick mode
(c) Slip mode
(1) Separation Mode. The individual traction at the pair of points and is zero, and the gap could be positive:
To maintain the consistency below, the relationship can be recast to .
(2) Stick Mode. The individual tangential traction at the pair of points and is zero, and the total normal tractions are equal to zero. The gap should be zero:
(3) Slip Mode. Because frictionless contact is considered, the individual tangential traction at the pair of points and is still zero, and the sum of normal traction is equal to zero. The gap should be zero:
It should be mentioned that the individual normal contact traction at the pair of points and is always less than zero (i.e., compressive). Here, subjected to the local coordinate system of point on the body , the normal contact pressure can be expressed asFurthermore, can be expressed aswhere defines the initial gap between the points and , and and define the displacement vector on the bodies of and , respectively. denotes the unit vector at the contacting point of body .
2.2. Complementarity Theory for Contact Problems
2.2.1. Complementarity Problem
The complementary problem [30] is an important optimization problem. It is widely employed in many problems, such as game theory, economy analysis, and traffic equilibrium. It is firstly proposed by Dantzig and Cottle in 1963. It can be stated as follows [31–33].
Definition 1. For a given map from the -dimensional Euclidean space into itself, find satisfyingWhen is the form , the above problem is referred to as the linear complementarity problem (LCP); otherwise, it is called the nonlinear complementarity problem (NCP).
In the past few decades, the complementary problem has attracted much attention because of its wide applications. Consequently, the algorithm has achieved fruitful results and mainly includes the Lemke algorithm [34], homotopy method [35], projection algorithm [36], interior point algorithm [37], and system of equations algorithm [22]. One of the most powerful and popular methods is to reformulate the complementary problem as the system of equations. To construct it, a class of functions, called NCP-functions, defined below, plays an important role.
Definition 2. A function is called an NCP-function if it satisfiesCommon NCP functions are as follows:(1),(2),(3).
The most wonderful feature of NCP function is that it transforms the problem containing two equalities and an inequality into a problem only containing an equality. Therefore, the complex contact problem can be solved by solution of the system of equations.
2.2.2. NCP-Function for Contact Problems
Using and to stand for the gap and contact pressure of any potential contacting pair of points, the normal contact conditions can be described as
Furthermore, the gap and contact pressure can be described by the following graph, as shown in Figure 4.
According to (8), the following relationships can be achieved:
Obviously, it is a complementary description between the gap and minus normal contact stress . Consequently, the above description can be expressed by the NCP function mentioned above. Using the function , we can express (9) as
Figure 5 shows the graph of the function above and it can be observed that the proposed function is smooth everywhere and well suitable for contact problems.
3. Complementarity Problem Formulation by BEM
The BEM formulation for an elastic continuum is well known and can be found in many classical texts such as Brebbia et al. [38] and Aliabadi [39]. The elastostatic boundary integral equation for a boundary point with no body force is given as follows:where is the source point and is the field point at the boundary. is the free coefficient of geometry. and represent the fundamental solutions for displacement and traction components, respectively. The boundary consists of displacement boundary , traction boundary , and potential contact region .
In order to perform numerical analysis, the boundary is discretized into linear elements. Equation (11) can be conveniently expressed in the following matrix form:
After each domain is treated separately to form the matrix above, the resulting matrices and are coupled together according to the relevant contact conditions. The total matrices for two contact bodies can be written as follows:where the superscripts and refer to the two bodies in contact. Note that under this arrangement the matrices are not fully populated. After the application of boundary conditions, (13) can be recast aswhere denotes the boundary unknowns and is the contribution of known boundary variables, that is, values prescribed by the boundary conditions. It is noted that the matrix is not a square matrix, so the contact conditions will be incorporated.
For any node pair, there are four unknowns for each node, namely, normal and tangential displacements and tractions, , and they are referred to a local coordinate system. Except for the two boundary integral equations for every contact node, other four complementary equations will be listed to solve the system of equations.
Take the pair of nodes and as an example; the following relationships hold:
So far, the numbers of unknowns are equal to the number of equations, and the system of equations can be solved by Newton method [40].
Consequently, the calculation process has been exhibited as shown in Figure 6, and the corresponding compute program by BEM is formulated. It should be noted that the proposed algorithm does not need to judge the contact modes in the calculating process. After solving the system of equations, the displacements and tractions for any potential node pair are shown, and the contact mode is clear at a glance.
4. Numerical Examples
It should be mentioned that contact problem widely exists in the static and dynamic problem. The difference between static and dynamic problem is that whether the effect of acceleration can be ignored or not. Although this paper mainly considers the static contact problem, the proposed method can be easily extended to dynamic problem.
A computer program by BEM is applied to two practical contact problems. The BEM results are compared to analytical solutions and FEM results to establish their accuracy.
Example 1 (cylinder and base contact problem). In this problem, a cylinder with a radius of 5 m is pressed against a 5 m deep base under plane strain conditions shown in Figure 7(a). Due to the symmetry of problem, only the quarter of cylinder and the half of base have been modeled as shown in Figure 7(b). Here, the uniform load along the cylinder’s top per unit axial length ( KN/M) can be expressed by the load divided by the diameter of cylinder ( KN/M2). Both the cylinder and base are given the same properties: elastic modulus MPa and Poisson’s ratio .
(a) Geometry
(b) Simplified model
Such problems are referred to as Hertz-type or Hertzian contact problems, and the analytical solution can be given as follows [41]:where is the half-width of the contact, is the maximum pressure at the center of the contact, and is the pressure distribution of the contact. and are given by
In the boundary element analysis, the cylinder and base consist of 31 and 37 linear elements as shown in Figure 8(a), respectively. In order to obtain more accurate stress, a very fine mesh of the potential contact region is designed in Figure 8(b). The boundary conditions are as follows: the vertical displacements are fixed at the bottom of the base, and the horizontal displacements are fixed for the left boundaries of the cylinder and base.
(a) Boundary discretization
(b) Fine discretization of the potential contact region
To compare the effects of discretization in potential contact region, three kinds of fine element size in potential contact region have been given, that is, element size = 0.03 m, 0.05 m, and 0.1 m. Consequently, the results of proposed method have been listed in Table 1. It is shown that the half-width of contact length is about 0.2 m, and the gap is not sensitive in element size. Because the analytic solutions exist, the result of proposed method for element size 0.03 m has been compared with them in Figure 9. According to the formula, the half-width of the contact is equal to 0.197 m.
Figure 10 plots the contact pressure distributions along the contact length, and clearly there is excellent agreement between the analytical solution and the proposed method. It is also found that when the element size in potential contact region is smaller, the result for contact pressure is more accurate.
In general, the proposed method can give the satisfactory results. To obtain accurate results, the small element size in potential contact region is recommended.
Example 2 (laminated beam problem). Considering a laminated beam problem [42] as shown in Figure 11, the size of every beam is 10 × 1 × 1 m. The parameters of material are as follows: Young’s modulus MPa and Poisson’s ratio . There are two cases to consider. In case 1, the point has a downward concentrate force KN, while in case 2 the force is upward. The problem is considered under plane stress condition.
In order to simulate the concentrate load , a distributed load is applied vertically along a very small element near the point . To investigate the accuracy, the displacements of points , , and have been selected. Figure 12 shows the deformation under different load, and the results have satisfactory agreement with the results by Zheng et al. [42] and Li [43]. Table 2 lists displacements of specified three points.
(a) Downward load
(b) Upward load
5. Conclusions
Frictionless contact problems in two-dimensional space are formulated by complementarity theory, where the system of equations is established by the nonlinear complementary functions and boundary integral equations. This algorithm by BEM is established. The accuracy and effectiveness of the method have been demonstrated by two numerical examples, and the effect of discretization has also been studied in the Hertzian contact problem. The results show that this technique is very competitive and elegant.
Several extensions of the current work are possible. The presented method can be easily extended to contact problems with friction or involving inelastic materials.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge the financial support of Natural Science Foundation of China (no. 41202226) and National 973 Project of China (no. 2011CB013505).