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

Research Article

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

Yijiang Peng, Qing Guo, Zhaofeng Zhang, and Yanyan Shan

The Key Laboratory of Urban Security and Disaster Engineering, Ministry of Education, Beijing University of Technology, Beijing 100124, China

Received 24 August 2014; Accepted 5 September 2014

Academic Editor: Song Cen

Copyright © 2015 Yijiang Peng 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.

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.

1. Introduction

The finite element method (FEM) has been playing a very important role in solving various problems in engineering and science. However, the conventional finite element method (FEM) based on the displacement model has some shortcomings, such as large deformation, treatment of incompressible materials, bending of thin plates, and moving boundary problems. In the past decades, numerous efforts techniques have been proposed for developing finite element models which are robust and insensitive to mesh distortion, such as the hybrid stress method [1–4], the equilibrium models [5, 6], the mixed approach [7], the integrated force method [8–11], the incompatible displacement modes [12, 13], the assumed strain method [14–17], the enhanced strain modes [18, 19], the selectively reduced integration scheme [20], the quasiconforming element method [21], the generalized conforming method [22], the Alpha finite element method [23], the new spline finite element method [24, 25], the unsymmetric method [26–29], the new natural coordinate methods [30–33], the smoothed finite element method [34], and the base force element method [35–43].

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–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–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 4-mid-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 large-scale computing problems in rock engineering problems.

2. 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, is the unit tensor, and is the dot product of radius vectors and at points and .

Figure 1: Four-mid-node plane element.

For a plane rectangular coordinate system, the radius vectors and of points and can be written as

in which ,   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 and by 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.

2.2. 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 and 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 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.

2.3. 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, and 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 of the node can be written as

where and are the components of the force vector 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]:

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

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

Figure 2: Equivalent node loads of gravity in an element.

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

Figure 3: An element without gravity.

Further, the stress tensors of an element can also be reduced as follows:

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

Figure 4: Considering gravity and a free boundary.

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 .

Figure 5: Considering gravity and a force boundary condition.

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

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

3.3. Simulation of Different Materials

We adopt two one-dimensional array variables to reflect the change of elastic modulus and Poisson’s ratio of different materials, respectively.

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

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

4.1.1. Condition of Elastic State at Contact Interface

When and or and , the contact interface is not cracked and in the elastic state. 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. is the tensile strength of an element, and it is positive in tension.

4.1.2. Condition of Positive Slip at Contact Interface

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

4.1.3. Condition of Negative Slip at Contact Interface

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

4.1.4. Condition of Pull Cracking at Contact Interface

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

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

Figure 7: Flow chart of the nonlinear BEFM for crack propagation problems.

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

6. Numerical Examples

6.1. 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 , Poisson ratio , density of rock , and acceleration of gravity  m/s². The calculation is considered into the plane stress problem.

Figure 8: A rock pillar under the action of gravity.

The calculation is done using three different element meshes with the center nodes of edges of elements as shown in Figure 9.

(a) 100 elements, 220 nodes

(b) 200 elements, 430 nodes

(c) 300 elements, 640 nodes

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

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 model are consistent with those of the Q4 model and Q4R model and have shown good computational stability.

Table 1: Displacement at the top of the pillar.

Meshes
BFEM model (×10⁻² m)
Q4 model (×10⁻² m)
Q4R model (×10⁻² m)

10 × 10
−1.1939
−1.1917
−1.2048

10 × 20
−1.1940
−1.1931
−1.1964

10 × 30
−1.1940
−1.1934
−1.1940

Table 2: Displacement at  m of the pillar.

Meshes
BFEM (×10⁻³ m)
Q4 model (×10⁻³ m)
Q4R model (×10⁻³ m)

10 × 10
−8.7592
−8.6975
−8.8045

10 × 20
−8.7630
−8.7118
−8.6726

10 × 30
−8.7637
−8.7149
−8.7363

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

Meshes
BFEM (kPa)
Q4 model (kPa)
Q4R model (kPa)

10 × 10
−12.016
−12.244
−12.062

10 × 20
−6.005
−6.0773
−6.013

10 × 30
−4.003
−4.0382
−4.002

Table 4: Stress at center of the elements at bottom of the pillar.

