Multibody Finite Element Method and Application in Hydraulic Structure Analysis

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 proposedmethod successfully. Furthermore, this method is adopted to study the stability issues of underground houses of a large hydropower station.


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][7][8], dams [9,10], dynamic block [11][12][13], and others.Although many modifications and improvements have been performed to the original DDA [14][15][16][17][18][19][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.

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

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 Ω 1 and Ω 2 , and Ω 2 is assumed as a rigid foundation.Applying (2) Interface node to the contact model as suggested in Figure 1, then equations can be obtained: 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 where   =  −1   is the rigid displacement matrix of Ω 1 and   is the additional displacement constraints imposed on Ω 1 .The relative displacement between Ω 1 and Ω 2 is defined by According to Newton's third law, it is obvious that  1 = − 2 = .Combining ( 5), (6), and ( 7), we have where It is assumed that the initial state of interface is continuous; that is,  * = 0.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.
When Ω 1 is completely free, the displacement   , V  of node (  ,   ) and the displacement V  of node (  ,   ) are known in advance.Let {  } = {   V  V  }  ; the displacement of any node (  ,   ) can be expressed as is expressed as When Ω 1 has a translational degree of freedom, the displacement of freedom degree   is V  .Let {  } = {V  }; we have When Ω 1 has a rotational freedom around the center (  ,   ), the displacement of node (  ,   ) perpendicular to the radius   is V  .Let {  } = {V  }; the displacement of any node (  ,   ) is then given by

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 Ω 1 .
When Ω 1 is completely free, we obtain In the case, Ω 1 has a translational degree of freedom; we can get In the case, Ω 1 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 with where C , R , 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 as Here the subscript  denotes the number of the node pairs.Equation ( 23) is also applicable to the problems which contain rigid displacement.

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.

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  1 and  2 , 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 2 )/ 2 , the maximum contact pressure   = 2/, and the approach distance  of theoretical solution refers to the help file of ANSYS.Here  1 ,  2 , V 1 , and V 2 are Young's modulus and Poisson's ratio for cylinder, respectively.
The material parameters of Example 1 are  1 = 30 GPa, V 1 = 0.25,  1 = 10 mm,  2 = 29.12GPa, V 2 = 0.3,  2 = 13 mm, and loading  = 3200 N/mm 2 .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  = 1.1609 mm,   = 1619.9N/mm 2 , and  = 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.
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  = 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  = 29.12GPa, Poisson's ratio V = 0.3, cohesion  = 0, and the coefficient of friction

Fixed state
Free state Slip state  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.

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

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-topoint contact model is used in the proposed method.Further studies and developments are going to be undertaken for multibody FEM.

Figure 1 :
Figure 1: Schematic of two blocks initially in contact.

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

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

Table 2 :
Results comparison of Example 1.