# Multibody Finite Element Method and Application in Hydraulic Structure Analysis

Authors: Chao Su; Yebin Zhao; Yusong Jiang
Journal: Mathematical Problems in Engineering
 (2015)

---

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


---

## Body

## 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. Section2 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(1)KU=R+F,where K is the global stiffness matrix, U is the displacement vector, and R and F are, respectively, the interface force vector and the external force vector. It is worth noting that (1) cannot be directly solved because U and R are both unknown. The key point of multibody FEM is to determine the interface force vector R and make it meet all requirements of the interface states. After the finial interface force R is obtained, vector U 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(2)U=CR+UFwith(3)C=K-1,UF=K-1F,where C is the flexibility matrix and UF is the gap vector resulting from the application of the external loads F.Figure1 depicts the contact model between two blocks Ω1 and Ω2, and Ω2 is assumed as a rigid foundation. Applying (2) to the contact model as suggested in Figure 1, then equations can be obtained:(4)U1=C1R1+UF1,(5)U2=C2R2+UF2.Figure 1
Schematic of two blocks initially in contact.It is important to be pointed out thatΩ1 may be in a stable equilibrium state even if it does not have enough boundary displacement constraints in blocky system. However, Ω1 may show rigid displacement if it is out alone. When there is a rigid displacement of Ω1, the global stiffness matrix of Ω1 is singular. Additional displacement constraint should be applied to eliminate the singularity of Ω1; (4) can be rewritten as(6)U1=C1R1+UF1+CeUe,where Ce=K-1Ke is the rigid displacement matrix of Ω1 and Ue is the additional displacement constraints imposed on Ω1.The relative displacement betweenΩ1 and Ω2 is defined by(7)U*=U1-U2.According to Newton’s third law, it is obvious thatR1=-R2=R. Combining (5), (6), and (7), we have(8)U*=CaR+ΔF+CeUe,where Ca=C1+C2, ΔF=UF1-UF2.It is assumed that the initial state of interface is continuous; that is,U*=0. The continuity equations of multibody FEM can be obtained:(9)CaR=-ΔF-CeUe.


### 2.2. The Situation of Different Rigid Displacement

The displacement of the interface nodei (Figure 2) caused by rigid displacement is defined by(10)ui=uo+w×r,where ui is the displacement vector of node i, uo is the rigid displacement vector, w is the rotation angle vector, and r is the radius vector.Figure 2
Displacement of interface node caused by rigid displacement.WhenΩ1 is completely free, the displacement uA, vA of node AxA,yA and the displacement vB of node BxB,yB are known in advance. Let Ue=uAvAvBT; the displacement of any node ixi,yi can be expressed as(11)ui=uA-vB-vAxB-xAyi-yA,vi=vA+vB-vAxB-xAxi-xA.Ce is expressed as(12)Ce=⋯⋯⋯1yi-yAxB-xA-yi-yAxB-xA01-xi-xAxB-xAxi-xAxB-xA⋯⋯⋯.WhenΩ1 has a translational degree of freedom, the displacement of freedom degree y′ is vA. Let Ue=vA; we have(13)Ce=-sin⁡θcos⁡⁡θ⋯-sin⁡θcos⁡⁡θT.WhenΩ1 has a rotational freedom around the center Oxo,yo, the displacement of node AxA,yA perpendicular to the radius OA is vA. Let Ue=vA; the displacement of any node ixi,yi is then given by(14)ui=-cos⁡⁡θxA-xoyi-yovA,vi=cos⁡⁡θxA-xoxi-xovA,Ce=⋯cos⁡⁡θxA-xoyi-yocos⁡⁡θxA-xoxi-xo⋯T.


### 2.3. Coordinate Equations

In order to solveUe, the equilibrium equations are given by(15)∑Ω1X=0;∑Ω1Y=0,∑Ω1M=0.The matrix form of (15) can be written as(16)QR=Fe,where Q is the equilibrium matrix and Fe is the external force matrix of Ω1.WhenΩ1 is completely free, we obtain(17)Q=⋯10⋯⋯01⋯⋯yi-yA-xi-xA⋯,Fe=∑Ω1Fxi∑Ω1Fyi-∑Ω1Fxiyi-yA-Fyixi-xA.In the case,Ω1 has a translational degree of freedom; we can get(18)Q=⋯-sin⁡βi-θcos⁡⁡βi-θ⋯,Fe=-∑Ω1Ryicos⁡⁡θ-Rxisin⁡θ.In the case,Ω1 has a rotational freedom around the center Oxo,yo; we have(19)Q=⋯yi-yo-xi-xo⋯,Fe=-∑Ω1Ryiyi-yo-Rxixi-xo.The coordinate equation of multibody FEM combines the rigid displacement and overall balance conditions can be expressed as(20)CaCeQ0RUe=-ΔFFe.