Meshes
BFEM (kPa)
Q4 model (kPa)
Q4R model (kPa)

10 × 10
−244.194
−242.945
−244.349

10 × 20
−262.976
−260.116
−263.387

10 × 30
−272.102
−268.002
−272.552

6.2. Example  2: Analysis for a Rock Pillar with Four Materials under Pure Shear Effect

Consider a rock pillar with four materials under pure shear effect shown in Figure 10. And its modulus of elasticity , , , and , respectively. And Poisson ratio , the shear stress on surface of the structure . The calculation is considered into the plane stress problem and the dimensionless values.

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

The calculation is done using the mesh with the center nodes of edges of elements as shown in Figure 11.

Figure 11: Mesh with 200 elements, 430 nodes.

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.

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

BFEM
0.000
1.000
0.000

Exact
0.000
1.000
0.000

6.3. 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 , Poisson’ ratio , density of rock , density of water , and acceleration of gravity  m/s². The calculation is considered into the plane strain problem.

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

The calculation is done using four different element meshes with the center nodes of edges of elements as shown in Figure 13.

(a) 800 elements, 1660 nodes

(b) 400 elements, 850 nodes

(c) 320 elements, 688 nodes

(d) 240 elements, 526 nodes

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

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.

Table 6: Displacement at  m of the pillar top.

Meshes
BFEM (×10⁻² m)
Q4 model (×10⁻² m)
Q4R model (×10⁻² m)

20 × 40
4.0219
4.0056
4.0190

10 × 40
4.0415
4.0094
4.0398

8 × 40
4.0563
4.0135
4.0532

6 × 40
4.0891
4.0239
4.0863

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

Meshes
BFEM (×10⁻² m)
Q4 model (×10⁻² m)
Q4R model (×10⁻² m)

20 × 40
2.0837
2.0758
2.0830

10 × 40
2.0816
2.0667
2.0815

8 × 40
2.0820
2.0626
2.0823

6 × 40
2.0849
2.0566
2.0862

Table 8: Displacement at  m of the pillar top.

Meshes
BFEM (×10⁻³ m)
Q4 model (×10⁻³ m)
Q4R model (×10⁻³ m)

20 × 40
−9.8891
−9.8869
−9.8895

10 × 40
−9.8888
−9.8856
−9.8898

8 × 40
−9.8889
−9.8850
−9.8943

6 × 40
−9.8889
−9.8840
−9.8926

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

Meshes
BFEM model (×10⁻² m)
Q4 model (×10⁻² m)
Q4R model (×10⁻² m)

20 × 40
−1.5918
−1.5878
−1.5913

10 × 40
−1.5488
−1.5412
−1.5483

8 × 40
−1.5290
−1.51894
−1.5284

6 × 40
−1.4983
−1.4832
−1.4975

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

Meshes
BFEM (kPa)
Q4 model (kPa)
Q4R model (kPa)

20 × 40
−128.963
−179.706
−132.974

10 × 40
−142.729
−191.962
−150.839

8 × 40
−141.996
−193.919
−152.815

6 × 40
−134.954
−194.249
−150.941

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

Meshes
BFEM (kPa)
Q4 model (kPa)
Q4R model (kPa)

20 × 40
152.442
172.480
153.366

10 × 40
153.439
167.708
154.637

8 × 40
151.262
163.966
152.612

6 × 40
146.664
157.624
148.288

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

Meshes
BFEM (kPa)
Q4 model (kPa)
Q4R model (kPa)

20 × 40
−802.328
−796.173
−802.741

10 × 40
−707.932
−700.016
−708.582

8 × 40
−675.867
−666.975
−676.616

6 × 40
−634.179
−623.753
−635.039

6.4. Example  4: Stress Analysis of Concrete Gravity Dam

Consider a concrete gravity dam shown in Figure 14. And height of the dam is 65 m, bottom width is 49 m, the water level is 60 m, the elastic modulus of concrete , Poisson ratio of concrete , the elastic modulus of rock , Poisson ratio of rock , density of concrete is 2.45 t/m³, density of water is 1 t/m³, and acceleration of gravity  m/s². The calculation is considered into the plane strain problem and considered the effect of rock foundation of the dam.

