Mathematical Problems in Engineering
Volume 2015 (2015), Article ID 907620, 8 pages
http://dx.doi.org/10.1155/2015/907620

Research Article

Multibody Finite Element Method and Application in Hydraulic Structure Analysis

Chao Su, Yebin Zhao, and Yusong Jiang

College of Water Conservancy and Hydropower, Hohai University, Xikang Road, Nanjing 210098, China

Received 18 September 2014; Accepted 15 December 2014

Academic Editor: Chenfeng Li

Copyright © 2015 Chao Su et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 as

where 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 that

with

where 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:

Figure 1: Schematic of two blocks initially in contact.

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 as

where 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 have

where , .

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 by

where is the displacement vector of node , is the rigid displacement vector, is the rotation angle vector, and is the radius vector.

Figure 2: Displacement of interface node caused by rigid displacement.

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 as

where 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,

with

where , , and are, respectively, the flexibility submatrix, the interface force, and the displacement components of node in , coordinates which are rotated .

Figure 3: Coordinate transformation.

Coordinate transformation of the continuity equations can be expressed as

with

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.

Table 1: Determination conditions of interface state.

State
Determination conditions

Fix
and and

Free

Forward slip
and and

Reverse slip
and and

Figure 4: Interface states.

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/mm². 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/mm², 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.

Table 2: Results comparison of Example 1.

Calculation parameters
Hertz solution
ANSYS solution
Present

(N/mm²)
1698.1
1617.7
1619.9

(mm)
1.1997
1.1609
1.1069

(mm)
0.4181
0.4213
0.4202

Figure 5: Example 1: contact between two cylinders.

Figure 6: Example 1: finite element model.

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.

Figure 7: Example 2: frictional contact between cantilever beams (unit: m).

Figure 8: Example 2: gap distribution in normal direction of interface node pairs.

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.

Figure 9: Global finite element model.

Figure 10: Finite element model of chambers.

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)

Figure 11: Horizontal displacement vector maps: (a) unit 2 and (b) unit 9, respectively (unit: m).

(a)

(b)

Figure 12: Deformed mesh: (a) unit 2 and (b) unit 9, respectively (50 of deformed scaling).

(a)

(b)

Figure 13: Major primary stress vector maps: (a) unit 2: (b) unit 9, respectively (unit: Pa).

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).

References

Y. M. Cheng, “Advancements and improvement in discontinuous deformation analysis,” Computers and Geotechnics, vol. 22, no. 2, pp. 153–163, 1998.

G. H. Shi, Discontinuous deformation analysis-a new numerical model for the statics and dynamics of block system [Ph.D. thesis], Department of Civil Engineering, University of California, Berkeley, Calif, USA, 1988.

G. H. Shi and R. E. Goodman, “Generalization of two-dimensional discontinuous deformation analysis for forward modelling,” International Journal for Numerical & Analytical Methods in Geomechanics, vol. 13, no. 4, pp. 359–380, 1989.

Q. H. Jiang and M. R. Yeung, “A model of point-to-face contact for three-dimensional discontinuous deformation analysis,” Rock Mechanics and Rock Engineering, vol. 37, no. 2, pp. 95–116, 2004.

Y. Zhang, Q. Xu, G. Chen, J. X. Zhao, and L. Zheng, “Extension of discontinuous deformation analysis and application in cohesive-frictional slope analysis,” International Journal of Rock Mechanics & Mining Sciences, vol. 70, pp. 533–545, 2014.

N. Sitar, M. M. MacLaughlin, and D. M. Doolin, “Influence of kinematics on landslide mobility and failure mode,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 131, no. 6, pp. 716–728, 2005.

J.-H. Wu, “Seismic landslide simulations in discontinuous deformation analysis,” Computers and Geotechnics, vol. 37, no. 5, pp. 594–601, 2010.

Y. Zhang, G. Chen, L. Zheng, and Y. Li, “Numerical analysis of the largest landslide induced by the Wenchuan earthquake, May 12, 2008 using DDA,” in Earthquake-Induced Landslides, pp. 617–626, Springer, Berlin, Germany, 2013.

W. Zhu, Q. Zhang, and L. Jing, “Stability analysis of the shiplock slopes of the Three Gorges project by three-dimensional FEM and DEM techniques,” in Third International Conference on Analysis of Discontinuous Deformation—From Theory to Practice (ICADD-3), Vail, Colorado, June 3-4, 1999, B. Amadei, Ed., pp. 263–272, American Rock Mechanics Association, Alexandria, Va, USA, 1999.

J. T. Kottenstette, “DDA analysis of the RCC modification for Pueblo Dam,” in Proceedings of the 3rd International Conference on Analysis of Discontinuous Deformation—From Theory to Practice (ICADD-3 '99), B. Amadei, Ed., pp. 127–132, American Rock Mechanics Association, Vail, Colo, USA, 1999.

D. M. Doolin and N. Sitar, “Displacement accuracy of discounting deformation analysis method applied to sliding block,” Journal of Engineering Mechanics, vol. 128, no. 11, pp. 1158–1168, 2002.

D. Bakun-Mazor, Y. H. Hatzor, and S. D. Glaser, “Dynamic sliding of tetrahedral wedge: the role of interface friction,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 36, no. 3, pp. 327–343, 2012.

Y. Zhang, X. Fu, and Q. Sheng, “Modification of the discontinuous deformation analysis method and its application to seismic response analysis of large underground caverns,” Tunnelling and Underground Space Technology, vol. 40, pp. 241–250, 2014.

Y. M. Cheng, “Advancements and improvement in discontinuous deformation analysis,” Computers and Geotechnics, vol. 22, no. 2, pp. 153–163, 1998.

T. C. Ke, “The issues of rigid-body rotation in DDA,” in Proceedings of the 1st International Forumon Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, pp. 318–325, Berkeley, Calif, USA, 1996.

D. M. Doolin and N. Sitar, “Time integration in discontinuous deformation analysis,” Journal of Engineering Mechanics, vol. 130, no. 3, pp. 249–258, 2004.

A. R. Keneti, A. Jafari, and J.-H. Wu, “A new algorithm to identify contact patterns between convex blocks for three-dimensional discontinuous deformation analysis,” Computers and Geotechnics, vol. 35, no. 5, pp. 746–759, 2008.

S. Amir Reza Beyabanaki, A. Jafari, S. Omid Reza Biabanaki, and M. Ronald Yeung, “Nodal-based three-dimensional discontinuous deformation analysis (3-D DDA),” Computers and Geotechnics, vol. 36, no. 3, pp. 359–372, 2009.

J. Liu, Z. Nan, and P. Yi, “Validation and application of three-dimensional discontinuous deformation analysis with tetrahedron finite element meshed block,” Acta Mechanica Sinica, vol. 28, no. 6, pp. 1602–1616, 2012.

A. Jafari, M. Khishvand, and H. Rahami, “Developing an algorithm for reconstruction blocky systems in discontinuous media: three-dimensional analysis,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 37, no. 7, pp. 661–684, 2013.