### 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 nodei as shown in Figure 3,(21)C~i=TiCiTi-1;R~i=TiRi;U~i=TiUiwith(22)Ti=cos⁡⁡αi-sin⁡αisin⁡αicos⁡⁡αi,where C~i, R~i, and U~i are, respectively, the flexibility submatrix, the interface force, and the displacement components of node i in x′, y′ coordinates which are rotated αi.Figure 3
Coordinate transformation.Coordinate transformation of the continuity equations can be expressed as(23)C~a=TCaTT;R~=TR,Δ~F=TΔFwith(24)T=T10T2⋱0Tm.Here the subscriptm denotes the number of the node pairs. Equation (23) is also applicable to the problems which contain rigid displacement.


## 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(1)KU=R+F,where K is the global stiffness matrix, U is the displacement vector, and R and F are, respectively, the interface force vector and the external force vector. It is worth noting that (1) cannot be directly solved because U and R are both unknown. The key point of multibody FEM is to determine the interface force vector R and make it meet all requirements of the interface states. After the finial interface force R is obtained, vector U 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(2)U=CR+UFwith(3)C=K-1,UF=K-1F,where C is the flexibility matrix and UF is the gap vector resulting from the application of the external loads F.Figure1 depicts the contact model between two blocks Ω1 and Ω2, and Ω2 is assumed as a rigid foundation. Applying (2) to the contact model as suggested in Figure 1, then equations can be obtained:(4)U1=C1R1+UF1,(5)U2=C2R2+UF2.Figure 1
Schematic of two blocks initially in contact.It is important to be pointed out thatΩ1 may be in a stable equilibrium state even if it does not have enough boundary displacement constraints in blocky system. However, Ω1 may show rigid displacement if it is out alone. When there is a rigid displacement of Ω1, the global stiffness matrix of Ω1 is singular. Additional displacement constraint should be applied to eliminate the singularity of Ω1; (4) can be rewritten as(6)U1=C1R1+UF1+CeUe,where Ce=K-1Ke is the rigid displacement matrix of Ω1 and Ue is the additional displacement constraints imposed on Ω1.The relative displacement betweenΩ1 and Ω2 is defined by(7)U*=U1-U2.According to Newton’s third law, it is obvious thatR1=-R2=R. Combining (5), (6), and (7), we have(8)U*=CaR+ΔF+CeUe,where Ca=C1+C2, ΔF=UF1-UF2.It is assumed that the initial state of interface is continuous; that is,U*=0. The continuity equations of multibody FEM can be obtained:(9)CaR=-ΔF-CeUe.


## 2.2. The Situation of Different Rigid Displacement

The displacement of the interface nodei (Figure 2) caused by rigid displacement is defined by(10)ui=uo+w×r,where ui is the displacement vector of node i, uo is the rigid displacement vector, w is the rotation angle vector, and r is the radius vector.Figure 2
Displacement of interface node caused by rigid displacement.WhenΩ1 is completely free, the displacement uA, vA of node AxA,yA and the displacement vB of node BxB,yB are known in advance. Let Ue=uAvAvBT; the displacement of any node ixi,yi can be expressed as(11)ui=uA-vB-vAxB-xAyi-yA,vi=vA+vB-vAxB-xAxi-xA.Ce is expressed as(12)Ce=⋯⋯⋯1yi-yAxB-xA-yi-yAxB-xA01-xi-xAxB-xAxi-xAxB-xA⋯⋯⋯.WhenΩ1 has a translational degree of freedom, the displacement of freedom degree y′ is vA. Let Ue=vA; we have(13)Ce=-sin⁡θcos⁡⁡θ⋯-sin⁡θcos⁡⁡θT.WhenΩ1 has a rotational freedom around the center Oxo,yo, the displacement of node AxA,yA perpendicular to the radius OA is vA. Let Ue=vA; the displacement of any node ixi,yi is then given by(14)ui=-cos⁡⁡θxA-xoyi-yovA,vi=cos⁡⁡θxA-xoxi-xovA,Ce=⋯cos⁡⁡θxA-xoyi-yocos⁡⁡θxA-xoxi-xo⋯T.