Figure 14: A concrete gravity dam.

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.

Figure 15: Meshes of the dam and its foundation (1344 elements, 2809 nodes).

The values of stress components and displacement components of the dam are plotted in Figures 16–25, respectively. Comparisons of the results from the conventional quadrilateral isoparametric element (Q4 model) and quadrilateral reduced integration element (Q4R model) are also given in Figures 16–25, respectively. The numerical results of the present model are consistent with those of the Q4 model and Q4R model and have shown good computational stability.

Figure 16: The curves at the upstream face of the dam.

Figure 17: The curves at the upstream face of the dam.

Figure 18: The curves at the upstream face of the dam.

Figure 19: The curves at the upstream face of the dam.

Figure 20: The curves at the upstream face of the dam.

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

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

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

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

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

6.5. Example  5: Simulation and Analysis on the Horizontal Crack Propagation of Rock Block

Consider a rock block subjected by the horizontal thrust and vertical pressure shown in Figure 26. For the convenience of study, we do not consider the weight. And we use the dimensionless numerical analysis. Assuming elastic modulus , Poisson ratio , tensile strength of a large number, and the uniform load and . The calculation is considered into the plane stress problem and the dimensionless values.

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

The calculation is done using the mesh with the center nodes of edges of elements as shown in Figure 27.

Figure 27: Meshes of the rock block.

In order to check whether the interface between the rock block and the ground will crack, we assume the friction coefficient of interface 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.

(1)When , there is no cracks.(2)When , there is one element interface crack.(3)When , there are three elements’ interfaces cracks.(4)When , there are five elements’ interfaces cracks.(5)When , there are seven elements’ interfaces cracks.(6)When , 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 which is consistent with the result of the theoretical analysis using the rigid limit equilibrium method.

(a) In the first cycle calculation

(b) In the second cycle calculation

(c) In the third cycle calculation

(d) In the fourth cycle calculation

(e) In the fifth cycle calculation

Figure 28: Crack propagation path of the rock block.

6.6. Example  6: Analysis on the Crack Propagation of Concrete Gravity Dam

Consider a concrete gravity dam shown in Figure 14. And height of the dam is 65 m, bottom width is 49 m, the water level is 60 m, The elastic modulus of concrete , Poisson ratio of concrete , density of concrete is 2.45 t/m³, density of water is 1 t/m³, and acceleration of gravity  m/s². The calculation is considered into the plane strain problem and is not considered the effect of rock foundation of the dam.

The calculation is done using the mesh with the center nodes of edges of elements as shown in Figure 29.

Figure 29: Mesh of gravity dam (800 elements, 1660 nodes).

6.6.1. No Initial Crack in Dam Foundation

We assume that the foundation surface of dam is weak structural interface and analyze the crack propagation and the safety factor by adjusting the value of the friction coefficient and the 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.

(1)When and , there is no cracks.(2)When and , there are two elements’ interface cracks.(3)When and , there are eight elements’ interfaces cracks.

The case of crack propagation of dam is shown in Figure 30. When and , the safety factor of stability which is consistent with the result of the theoretical analysis. When and , the safety factor of stability which is consistent with the results of the theoretical analysis using the rigid limit equilibrium method.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 30: Crack propagation path of the gravity dam.

6.6.2. Existing an Initial Crack in Dam Foundation

There is an initial crack in the dam heel as shown in Figure 31. We assume that the dam foundation surface is weak structural interface and analyze the crack propagation and the safety factor by adjusting the value of the friction coefficient and the 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.

(1)When and , there is no cracks.(2)When and , there is one element interface cracks.(3)When and , there are three element’ interface cracks.(4)When and , there are seven elements’ interfaces cracks.

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

7. 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 large-scale 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.

Conflict of Interests

The authors declare that there is no conflict of interests.

Acknowledgments

This work is supported by the National Natural Science Foundation of China, nos. 10972015 and 11172015 and the preexploration project of the Key Laboratory of Urban Security and Disaster Engineering, Ministry of Education, Beijing University of Technology, no. USDE201404.

References

T. H. H. Pian, “Derivation of element stiffness matrices by assumed stress distributions,” AIAA Journal, vol. 2, no. 7, pp. 1333–1336, 1964.

