Application of Base Force Element Method on Complementary Energy Principle to Rock Mechanics Problems

The four-mid-node plane model of base force element method (BFEM) on complementary energy principle is used to analyze the rock mechanics problems. The method to simulate the crack propagation using the BFEM is proposed. And the calculation method of safety factor for rock mass stability was presented for the BFEM on complementary energy principle. The numerical researches show that the results of the BFEM are consistent with the results of conventional quadrilateral isoparametric element and quadrilateral reduced integration element, and the nonlinear BFEM has some advantages in dealing crack propagation and calculating safety factor of stability.

In recent years, some scholars are studying other types of numerical analysis methods, such as boundary element method [44,45] and meshless method [46,47].And some scholars still adhere to explore the finite element method based on complementary energy principle [48][49][50][51].However, these methods have not been widely applied in engineering.
In this paper, the base force element method (BFEM) on complementary energy principle is used to analyze the engineering problems of rock mechanics.The "base forces" was introduced by Gao [52], who used the concept to replace various stress tensors for the description of the stress state at a point.These base forces can be directly obtained from the strain energy.For large deformation problems, when the base forces were adopted, the derivation of basic formulae was simplified by Gao [53] and Gao et al. [54][55][56].Based on the concept of the base forces, precise expressions for stiffness and compliance matrices for the FEM were obtained by Gao [52].The applications of the stiffness matrix to the plane problems of elasticity using the plane quadrilateral element and the polygonal element were researched by Peng et al. [37].Using the concept of base forces as state variables, a three-dimensional formulation of base force element method (BFEM) on complementary energy principle was proposed by Peng and Liu [35] for geometrically nonlinear problems.And the new finite element method based on the concept of base forces was called as the Base Force Element Method (BFEM) by Peng and Liu [35].A three-dimensional model of base force element method (BFEM) on complementary energy principle was proposed by Liu and Peng [36] for elasticity problems.A 4-mid-node plane element model of the BFEM on complementary energy principle was proposed by Peng et al. [38] for geometrically nonlinear problem, which is derived by assuming that the stress is uniformly distributed on each edges of a plane element.In the paper [39], an arbitrary convex polygonal element model of the BFEM on complementary energy principle was proposed for geometrically nonlinear problem.In the paper [43], a 4mid-node plane model of BFEM on complementary energy principle was researched, and its computational performance was studied.The convex polygonal element model of BFEM on complementary energy principle was given by Peng et al. [40] for arbitrary mesh problems.In the paper [41], the concave polygonal element model of BFEM on complementary energy principle was proposed for the concave polygonal mesh problems.In the paper [42], the BFEM on potential energy principle was used to analyze recycled aggregate concrete (RAC) on mesolevel, in which the model of BFEM with triangular element was derived, and the simulation results of the BFEM agree with the test results of recycled aggregate concrete.In recently, the BFEM on damage mechanics has been used to analyze the compressive strength, the size effects of compressive strength, and fracture process of concrete at mesolevel, and the analysis method is the new way for investigating fracture mechanism and numerical simulation of mechanical properties for concrete.
The purpose of this paper is to survey the base forces element method on complementary energy principle for largescale computing problems in rock engineering problems.

Model of the BFEM
2.1.Compliance Matrix.Consider a 4-mid-node plane element as shown in Figure 1; the compliance matrix of a base force element can be obtained as [43] in which  is Young's modulus, ] is Poisson's ratio,  is the area of an element, U is the unit tensor, and   is the dot product of radius vectors r  and r  at points  and .
For a plane rectangular coordinate system, the radius vectors r  and r  of points  and  can be written as in which e  , e  are the unit vectors.Further, the compliance matrix of an element can be reduced as follows: For a plane strain problem, it is necessary to replace  by /(1 − ] 2 ) and ] by ]/(1 − ]) in ( 1) and (3).
The characteristics of the BFEM on complementary energy principle are that the model does not introduce an interpolating function and is not necessary to introduce the Gauss integral for calculating the compliance coefficient at a point.