## 2.3. Coordinate Equations

In order to solveUe, the equilibrium equations are given by(15)∑Ω1X=0;∑Ω1Y=0,∑Ω1M=0.The matrix form of (15) can be written as(16)QR=Fe,where Q is the equilibrium matrix and Fe is the external force matrix of Ω1.WhenΩ1 is completely free, we obtain(17)Q=⋯10⋯⋯01⋯⋯yi-yA-xi-xA⋯,Fe=∑Ω1Fxi∑Ω1Fyi-∑Ω1Fxiyi-yA-Fyixi-xA.In the case,Ω1 has a translational degree of freedom; we can get(18)Q=⋯-sin⁡βi-θcos⁡⁡βi-θ⋯,Fe=-∑Ω1Ryicos⁡⁡θ-Rxisin⁡θ.In the case,Ω1 has a rotational freedom around the center Oxo,yo; we have(19)Q=⋯yi-yo-xi-xo⋯,Fe=-∑Ω1Ryiyi-yo-Rxixi-xo.The coordinate equation of multibody FEM combines the rigid displacement and overall balance conditions can be expressed as(20)CaCeQ0RUe=-ΔFFe.


## 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 nodei as shown in Figure 3,(21)C~i=TiCiTi-1;R~i=TiRi;U~i=TiUiwith(22)Ti=cos⁡⁡αi-sin⁡αisin⁡αicos⁡⁡αi,where C~i, R~i, and U~i are, respectively, the flexibility submatrix, the interface force, and the displacement components of node i in x′, y′ coordinates which are rotated αi.Figure 3
Coordinate transformation.Coordinate transformation of the continuity equations can be expressed as(23)C~a=TCaTT;R~=TR,Δ~F=TΔFwith(24)T=T10T2⋱0Tm.Here the subscriptm 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 n, s denote the unit normal and tangential vector, respectively. The attention can be focused on a single node pair constituted by nodes i and i′ 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
R
n
i
≤
0 and Rsi≤μRni and Rsi>0
Free
R
n
i
>
0
Forward slip
R
n
i
≤
0 and Rsi>μRni and Rsi>0
Reverse slip
R
n
i
≤
0 and Rsi>μRni and Rsi≤0Figure 4
Interface states.At each iteration, the interface forceR 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, Ca and ΔF in (9) should be modified until all the conditions are satisfied. The finial interface force R can be obtained after all the steps are performed. Using (1), the displacement U 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 radiusR1 and R2, 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 b=4FR1R2/πLE′R1+R2, 1/E′=1-v12/E1+1-v22/E2, the maximum contact pressure po=2F/πbL, and the approach distance d of theoretical solution refers to the help file of ANSYS. Here E1, E2, v1, and v2 are Young’s modulus and Poisson’s ratio for cylinder, respectively.
The material parameters of Example1 are E1=30 GPa, v1=0.25, R1=10 mm, E2=29.12 GPa, v2=0.3, R2=13 mm, and loading F=3200 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 areb=1.1609 mm, po=1619.9 N/mm2, and d=0.4202 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 Example1.
Calculation parameters
Hertz solution
ANSYS solution
Present
P
0 (N/mm2)
1698.1
1617.7
1619.9
b (mm)
1.1997
1.1609
1.1069
d (mm)
0.4181
0.4213
0.4202Figure 5
Example1: contact between two cylinders.Figure 6
Example1: 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 forceP=1500 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 E=29.12 GPa, Poisson’s ratio v=0.3, cohesion c=0, and the coefficient of friction μ=0.2.
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−3mm. 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 Section4 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
Example2: frictional contact between cantilever beams (unit: m).Figure 8
Example2: 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 Figure9, and the structural planes are assumed to have no thickness. The material parameters of rock and structural planes are Young’s modulus E=15 GPa, Poisson’s ratio v=0.25, density ρ=2400kg/m3, coefficient of friction μ=0.32, cohesive strength c=0.04 MPa, and initial tensile strength σt=0. 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. Figures11 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.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).
(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.



---

*Source: 907620-2015-05-19.xml*