T. H. Pian and D. P. Chen, “Alternative ways for formulation of hybrid stress elements,” International Journal for Numerical Methods in Engineering, vol. 18, no. 11, pp. 1679–1684, 1982.

T. H. H. Pian and K. Sumihara, “Rational approach for assumed stress finite elements,” International Journal for Numerical Methods in Engineering, vol. 20, no. 9, pp. 1685–1695, 1984.

C. Zhang, D. Wang, J. Zhang, W. Feng, and Q. Huang, “On the equivalence of various hybrid finite elements and a new orthogonalization method for explicit element stiffness formulation,” Finite Elements in Analysis and Design, vol. 43, no. 4, pp. 321–332, 2007.

B. Fraeijs de Veubeke, “Displacement and equilibrium models in the finite element method,” in Stress Analysis, O. C. Zienkiewicz and G. S. Holister, Eds., pp. 145–197, John Wiley & Sons, New York, NY, USA, 1965.

B. F. de Veubeke, “A new variational principle for finite elastic displacements,” International Journal of Engineering Science, vol. 10, no. 9, pp. 745–763, 1972.

R. L. Taylor and O. C. Zienkiewicz, “Complementary energy with penalty functions in finite element analysis,” in Energy Methods in Finite Element Analysis, R. Glowinski, Ed., pp. 153–174, John Wiley & Sons, New York, NY, USA, 1979.

S. N. Patniak, “An integrated force method for discrete analysis,” International Journal for Numerical Methods in Engineering, vol. 6, no. 2, pp. 237–251, 1973.

S. N. Patnaik, “The integrated force method versus the standard force method,” Computers and Structures, vol. 22, no. 2, pp. 151–163, 1986.

S. N. Patnaik, “The variational energy formulation for the integrated force method,” AIAA Journal, vol. 24, no. 1, pp. 129–137, 1986.

S. N. Patnaik, L. Berke, and R. H. Gallagher, “Integrated force method versus displacement method for finite element analysis,” Computers and Structures, vol. 38, no. 4, pp. 377–407, 1991.

E. L. Wilson, R. L. Tayler, W. P. Doherty, and J. Ghaboussi, “Incompatible displacement models,” in Numerical and Computational Methods in Structural Mechanics, S. J. Fenves, N. Perrone, A. R. Robinson, and W. C. Schnobrich, Eds., pp. 43–57, Academic Press, New York, NY, USA, 1973.

R. L. Taylor, P. J. Beresford, and E. L. Wilson, “A non-conforming element for stress analysis,” International Journal for Numerical Methods in Engineering, vol. 10, no. 6, pp. 1211–1219, 1976.

J. C. Simo and T. J. R. Hughes, “On the variational foundations of assumed strain methods,” Journal of Applied Mechanics, vol. 53, no. 1, pp. 51–54, 1986.

J. C. Simo and M. S. Rifai, “Class of mixed assumed strain methods and the method of incompatible modes,” International Journal for Numerical Methods in Engineering, vol. 29, no. 8, pp. 1595–1638, 1990.

R. H. Macneal, “Derivation of element stiffness matrices by assumed strain distributions,” Nuclear Engineering and Design, vol. 70, no. 1, pp. 3–12, 1982.

T. Belytschko and L. P. Bindeman, “Assumed strain stabilization of the 4-node quadrilateral with 1-point quadrature for nonlinear problems,” Computer Methods in Applied Mechanics and Engineering, vol. 88, no. 3, pp. 311–340, 1991.

R. Piltner and R. L. Taylor, “A quadrilateral mixed finite element with two enhanced strain modes,” International Journal for Numerical Methods in Engineering, vol. 38, no. 11, pp. 1783–1808, 1995.

R. Piltner and R. L. Taylor, “A systematic constructions of B-bar functions for linear and nonlinear mixed-enhanced finite elements for plane elasticity problems,” International Journal for Numerical Methods in Engineering, vol. 44, no. 5, pp. 615–639, 1997.

T. J. R. Hughes, “Generalization of selective integration procedures to anisotropic and nonlinear media,” International Journal for Numerical Methods in Engineering, vol. 15, no. 9, pp. 1413–1418, 1980.

