Abstract
Multibody finite element method is proposed for analysis of contact problems in hydraulic structure. This method is based on the block theory of discontinuous deformation analysis (DDA) method and combines advantages of finite element method (FEM) and the displacement compatibility equation in classical elastic mechanics. Each single block is analyzed using FEM in corresponding local coordinate system and all contacting blocks need to satisfy the displacement compatibility requirement between any two blocks in a blocky system. It is proved that this method is very efficient and practical to overcome the limitations in DDA method when tackling contact problems, such as the overlap problem and the equal strain assumption. In this paper, detailed theoretical basis and formulations are given. Two numerical examples are performed to verify the proposed method successfully. Furthermore, this method is adopted to study the stability issues of underground houses of a large hydropower station.
1. Introduction
The analysis of contact problems is of great importance in engineering applications. Much attention has been devoted to the development of various numerical methods for analyzing contact problems. The essence of this problem belongs to discontinuous deformation problem, where discontinuities are dominant in the analysis. Most of traditional numerical methods used in engineering are based on continuous model and do not fit contact problems apparently [1]. Discontinuous deformation analysis (DDA) is a kind of discrete numerical method for discontinuous blocky system that was proposed by Shi [2, 3]. DDA has the advantages of finite element method (FEM) and distinct element method (DEM); it can be used to analyze jointed rock mass behavior, such as dislocations, sliding, and rotations [4]. Since it was proposed, many developments and applications have been proposed to make it more practical in engineering applications [5]. In the past twenty years, DDA has been widely used to analyse static and dynamic problems, including landslides [6–8], dams [9, 10], dynamic block [11–13], and others. Although many modifications and improvements have been performed to the original DDA [14–20], two limitations still exist. First, normal and tangential stiff springs are applied to contact problems in the original DDA. Normal or tangential stiffness is difficult to be determined and directly affects the results. This kind of approximation method is not rigorous in theory and may lead to interpenetration between blocks and convergence problems in algorithm. Secondly, the applicability is greatly weakened because of the equal strain assumption. Due to this the ability of the original DDA to simulate complex deformation of blocky system is limited.
The purpose of this paper is to present a new numerical method—multibody finite element method (multibody FEM) that can effectively simulate the contact problems in hydroengineering. Based on the block theory of DDA, each block can move and deform independently, and blocks contact with each other on their boundary. The discontinuous deformation problem is changed to be contact problems between blocks in the proposed method. Furthermore, multibody FEM combines the basic theory of FEM and the displacement compatibility in classical mechanics. Each single block is analyzed with FEM in local coordinate system and the contacting blocks satisfy the displacement compatibility in blocky system. The interface force between blocks, deformation of each block, and the whole motion of blocks are considered as variables governed by the displacement compatibility and the equilibrium equations of blocky system. In this way, normal and tangential stiffness are not required in algorithm and the overlap problem is successfully avoided. Moreover, the accuracy of the block internal deformation is effectively improved.
The outline of the paper is organized as follows. Section 2 presents the theoretical basis and fundamental formulas of multibody FEM for the analysis of contact problems. The judgment of the interface states used in this study will be given in Section 3. Two numerical examples are performed to show the validity of the proposed method in Section 4. In Section 5, this method is used to study the stability of underground houses of a large hydropower station. Finally in Section 6, a conclusion concerning the proposed method is given.
2. The Multibody FEM
In this section, discussion is limited to two-dimensional contact blocks with linearly elastic materials and static formulations. Of course, the contact model can be extended to three-dimensional multiple blocks contact problems. We assume that blocks yield the assumption of small deformation, and the initial state of the interface is continuous without external force; that is, there is no gap between blocks. Further, it is assumed that any interface can be modelled as a sequence of node pairs, and the interface node can be scheduled in advance.
2.1. Finite Element Formulations
The popular FEM has been used extensively for scientific computing in engineering. The global set of equilibrium equations of two contacting elastic blocks can be simply expressed aswhere is the global stiffness matrix, is the displacement vector, and and are, respectively, the interface force vector and the external force vector. It is worth noting that (1) cannot be directly solved because and are both unknown. The key point of multibody FEM is to determine the interface force vector and make it meet all requirements of the interface states. After the finial interface force is obtained, vector and other quantities such as strains and stresses can be computed in the blocky system. By using an elementary calculation for (1), it is obtained thatwithwhere is the flexibility matrix and is the gap vector resulting from the application of the external loads .
Figure 1 depicts the contact model between two blocks and , and is assumed as a rigid foundation. Applying (2) to the contact model as suggested in Figure 1, then equations can be obtained:
It is important to be pointed out that may be in a stable equilibrium state even if it does not have enough boundary displacement constraints in blocky system. However, may show rigid displacement if it is out alone. When there is a rigid displacement of , the global stiffness matrix of is singular. Additional displacement constraint should be applied to eliminate the singularity of ; (4) can be rewritten aswhere is the rigid displacement matrix of and is the additional displacement constraints imposed on .
The relative displacement between and is defined by
According to Newton’s third law, it is obvious that . Combining (5), (6), and (7), we havewhere , .
It is assumed that the initial state of interface is continuous; that is, . The continuity equations of multibody FEM can be obtained:
2.2. The Situation of Different Rigid Displacement
The displacement of the interface node (Figure 2) caused by rigid displacement is defined bywhere is the displacement vector of node , is the rigid displacement vector, is the rotation angle vector, and is the radius vector.
When is completely free, the displacement , of node and the displacement of node are known in advance. Let ; the displacement of any node can be expressed as is expressed as
When has a translational degree of freedom, the displacement of freedom degree is . Let ; we have
When has a rotational freedom around the center , the displacement of node perpendicular to the radius is . Let ; the displacement of any node is then given by
2.3. Coordinate Equations
In order to solve , the equilibrium equations are given by
The matrix form of (15) can be written aswhere is the equilibrium matrix and is the external force matrix of .
When is completely free, we obtain
In the case, has a translational degree of freedom; we can get
In the case, has a rotational freedom around the center ; we have
The coordinate equation of multibody FEM combines the rigid displacement and overall balance conditions can be expressed as
2.4. Coordinate Transformation
Normal and tangential local coordinates of the contact boundary are introduced to modify the flexibility matrix in the proposed method. For interface node as shown in Figure 3,withwhere , , and are, respectively, the flexibility submatrix, the interface force, and the displacement components of node in , coordinates which are rotated .
Coordinate transformation of the continuity equations can be expressed aswith
Here the subscript denotes the number of the node pairs. Equation (23) is also applicable to the problems which contain rigid displacement.
3. Interface States Judgment
The judgment of the interface state is an important step in contact problems. Assume that all the interface node pairs are initially fixed. According to this assumption, the continuity equation (9) should be modified when the interface node pairs are in discontinuous status. Let , denote the unit normal and tangential vector, respectively. The attention can be focused on a single node pair constituted by nodes and on and as described in Figure 4; the interface state judgment is given in Table 1.
At each iteration, the interface force is chosen as the iterative variable, and according to the determination conditions in Table 1, the state of the interface node pairs can be determined. If the current interface state is not in accordance with the determination condition, and in (9) should be modified until all the conditions are satisfied. The finial interface force can be obtained after all the steps are performed. Using (1), the displacement can be computed in the blocky system.
4. Numerical Examples
Based on the theory discussed above, a program is developed in this paper. Two typical numerical examples are solved to verify the accuracy of multibody FEM. The results of the proposed method are compared with the results obtained using the general purpose finite element software ANSYS in which the Lagrange multiplier method is used. The first example is a classical Hertz contact problem without friction. The second problem deals with frictional contact between two cantilever beams.
Example 1 (Hertz contact problem). The first validation example is a classical Hertz contact between two cylinders. This frictionless contact problem is selected from the help file of ANSYS. Two long individual cylinders of radius and , respectively, in frictionless contact with their axes parallel to each other are pressed together with a force per unit length as shown in Figure 5. According to the Hertz theory, the semicontact length , , the maximum contact pressure , and the approach distance of theoretical solution refers to the help file of ANSYS. Here , , , and are Young’s modulus and Poisson’s ratio for cylinder, respectively.
The material parameters of Example 1 are GPa, , mm, GPa, , mm, and loading N/mm2. Because of symmetry, only a quarter of the body was modeled as a plane strain problem, the finite element model of this example as shown in Figure 6.
According to the proposed method, the results of semicontact length, maximum contact pressure, and approach distance are mm, N/mm2, and mm, respectively. Table 2 is the comparison of Hertz theory, commercial ANSYS, and the presented method. It is worth noting that Lagrange multiplier method is a more accurate method. From this table it is shown that the solution agrees very well with those obtained by ANSYS.
Example 2 (frictional contact between cantilever beams). This example is a well-known problem about the frictional contact between beams. A cantilever beam is superposed with two identical beams; the model size is 10 m × 2 m × 1 m. A distributed force N/m is applied at the bottom of the beam as shown in Figure 7. The material properties used for the beams are Young’s modulus GPa, Poisson’s ratio , cohesion , and the coefficient of friction .
The potential contact zone covers 105 node pairs and the initial state of all the node pairs was assumed to be already in contact. The result of maximum normal gap is 5.36 × 10−3 mm. Comparison has been made with normal gap distribution of interface node pairs along the contact length in Figure 8. It can be seen that the solution agrees well with ANSYS.
Two numerical examples have been discussed in Section 4 to verify the accuracy of the proposed method. The results are in a good agreement with the Lagrange multiplier method used in commercial ANSYS. It is to be noted that no parameters such as normal and tangential stiffness are required in this algorithm. Furthermore, the proposed method has a faster convergence speed; generally, each incremental step can be obtained after about 3–7 iterations.
5. Application
Multiblock contact problems exist widely in hydraulic structure, such as dams with structure joints, soil-structure interaction, and structure plane of rock. In this study, the multibody FEM is applied to the numerical analysis of an underground powerhouse with structural planes. The model contains the main powerhouse, main transformer chamber, and tail surge chamber, and there are 9 hydropower units sequentially along the powerhouse’s axis. The axis of the underground powerhouse is located at N35°E. Around the house, there is a layer fault zone located at N40°~50°/SE∠15°~25°; the thickness is about 20~30 cm. The global finite element model is modeled by 68400 8-node hexahedron elements as shown in Figure 9, and the structural planes are assumed to have no thickness. The material parameters of rock and structural planes are Young’s modulus GPa, Poisson’s ratio , density , coefficient of friction , cohesive strength MPa, and initial tensile strength . Figure 10 shows the finite element model for chambers; the potential contact zone covers 1715 interface node pairs. The thick black line in Figures 9 and 10 shows the position of the structural plane.
The maximum horizontal displacement of rock masses surrounding the caverns after excavation occurs at unit 2 of the upstream side wall and unit 9 of the downstream side wall. Figures 11 and 12, respectively, show the horizontal displacement vector and deformed mesh of units 2 and 9. Due to the layer fault zone, the stress concentration appears near the structural plane. The major primary stress vector of units 2 and 9 is shown in Figure 13. The results obtained by the presented method are reasonable and believable.
(a)
(b)
(a)
(b)
(a)
(b)
6. Conclusions
A novel multibody finite element method is proposed for the challenging contact problems in hydraulic structure. Based on the block theory of DDA and FEM, discontinuous deformation problem is transformed to contact problems between blocks to be simulated and solved in this paper. It is proved that this method is very efficient and practical to overcome the limitations of DDA method when tackling the contact problems. The theoretical basis and fundamental formulas of the proposed method are set up. Two numerical examples and application have been studied in detail. Performance analysis and simulation results show the feasibility and effectiveness of the proposed method.
It should be noted that discussion about contact blocks is limited to linearly elastic materials and only point-to-point contact model is used in the proposed method. Further studies and developments are going to be undertaken for multibody FEM.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The present work is supported by the National Natural Science Foundation of China (nos. 51279053, 51109067) and the Fundamental Research Funds for the Central Universities (no. 2014B36514).