Governing Equations.
The total complementary energy of the elastic system which has  elements can be written as where    is the complementary energy of an element and T  and u  are the resultant force vectors and the given displacement acting on the center node  of the edge , respectively.
The equilibrium conditions can be released by the Lagrange multiplier method, and a new complementary energy function for an element can be introduced as follows: where  =   e  +   e  and   are the Lagrange multipliers.For the elastic system, the modified total complementary energy function of the elastic system which contains  elements should meet the following equation by means of the modified complementary energy principle: Further, (6) can be expressed as The first of ( 7) is the compatibility equations and displacement boundary conditions for the elastic system.According to this equation, the displacement boundary conditions in this paper can be implemented in the BFEM.The second of (7) is the force equilibrium equation of each element.The third of (7) is the moment equilibrium equation of each element.These are the governing equations of the BFEM.From the equations, we can obtain the resultant forces acting on the center points of the edges of all elements.

Stress
Tensor of an Element.Consider the 4-mid-node plane element as shown in Figure 1; the real stress  of the element can be replaced by the average stress  if the element is small enough.According to the definitions of Cauchy stress tensors, the stress expressions of an element can be obtained as where ⊗ is the dyadic symbol, T  and r  are the resultant force vectors acting on the center node  of the edge  and the radius vector of the node , respectively, and the summation rule is implied.
For a plane rectangular coordinate system, the force vectors T  of the node  can be written as where    and    are the components of the force vector T  along coordinates  and , respectively.
Further, the stress tensors of an element can be reduced as follows: 2.4.Displacement Vector of Nodes.According to the governing equation of an element of the BFEM, the explicit expression of displacement can be obtained as in which  is the alternating tensor,  and   are the Lagrange multipliers [43], and  and  can be expressed as Further, the displacement vectors of an element can be reduced as follows [43]: Figure 2: Equivalent node loads of gravity in an element.

Simulation of Gravity and Material
In engineering problems, regardless of the dam, rock slope, or other structure of rock mass, the gravity of structure should be considered in the numerical calculation.We did not take into account the gravity problem when we previously prepared the program of BFEM on complementary energy principle.In order to consider the gravity of structure, three problems must be solved, including the calculation problem of equivalent node loads in the BFEM on complementary energy principle, the problem of exerting gravity in the software of the BFEM and the calculation problem of stress tensor in an element when the gravity is added.

Equivalent Node Loads of Gravity.
For the BFEM on complementary energy principle, the equivalent node loads of gravity in an element can be calculated according to the principle of virtual work, as shown in Figure 2, and the expression can be given as or where  is the thickness of an element,  is the density of material, and  is acceleration of gravity.

Stress Calculation of an Element.
When the gravity of an element is not taken into account, as shown in Figure 3, we calculate the force acting on the edges of an element first.
Then, the stress tensor of the element can be calculated by (8) or (10).Further, the stress tensors of an element can also be reduced as follows: where T  =    e  +    e  , r  =   e  +   e  , ( = 1, 2, 3, 4),   and   are the coordinates of node , respectively.When the gravity of an element is taken into account and there is a free boundary, as shown in Figure 4, the stress tensor of the element can be calculated as When the gravity of an element is taken into account and there is a force boundary condition, as shown in Figure 5, the stress tensor of the element can be calculated as in which F =   e  +   e  .
When the gravity of an element is not taken into account and there is a force boundary condition, as shown in Figure 6, the stress tensor of the element can be calculated as

Simulation of Different Materials.
We adopt two onedimensional array variables to reflect the change of elastic modulus and Poisson's ratio of different materials, respectively.