T. Limin, C. Wanji, and L. Yingxi, “Formulation of quasi-conforming element and Hu-Washizu principle,” Computers and Structures, vol. 19, no. 1-2, pp. 247–250, 1984.

L. Yu-qiu and H. Min-feng, “A generalized conforming isoparametric element,” Applied Mathematics and Mechanics, vol. 9, no. 10, pp. 929–936, 1988.

G. R. Liu, T. Nguyen-Thoi, and K. Y. Lam, “A novel alpha finite element method (αFEM) for exact solution to mechanics problems using triangular and tetrahedral elements,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 45–48, pp. 3883–3897, 2008.

J. Chen, C.-J. Li, and W.-J. Chen, “A 17-node quadrilateral spline finite element using the triangular area coordinates,” Applied Mathematics and Mechanics (English Edition), vol. 31, no. 1, pp. 125–134, 2010.

J. Chen, C.-J. Li, and W.-J. Chen, “A family of spline finite elements,” Computers and Structures, vol. 88, no. 11-12, pp. 718–727, 2010.

S. Rajendran and K. M. Liew, “A novel unsymmetric 8-node plane element immune to mesh distortion under a quadratic displacement field,” International Journal for Numerical Methods in Engineering, vol. 58, no. 11, pp. 1713–1748, 2003.

S. Rajendran, “A technique to develop mesh-distortion immune finite elements,” Computer Methods in Applied Mechanics and Engineering, vol. 199, no. 17–20, pp. 1044–1063, 2010.

E. T. Ooi, S. Rajendran, and J. H. Yeo, “A 20-node hexahedron element with enhanced distortion tolerance,” International Journal for Numerical Methods in Engineering, vol. 60, no. 15, pp. 2501–2530, 2004.

E. T. Ooi, S. Rajendran, and J. H. Yeo, “Remedies to rotational frame dependence and interpolation failure of US-QUAD8 element,” Communications in Numerical Methods in Engineering, vol. 24, no. 11, pp. 1203–1217, 2008.

Y. Long, L. Juxuan, Z. Long, and C. Song, “Area co-ordinates used in quadrilateral elements,” Communications in Numerical Methods in Engineering, vol. 15, no. 8, pp. 533–543, 1999.

Y. Q. Long, S. Cen, and Z. F. Long, Advanced Finite Element Method in Structural Engineering, Springer/Tsinghua University Press, Berlin, Germany, 2009.

Z. F. Long, J. X. Li, S. Cen, and Y. Q. Long, “Some basic formulae for area coordinates used in quadrilateral elements,” Communications in Numerical Methods in Engineering, vol. 15, no. 12, pp. 841–852, 1999.

Z.-F. Long, S. Cen, L. Wang, X.-R. Fu, and Y.-Q. Long, “The third form of the quadrilateral area coordinate method (QACM-III): theory, application, and scheme of composite coordinate interpolation,” Finite Elements in Analysis and Design, vol. 46, no. 10, pp. 805–818, 2010.

G. R. Liu, K. Y. Dai, and T. T. Nguyen, “A smoothed finite element method for mechanics problems,” Computational Mechanics, vol. 39, no. 6, pp. 859–877, 2007.

Y. Peng and Y. Liu, “Base force element method of complementary energy principle for large rotation problems,” Acta Mechanica Sinica, vol. 25, no. 4, pp. 507–515, 2009.

Y. Liu and Y. Peng, “Base force element method (BFEM) on complementary energy principle for linear elasticity problem,” Science China: Physics, Mechanics and Astronomy, vol. 54, no. 11, pp. 2025–2032, 2011.

Y. Peng, Z. Dong, B. Peng, and Y. Liu, “Base force element method (BFEM) on potential energy principle for elasticity problems,” International Journal of Mechanics and Materials in Design, vol. 7, no. 3, pp. 245–251, 2011.

Y. Peng, Z. Dong, B. Peng, and N. Zong, “The application of 2D base force element method (BFEM) to geometrically non-linear analysis,” International Journal of Non-Linear Mechanics, vol. 47, no. 3, pp. 153–161, 2012.

