Abstract
A novel mixed finite element method is proposed for static and dynamic contact problems with friction and initial gaps. Based on the characteristic of local nonlinearity for the problem, the system of forces acting on the contactor is divided into two parts: external forces and contact forces. The displacement of structure is chosen as the basic variable and the nodal contact force in contact region under local coordinate system is selected as the iteration variable to confine the nonlinear iteration process in the potential contact surface which is more numerically efficient. In this way, the sophisticated contact nonlinearity is revealed by the variety of the contact forces which are determined by the external load and the contact state stick, slip, or separation. Moreover, in the case of multibody contact problem, the flexibility matrix is symmetric and sparse; thus, the iterative procedure becomes easily carried out and much more economical. In the paper, both the finite element formulations and the iteration process are given in detail for static and dynamic contact problems. Four examples are included to demonstrate the accuracy and applicability of the presented method.
1. Introduction
Various types of static or dynamic contact problems are widely encountered in many fields, such as mechanical, civil, and hydraulic engineering [1, 2]. The chief difficulty in simulating those problems lies in the inherent and strong nonlinearity due to the fact that both the region of contact and the contact forces distribution are unknown prior to the analysis. The friction effects and initial gaps make the problem more difficult to handle. So problems involving contact with friction and initial gaps are among the most challenging issues in solid mechanics and at the same time of crucial practical importance in many engineering branches.
Although several formulation procedures have been developed for contact problems, such as Lagrange multiplier method [3–5], penalty method [6–9], mathematical programming technique [10–12], or complementary method [13, 14], there is still a great need for a reliable, accurate, and more effective method. Since the area of contact is not known prior to the application of loads, the contact problems inherently impose nonlinear boundary conditions. Finite element technologies based on above methods have been widely used [15]; however, conventional methods have a fundamental drawback from a phenomenological point of view. Notice that the actual contact occurs only in a local part of the boundary, and its state is changed during time increments or iterations; therefore, deriving the governing equations on the basis of the entire domain and performing the related analyses may raise the computational inefficiency. Generally, the forgoing methods treat all the freedom of degrees (DOFs) of the system as the iterative variables and they are thus computationally expensive. Numerous stiffness matrices decomposition is required, especially for a dynamic analysis to track down the long-term structural response. But the important characteristic of contact problem is it belongs to local nonlinear problem; that is, the contact nonlinearity only comes from the potential contact region which is always only a small part relatived to the whole model, so those kinds of method are obviously very uneconomical in a computational sense.
Some research efforts have been focused on finding effective numerical techniques for those complex contact problems, but it is still under development. For instance, based on the characteristic of local nonlinearity of this problem, the nonlinear iterative process can be limited on the potential contact region through condensation of DOFs, such that it is computational saving in some sense [16, 17]. However, it is dissatisfactory because the flexible matrix obtained through this kind of methods is a full matrix which is not easy to handle. In the case of multibody contact problem or a lot of nodes which likely come into contact, efficiency is still a problem and remains a challenge.
The objective of this study is to develop a new and highly effective iterative method—mixed finite element method, a node-to-node contact algorithm for both static and dynamic contact problems with friction and initial gaps. Based on the characteristic of local nonlinearity of this problem, only the nodal contact forces on potential contact region are chosen as the iterative variables in local coordinate system. Therefore, the nonlinear iterative process is only limited on the potential contact surfaces and could be very economic. In this way, the sophisticated contact nonlinearity is revealed by the variety of the contact forces which are determined by the external load and the contact state stick, slip, or separation. The iterative process both for contact forces and contact state is carried out during the simulation. Some parameters which is difficult to determine or sensitive to the results, such as normal or tangential stiffness and penalty factors, are not required in this algorithm. The overlap problem which often appeared in other methods is successfully avoided. Moreover, in the case of multibody contact problems, the flexibility matrix is symmetric and sparse; thus, the iterative procedure becomes easier to be carried out and much more economical.
At the beginning of the analysis, the global stiff matrix is decomposed only once, and then each component of the flexibility matrix is calculated by the back-substitution in the local coordinate system defined in the surface where contact is likely to take place. The flexibility matrix is obtained by assembling its valid components according to the contact status during the iteration, while the contact status will be modified by the directly solved contact forces in terms of the contact principles. After convergences of both contact status and contact forces are achieved, the algorithm is started again with a new increment of external information. The problems considered here involve the contact of linearly elastic bodies on potential contact surface for which a generalized Mohr-Coulomb’s law of friction is assumed. In dynamic analysis, the time integration of the equilibrium equations is performed using the generalized Newmark method [18].
The rest of the paper is organized as follows. In the next section, the mechanical model for contact problems is described and the contact principles are summarized. Sections 3 and 4 present the implementations of static and dynamic contact problems respectively, including the finite element formulations and the iteration process. Some numerical examples and applications in engineering will be given in Section 5. In Section 6, a concise conclusion concerning the above method is given.
2. Mechanical Description
This paper is limited to contact of two or more bodies with linearly elastic materials and small deformations. With no loss of generality, consider two elastic bodies , (shown in Figure 1) which are brought into contact by the external force . Each boundary of the two bodies is divided into three disjoint parts : where the traction boundary condition is given, , where the displacement is prescribed, and , which is the potential contact region, taken as sufficiently large to include the real contact surface. Here the superscript denotes the two bodies, respectively. We now focus our attention on an arbitrary node pair constituted by points 1 and 2 on and , respectively. The two nodes of the pair are assumed to be close to each other in any case that common normal and tangential directions can be determined. Furthermore, a right-handed orthogonal local coordinate system centered at point 1 or point 2 in terms of undeformed geometry is defined. is taken to be normal to the contact surface and and are two tangent directions. It is convenient to work in this local coordinates so that both displacement and contact forces are defined with respect to the boundary normal. Other notations used in this section are given as follows. : contact forces on . According to Newton’s third law, it is obvious that ; , , : component of contact forces in local coordinate , respectively; : displacement vector on ; : the gaps between and measured in the normal direction and , where is the normal direction to the contact region; is the initial normal gaps between and .
Contact problem is a fairly complex mechanical phenomenon in practice. When it is described as a mechanical model, three principles must always be satisfied to obtain the correct solution.
(1) Nonpenetration Condition. Following this condition, no material particle will be allowed to penetrate into the surface of the opposing body. The geometric compatibility condition has been generally described in [19]. The nonpenetration condition can be stated as
(2) Normal Traction Condition. The normal traction between two bodies must be compressive or zero once it exceeds the specified tensile strength, , for the first times. If is regularized when the specified tensile strength has been exceeded once or there is no specified tensile strength at all, this condition can be simply given as
(3) Frictional Condition. This condition represents the relationship between the normal traction and tangential traction. When the total tangential contact forces, denoted as , reach or exceed the maximum resistance in a node pair, slipping will take place in the same direction as . When the contact state is stick, the friction between the two contact boundaries becomes effective in resisting the slipping motion. In the presented paper, this is described by the Mohr-Coulomb friction law: where is the tributary area of point 1 or point 2; is the friction coefficient; is the cohesive strength.
When the three contact statuses are involved, the contact conditions presented above can be summarized as shown in Table 1.
In the following two sections, the algorithm for the numerical implementation of the contact problem will be discussed separately for the static and the dynamic cases; both the mixed finite element formulations and the iteration process are included. In static analysis, the time step is the pseudovariable used to determine the current level of the applied load, while for dynamic case the time step refers to the physical time value.
3. Static Contact Problem
3.1. Finite Element Formulations
Assuming that the analysis of th step has been finished, the incremental static equilibrium equation for th step is stated as follows by performing a finite element discretization of each body and . Consider where is the global stiffness matrix, is the vector of displacement increment at step , is the vector of total external load at step , is the vector of total contact force at step , is the strain-displacement matrix, is the Cauchy stress tensor, and is the vector of incremental contact force. Clearly, the first term of the right-hand side (RHS) of above equation, that is, content in the bracket, stands for the external load increment of current load step, since the contribution of to the system is already included into .
If a typical decomposition of global stiffness matrix is performed at the beginning of the analysis for one and only one time, (4) could be rewritten as
If is defined as
and the matrix is introduced into above equation, (6) is rewritten as where the matrix is the flexibility matrix defined on the potential contact boundary . An arbitrary component of , , represents flexibility coefficient corresponding to the displacement at the freedom due to a unit force at the freedom . Here or is only limited in the freedom of the region where contact is likely to take place.
It is important to be pointed out that (6) and (7) and the following equations in this section are only performed for the DOFs on the potential contact surface, not for all the DOFs of the whole system.
Applying (7) to points 1 and 2 of a given node pair, respectively, gives
According to Newton’s third law, it is obvious that , and moreover, incorporating flexibility matrices and into , and subtracting (8) from (9), we obtain
Equation (10) is the finite element compliance equation of the mixed finite element method for the static contact problems with friction and initial gaps. In this equation, the second term in RHS, , stands for the difference of incremental displacement only due to the external load increment, as can be seen easily from (6). Therefore, it is nothing to do with the current contact state and can be obtained by back-substitution directly. While as for the first term in RHS, , it denotes the difference of incremental displacement induced by both the external load increment and the contact force increment. From this point of view, it can be seen that both the RHS and the left-hand side (LHS) of the above equation are associated with the contact state. An iterative method taking into account different contact state is necessary to solve (10), and this will be given next in detail.
3.2. Iteration Process
Considering the contact conditions stated in Section 2 and the finite element formulations derived in first part of this section, the algorithm for static contact problems can be summarized as follows.(1)Assemble the global stiffness matrix from the element level and decompose it. Obtain each component of flexibility matrix that might be used during the iteration through the unit force algorithm.(2)Assuming that the analysis of the th step has been finished, one comes into next time step . Update the contact forces and contact state; that is, set the results obtained in last step as the initial condition. can be solved simply by back-substitution of (6) with the known external force increment.(3)Assuming the th iteration is finished, one comes into next to perform iteration both for contact state and contact forces. Transmit the contact state and contact forces at th iteration to current iteration step. Different treatments will be carried out according to different contact state as follows.(a)If the current contact state is separated, calculate the value of gaps between two bodies following the expression , where comes from the results of last iteration step. Here and the remainders of this paper, the superscript denotes the iteration step while the subscript stands for the loading increment step. If is less than or equal to zero, it implies that the bodies are brought into contact now and produce an overlap; therefore, change the contact state from separation to stick. Moreover, in the case of is less than zero, set to zero to satisfy the nonpenetration condition stated in Section 2.(b)If the current contact state is slip or stick, obtain the normal stress component by . If is larger than , separation will occur according to the normal traction condition, so the current contact state is changed to separation and is set. At last is set since the specified tensile strength can be exceeded only once.(c)If the current contact state is slip or stick after (a) and (b), get the total tangential contact forces first and then employ the Mohr-Coulomb friction law to determine whether slipping will occur or not, according to the frictional condition. Change the contact state from stick to slip if slipping will take place. At last, both three freedoms of each node pair which is under stick state and normal freedom of each node pair which is under slip state are marked to determine the size of the flexibility matrix in the current iteration step.(4)Perform the convergence check for contact state. icsta is used as the sign to indicate whether the contact state iteration is convergence or not, 0 for not and 1 for yes. If the contact state of all the pair nodes is identical with the state obtained in last iteration step, icsta is set to 1 and turn to step ; otherwise, icsta is set to 0 and go to next step. It is to be noted that if icsta = 0 is satisfied in current iteration step, assembling and decomposition of matrix are not required when we continue our analysis for contact forces convergence since the convergence for contact state has been achieved.(5)Assemble the flexibility matrix and decompose it. Only the freedoms marked in step 3 are included.(6)Solve the finite element compliance equation (10) by back-substitution to get iterative increment of contact forces : (7)The contact forces increment for th loading step after current iteration can be obtained: (8)Perform the convergence check for contact forces. The contact forces convergence criterion employed here is where is a predefined tolerable limit. Set icfor, a sign used to indicate whether the contact force iteration is convergent or not, to 1 for (13) is satisfied and icfor = 0 for others. If both the convergence for contact state and contact forces are achieved, that is, icsta = 1 and icfor = 1, go to step 9, otherwise continue.(9)Obtain the new iterative increment of displacement according to equilibrium equation (5): Then turn to step 3 to continue next iteration for contact state and contact forces.
Apply the actual contact forces increment of th loading step into equilibrium equation (4) to solve the actual displacement increment. Go to step 2 until all the time steps are finished.
4. Dynamic Contact Problem
4.1. Finite Element Formulations
The basic ideas adopted here are quite similar to those stated for static problem except for some little difference. The previously derived algorithm is now modified to account for dynamic contact problem using the popular predictor-corrector Newmark time integration scheme.
The incremental dynamic equilibrium equation for th step can be expressed as where is the global mass matrix; is the global damping matrix; here the proportional Rayleigh damping assumption is adopted; , , , is the vector of displacement increment, velocity increment, and acceleration increment at time , respectively. Other notations used here and below are the same with those in last section.
Following the generalized Newmark time integration scheme, the predicting and correcting stage can be written in following two expressions, respectively: where the superscript refers to the predicted quantities and , are Newmark parameters for time integration that control the accuracy and stability of the solution. In the subsequent calculations, and are employed. For this condition, the Newmark method is unconditionally stable and it introduces no artificial damping in linear vibration analysis.
By plugging (16) and (17) into (15), the equation with the variables of the acceleration increment is obtained as where , are the effective stiffness matrix and effective external load increment and can be expressed as follows:
It is to be noted that (18) is similar to (4) except for the fact that the basic unknown variable is the increment of acceleration instead of increment of displacement. So the previously derived formulations for the static contact algorithm remain valid if some necessary modifications are carried out.
Similar to static case, if the effective stiffness matrix is decomposed at the beginning of the analysis, for one and only one time, (18) could be rewritten in
Introducing the matrix into above equation results in where the matrix is the flexibility matrix which is defined on the potential contact boundary . An arbitrary component of , , represents flexibility coefficient corresponding to the acceleration, instead of displacement for static case, at the freedom due to a unit force at the freedom .
Considering (16), (17), and , if is defined as
one obtains
It is important to be noted again that (23), (24), and the following equations in this section are only performed for the DOF on the potential contact surface, not for all the DOFs of the whole system.
Applying above equation (24) to points 1 and 2 of a given node pair gives
Same deduces as static case result in the element compliance equation for the dynamic contact problems
It is exactly the same with (10) except for that there is a constant coefficient in RHS. A similar iteration process to solve the dynamic contact problems with friction and initial gaps will be given below briefly.
4.2. Iteration Process
(1)Assemble the effective stiffness matrix according to (19) from the element level and perform decomposition. Obtain each component that might be used during the iteration of flexibility matrix through the unit force algorithm.(2)Assuming that the analysis of the th step was achieved, one comes into next time step . Update the contact forces and contact status. Apply the external force increment to the system to obtain following (23) by back-substitution.(3)Assuming that iteration was finished, one comes into next to perform iteration both for contact state and contact force. Transmit the contact state and contact force to current iteration step. Different treatments will be carried out according to different contact status as stated in Section 3.2.(4)Perform convergence check for the contact state. If the convergence for contact state is achieved, turn to step 6, otherwise continue.(5)Assemble the flexibility matrix and decompose it.(6)Solve the finite element compliance equation (26) by back-substitution to get iterative increment of contact force : (7)The contact force increment for th loading step is obtained as (8)Perform convergence check for contact forces. If both the convergence check for contact state and contact forces are achieved, go to step 9, otherwise continue.(9)Obtain the new iterative increment of acceleration according to following equilibrium equation: Then turn to step 3 to continue next iteration for contact state and contact forces.
Put the actual contact forces increment of th loading step into equilibrium equation (18) to get the actual acceleration increment . Obtain the displacement, velocity, and acceleration at time according to (17). Go to step 2 until all the time steps are finished.
5. Numerical Experiments
Several numerical experiments are performed to verify the ability and accuracy of the proposed method. These are as follows:(i)contact between cylinder and semispace (static Hertz problem),(ii)contact between two blocks,(iii)behavior of concrete block with two joints under dynamic concentrated load,(iv)seismic response of high arch dams with contraction joints.
Example 1 (contact between cylinder and semispace). In the first example, a classical static Hertz problem with no friction is included to assess the accuracy and reality of the algorithm proposed in this paper. An infinitely long circular cylinder with radius is pushed against semispace by a uniform vertical load , which is acting on in its neutral axis, as shown in Figure 2. Following the Hertz theory [20], the analytical results are half contact width and distribution of contact pressure along direction . Here , , , , is the modulus and Poisson’s ratio for cylinder and semispace, respectively. Owing to the symmetry, only half of the model was modeled as a plane strain problem. The model is discretized by four-node isoparametric elements. The potential contact zone covers 20 node pairs.
For computational purposes the following data will be adopted:
According to the presented method, the results of half contact width and maximum contact pressure are mm and N/mm2, respectively. Comparing these with the corresponding exact solution, mm, N/mm2, the difference is very small. Comparison made with contact pressure distribution along direction in Figure 3 shows good agreements to the exact Hertz solution.
Example 2 (contact between two blocks). Two rectangular blocks are superimposed onto each other. A concentrated force is applied at the center of the upper block. For convenient of comparison, the setup of this example, as shown in Figure 4, is the same as considered in [17, 21]. The material parameters of the two structures are identical and have the following values: Young’s modulus GPa, Poisson’s ratio , friction coefficient varies from 0.0 to 0.7, loading N, and thickness cm. Since the problem is symmetric about the central plane, only one-half of each body is modeled by 321 4-node elements as shown in Figure 5. Twenty contact node pairs were laid out to represent the interface between the two blocks. Initially, all the node pairs were assumed to be already in contact. The total load was applied to the system in 5 equal increments to reach solution convergence.
The same problems were previously analyzed in [17, 21] for . To verify the algorithm presented in this paper, the same problem is conducted with the same conditions. Comparison has been made with the normal pressure distribution along the contact surface in Figure 6. It can be seen that the solution agrees very well with those obtained by others.
The distribution of normal pressure along the contact surface for the contact system with different friction coefficient is shown in Figure 7. It is evident that the contact zone increases with increase in friction coefficient, while the maximum normal pressure decreases with increase in friction coefficient. This remark agrees well with reference [21].
Example 3 (behavior of concrete block with two joints under dynamic concentrated load). In this example, the behavior of concrete block with two vertical joints is analyzed under polygonal linear dynamic concentrated load :
The problem was meshed by 287 4-node isoparametric elements. The finite element model and its size were shown in Figure 8. For computations, the following data is assumed: amplitude of load N, Young’s modulus GPa, Poisson’s ratio , coefficient of friction , cohesive strength , initial normal gaps .
The time histories of displacement with respect to loading for the top node pairs of two joints were plotted in Figures 9 and 10. Compared with the results of the system for and , it is founded that the effect of friction not only decreases the amplitude of displacement but also deteriorates the nonlinearity of the relationship between displacement and loading. Although there is no exact continuum solution to be used for comparison, the results obtained following the presented method are reasonable and believable in a qualitative sight.
Example 4 (seismic response of high arch dams with contraction joints). The previously developed model has been applied to the seismic response of a high arch dam with contraction joints. The material properties are modulus of elasticity GPa, Poisson’s ratio , coefficient of friction , density , cohesive strength , initial tensile strength , and initial normal gaps . The dam is a 235 m high double-curvature arch dam. It has a crest length of 304.8 m, and the thickness of the center cantilever varies from 48 m at the base to 10 m at the crest. The dam consists of fourteen cantilevers separated by thirteen contraction joints. Finite element model of the dam and the layout of transverse joints are shown in Figure 11; hexahedral isoparametric element has been used in it. An artificial ground motion, with the peak acceleration 0.27 g and duration of 10 s, as shown in Figure 12, using the massless foundation with fixed boundary, is applied both in the stream and cross-stream direction. The opening and closing of the joints and its effects to the response of dam under the excitation of ground motion are studied. Considering free surface boundary, fluid-structure-interaction boundary, fluid absorb boundary, and bottom absorb boundary, the governing equation of the fluid domain can be discretized by the standard Galerkin weighted residual method. Combining with the general dynamic equilibrium equation for the solid domain, the static response of the dam-water foundation with upstream water level 202 m is obtained before determining that the nonlinear earthquake response can be solved.
Time histories of joint openings of points A and B, as marked in Figure 11, are plotted in Figure 13, while Figure 14 gives the time histories of the major principal stress of point C in the condition of contraction joints are considered or not.
It is found that once the tensile stress exceeds the tensile strength during the earthquake, the joints will open and consequently the tensile stress will be released. This is because large arch tensile stresses cannot be transferred across the contraction joints; hence, arch tensile stresses releasing and the internal forces redistributing from the arch action to cantilever action will take place. Some researchers have drawn the same conclusion about this, as in [22–24].
6. Conclusions
A novel mixed finite element method is proposed for static and dynamic contact problems with friction and initial gaps. Firstly, the mechanical model for three-dimensional static and dynamic frictional contact problems with initial gaps is presented and then the finite element compliance equation is derived. The iteration process is given in detail at last. The proposed method is first assessed through several selected sample numerical solutions to static and dynamic contact problems and is subsequently applied to the seismic response of high arch dams with contraction joints. It compares well to some exact solution or results obtained by others.
It is to be noted that, for a structure with several joints or cracks, for instance in Example 3 or 4, the system is divided into several substructures by joints. The substructure keeps its linearity during loading if the material is linear elastic one. Some methods couple all the DOFs of the joints together according to the condensation of DOF; thus, the flexible matrix obtained through those methods is a full matrix, which is not easy to handle. However, numerical studies have shown that it is only necessary to retain the coupling between two adjacent joints. In effect, each linear substructure is only coupled to its adjacent substructures through foundation. From this point of view, the flexibility matrix obtained in this paper is symmetric and sparse; thus, the iterative procedure becomes easier to be carried out and much more economical.
Conflict of Interests
The authors do not have any conflict of interests regarding the content of the paper.
Acknowledgments
The present work is supported by the National Natural Science Foundation of China (no. 51279050), National High-tech R&D Program of China (863 Program) (no. 2012BAK10B04), National Key Technology R&D Program in 12th Five-Year Plan (no. SS2012AA112507), and Non-Profit Industry Financial Program of MWR (Ministry of Water Resources of China) (no. 201301058).