The Nonlinear Model of BFEM for Crack Propagation Problems
The conventional displacement model of FEM requires the mesh reconstruction for the crack propagation problems.Therefore, the conventional FEM has deficiencies for the crack propagation problems.Because the contact forces between each element are used, the BFEM on complementary energy principle has advantage in the simulation of crack propagation problems.In the BFEM on complementary energy principle, we only need to deal with the constraint conditions and compatibility equations of displacement.

Failure Criteria of the Contact Interface between Two
Elements.According to the control equations of the BFEM on complementary energy principle, the forces acting on the edge on an element can easily be calculated.According to the failure criteria expressed by the forces acting on the edge of an element, we can judge whether the interface of elements is cracking.If there is cracking of element interface, the nonlinear processing must be carried out.Here,   and   are the normal interface force and tangential interface force at contact interface between two elements, respectively. and  are the cohesion and the internal friction coefficient of the contact interface, respectively. is the length of an element. 0 is the tensile strength of an element, and it is positive in tension.

Condition of Positive Slip at Contact
Interface.When   < 0 and   ≥ − ⋅   +  ⋅ , the contact interface of the two elements began to slip.We need to get rid of the tangential displacement constraint at the contact interface between the two elements that it has been cracked and put the frictional force as load that is used to an initial force condition.The new load acting on the contact interface between the two elements in the second loop of calculation can be written as follows: When 0 <   <  0 and   ≥  ⋅ , the contact interface of the two elements begins to slip.We need to get rid of the tangential displacement constraint and the normal displacement constraint at the contact interface of the two elements.And the contact interface will crack after slipping.

Condition of Negative Slip at Contact Interface. When
< 0 and   ≤  ⋅   −  ⋅ , the contact interface of the two elements began to slip.We need to get rid of the tangential displacement constraint at the contact interface between the two elements that it has been cracked, and put the frictional force as load that is an initial force condition.The new load acting on the contact interface between the two elements in the second loop of calculation is When 0 <   <  0 and   ≤ − ⋅ , the contact interface of the two elements begins to slip.We need to get rid of the tangential displacement constraint and the normal displacement constraint at the contact interfaces of elements.And the contact interface will crack after slipping.

Condition of Pull Cracking at Contact Interface.
When   ≥  0 , the contact interface of the two elements begins to crack.We need to get rid of the tangential displacement constraint and the normal displacement constraint at the contact interfaces of elements.
After the above checks, the computer program of BFEM uses the new loads and constraint conditions to solve the governing equations of BFEM on complementary energy principle and obtains the new forces acting on the contact interfaces of elements.The program repeats the above checks until there are no new cracking at the contact interfaces or the solutions of nonlinear equations cannot be convergence since the interface cracks are too long.

Flow Chart of the Nonlinear BEFM for Crack Propagation
Problems.The flow chart of the nonlinear BEFM of crack propagation problems can be shown in Figure 7.

Calculation Method on Safety Factor of Stability in the BFEM
There are many methods to calculate the safety factor in engineering.The traditional finite element method usually used the rock joint elements to calculate factor of safety along the joint path.When the base force element method is used to analyze the stability of the rock mass in order to get the safety factor along the joint path, it is very easy.First, we calculate the surface forces of all elements according to the different load combinations.Then, we accumulate the sliding resistances and the sliding forces along the sliding path, respectively.Further, the safety factor of stability along the sliding path can obtained by the following equation: where   and   are the normal interface force and tangential interface force at contact interface between two elements along the sliding path, respectively.  and   are the cohesion and the internal friction coefficient of the contact interface between two elements along the sliding path, respectively.  is the length of the interface in the element along the sliding path.