Y.-J. Peng, J.-W. Pu, B. Peng, and L.-J. Zhang, “Two-dimensional model of base force element method (BFEM) on complementary energy principle for geometrically nonlinear problems,” Finite Elements in Analysis and Design, vol. 75, pp. 78–84, 2013.

Y. J. Peng, N. N. Zong, L. J. Zhang, and J. W. Pu, “Application of 2D base force element method with complementary energy principle for arbitrary meshes,” Engineering Computations, vol. 31, no. 4, pp. 1–15, 2014.

Y. Peng, L. Zhang, J. Pu, and Q. Guo, “A two-dimensional base force element method using concave polygonal mesh,” Engineering Analysis with Boundary Elements, vol. 42, pp. 45–50, 2014.

Y. Peng, Y. Liu, J. Pu, and L. Zhang, “Application of base force element method to mesomechanics analysis for recycled aggregate concrete,” Mathematical Problems in Engineering, vol. 2013, Article ID 292801, 8 pages, 2013.

Y. Liu, Y. Peng, L. Zhang, and Q. Guo, “A 4-mid-node plane model of base force element method on complementary energy principle,” Mathematical Problems in Engineering, vol. 2013, Article ID 706759, 8 pages, 2013.

C. Y. Dong and G. L. Zhang, “Boundary element analysis of three dimensional nanoscale inhomogeneities,” International Journal of Solids and Structures, vol. 50, no. 1, pp. 201–208, 2013.

C. Y. Dong and E. Pan, “Boundary element analysis of nanoinhomogeneities of arbitrary shapes with surface and interface effects,” Engineering Analysis with Boundary Elements, vol. 35, no. 8, pp. 996–1002, 2011.

S. S. Chen, Q. H. Li, Y. H. Liu, and Z. Q. Xue, “A meshless local natural neighbour interpolation method for analysis of two-dimensional piezoelectric structures,” Engineering Analysis with Boundary Elements, vol. 37, no. 2, pp. 273–279, 2013.

S. Chen, Y. Liu, J. Li, and Z. Cen, “Performance of the MLPG method for static shakedown analysis for bounded kinematic hardening structures,” European Journal of Mechanics, A/Solids, vol. 30, no. 2, pp. 183–194, 2011.

S. Cen, X.-R. Fu, and M.-J. Zhou, “8- and 12-node plane hybrid stress-function elements immune to severely distorted mesh containing elements with concave shapes,” Computer Methods in Applied Mechanics and Engineering, vol. 200, no. 29–32, pp. 2321–2336, 2011.

S. Cen, G.-H. Zhou, and X.-R. Fu, “A shape-free 8-node plane element unsymmetric analytical trial function method,” International Journal for Numerical Methods in Engineering, vol. 91, no. 2, pp. 158–185, 2012.

H. A. F. A. Santos, “Complementary-energy methods for geometrically non-linear structural models: an overview and recent developments in the analysis of frames,” Archives of Computational Methods in Engineering, vol. 18, no. 4, pp. 405–440, 2011.

H. A. F. A. Santos and C. I. Almeida Paulo, “On a pure complementary energy principle and a force-based finite element formulation for non-linear elastic cables,” International Journal of Non-Linear Mechanics, vol. 46, no. 2, pp. 395–406, 2011.

Y. C. Gao, “A new description of the stress state at a point with applications,” Archive of Applied Mechanics, vol. 73, no. 3-4, pp. 171–183, 2003.

Y. C. Gao, “Asymptotic analysis of the nonlinear Boussinesq problem for a kind of incompressible rubber material (compression case),” Journal of Elasticity, vol. 64, no. 2-3, pp. 111–130, 2001.

Y. C. Gao and T. J. Gao, “Large deformation contact of a rubber notch with a rigid wedge,” International Journal of Solids and Structures, vol. 37, no. 32, pp. 4319–4334, 2000.

Y. C. Gao and S. H. Chen, “Analysis of a rubber cone tensioned by a concentrated force,” Mechanics Research Communications, vol. 28, no. 1, pp. 49–54, 2001.

Y.-C. Gao, M. Jin, and G.-S. Dui, “Stresses, singularities, and a complementary energy principle for large strain elasticity,” Applied Mechanics Reviews, vol. 61, no. 3, Article ID 030801, 16 pages, 2008.