Example 1: A Rock Pillar under the Action of Gravity.
Consider a thick pillar of rock under the action of gravity shown in Figure 8.And its width is 5 m, its height is 10 m, modulus of elasticity  = 1 × 10 8 Pa, Poisson ratio ] = 0.3, density of rock  = 2.45 t/m 3 , and acceleration of gravity  = 9.8 m/s 2 .The calculation is considered into the plane stress problem.The calculation is done using three different element meshes with the center nodes of edges of elements as shown in Figure 9.
The values of stress components and displacement components of the rock pillar are listed in Tables 1-4, respectively.Comparisons of the results from the conventional quadrilateral isoparametric element (Q4 model) and quadrilateral reduced integration element (Q4R model) are also given in Tables 1-4, respectively.The numerical results of the present     10.And its modulus of elasticity  1 = 10 9 ,  2 = 10 8 ,  3 = 10 7 , and  4 = 10 6 , respectively.And Poisson ratio ] = 0.3, the shear stress on surface of the structure  = 1.The calculation is considered into the plane stress problem and the dimensionless values.The calculation is done using the mesh with the center nodes of edges of elements as shown in Figure 11.
The values of stress components of the rock pillar are listed in Table 5, respectively.Comparisons of the results from the theoretical analysis are also given in Table 5, respectively.The numerical results of the present model are consistent with those of the theoretical analysis.

Example 3: A Rock Pillar under Water Pressure and
Gravity.Consider a rock pillar under water pressure and gravity shown in Figure 12.And its modulus of elasticity  = 10 8 Pa, Poisson' ratio ] = 0.3, density of rock  1 = 2.4 t/m 3 , density of water  2 = 1.0 t/m 3 , and acceleration of gravity  = 9.8 m/s 2 .The calculation is considered into the plane strain problem.The calculation is done using four different element meshes with the center nodes of edges of elements as shown in Figure 13.
The values of stress components and displacement components of the rock pillar are listed in Tables 6-12, respectively.Comparisons of the results from the conventional quadrilateral isoparametric element (Q4 model) and quadrilateral reduced integration element (Q4R model) are also given in Tables 6-12, respectively.The numerical results of the present model are consistent with those of Q4R model,  and the 4-mid-node element of BFEM has given good performance compared with Q4 model for the large aspect ratio of elements.Due to the different method calculating the equivalent node loads of water pressure between the base force elements and the element of traditional FEM, there is a slight error about the calculation results of the horizontal stresses   of the element in lower right corner of rock pillar between the BFEM and the Q4R model, as shown in Table 10.The calculation is done using the mesh with the center nodes of edges of elements as shown in Figure 15.In this calculation, we do not consider the initial geostress field.The boundaries of the foundation are used the fixed constraint.The origin of coordinates is located at the bottom of the dam, and is 25 meters away from the dam heel.The calculation is done using the mesh with the center nodes of edges of elements as shown in Figure 27.
In order to check whether the interface between the rock block and the ground will crack, we assume the friction coefficient of interface  = 0.5 and change the value of interface cohesion .The results of calculation using the computer program of the nonlinear BFEM shown in Section 4 and the failure criteria in Section 4.1 are as follows.
(6) When  = 1.0, the cracks are too long, and too little structural constraints have been insufficient to solve the equations.
When the interface cohesion  is 1.5, the case of crack propagation of rock block is shown in Figure 28, and the safety factor of stability is  = 2 which is consistent with the result of the theoretical analysis using the rigid limit equilibrium method.The calculation is done using the mesh with the center nodes of edges of elements as shown in Figure 29.(1) When  = 10 6 Pa and  = 0.95, there is no cracks.
The case of crack propagation of dam is shown in Figure 30.When  = 10 6 Pa and  = 0.5, the safety factor of stability  = 4.0 which is consistent with the result of the theoretical analysis.When  = 0.5 × 10 6 Pa and  = 0.5,

Conclusions
In this paper, the base force element method (BFEM) on complementary energy principle is used to analyze the rock mechanics problems.The methods to simulate the gravity of an element, the crack propagation and the safety factor of stability are proposed for the BFEM on complementary energy principle.The following conclusions can be drawn.
(1) The calculation results of the BFEM on complementary energy principle show that the numerical results of the present method coincide with the theoretical solution, the results of conventional quadrilateral isoparametric element (Q4 model) and quadrilateral reduced integration element (Q4R model).The correctness of the present method and its computer program is verified.(2) The research results show that the BFEM on complementary energy principle has a good computational precision and stability is not sensitive to the effects on the aspect ratio of element and can be used for largescale scientific and engineering computing.(3) The results of the BFEM for crack propagation problems show that the nonlinear BFEM can solve the cracking problem and simulate the crack propagation of the interface in rock mechanics engineering.(4) The BFEM on complementary energy principle was applied to analyze the stability of rock mass and dam, and the results of safety factor are consistent with the results of the theoretical solutions using the rigid limit equilibrium method.The research results show the present method can be easily used to calculate the safety factor in rock engineering.(5) This paper researched only the cracking problems with horizontal crack in rock mass and calculated only the safety factor of a single slip channel in rock mass.(6) The cracking problems of inclined cracks and the safety factor of multiple sliding channels in rock mass are studying, and the further research results will be published in the future.

Figure 4 :
Figure 4: Considering gravity and a free boundary.

Figure 5 :
Figure 5: Considering gravity and a force boundary condition.

Figure 6 :
Figure 6: Considering a force boundary condition and no gravity.

Figure 9 :
Figure 9: Three kinds of meshes for a rock pillar.

Figure 10 :
Figure 10: A rock pillar with four kinds of materials.

Figure 12 :
Figure 12: A rock pillar subjected to water pressure and gravity.

Figure 13 :
Figure 13: Four kinds of meshes for a rock pillar.

Figure 21 :Figure 22 :
Figure 21: The ℎ-  curves at the half height of the dam.

Figure 23 :Figure 24 :
Figure 23: The ℎ-  curves at the half height of the dam.

Figure 25 :
Figure 25: The ℎ-  curves at the half height of the dam.

Figure 26 :
Figure 26: A rock block subjected by the horizontal thrust and vertical pressure.

Figure 27 :
Figure 27: Meshes of the rock block.

Crack
(a) In the first cycle calculation Crack (b) In the second cycle calculation Crack (c) In the third cycle calculation Crack (d) In the fourth cycle calculation Crack (e) In the fifth cycle calculation

Figure 28 :
Figure 28: Crack propagation path of the rock block.

Figure 31 :
Figure 31: A gravity dam with a crack at the dam heel.

Table 1 :
Displacement   at the top of the pillar.

Table 2 :
Displacement   at ℎ = 5 m of the pillar.

Table 3 :
Stress   at center of the elements at top of the pillar.

Table 4 :
Stress   at center of the elements at bottom of the pillar.
model are consistent with those of the Q4 model and Q4R model and have shown good computational stability.

Table 5 :
Stress solution of the pillar with four kinds of materials under uniform shearing forces.

Table 6 :
Displacement   at  = 2.5 m of the pillar top.

Table 7 :
Displacement   at ℎ = 5 m on the right of the pillar.

Table 8 :
Displacement   at  = 2.5 m of the pillar top.

Table 9 :
Displacement   at ℎ = 5 m on the right of the pillar.

Table 10 :
Stress   at center of the element at lower right of pillar.

Table 11 :
Stress   at center of the element at lower right of pillar.

Table 12 :
Stress   at center of the element at lower right of pillar.Consider a concrete gravity dam shown in Figure14.And height of the dam is 65 m, bottom width is 49 m, the water level is 60 m, the elastic modulus of concrete  1 = 15 GPa, Poisson ratio of concrete ] 1 = 0.2, the elastic modulus of rock  2 = 30GPa, Poisson ratio of rock ] 2 = 0.3, density of concrete is 2.45 t/m 3 , density of water is 1 t/m 3 , and acceleration of gravity  = 9.8 m/s 2 .The calculation is considered into the plane strain problem and considered the effect of rock foundation of the dam.