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 [14], the equilibrium models [5, 6], the mixed approach [7], the integrated force method [811], the incompatible displacement modes [12, 13], the assumed strain method [1417], 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 [2629], the new natural coordinate methods [3033], the smoothed finite element method [34], and the base force element method [3543].

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 [4851]. 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. [5456]. 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 .

For a plane rectangular coordinate system, the radius vectors and of points and can be written asin 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 aswhere 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 aswhere 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 aswhere 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 asin 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 asorwhere is the thickness of an element, is the density of material, and is acceleration of gravity.

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

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

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 asin which .

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

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.

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/s2. 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 14, 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 14, 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.

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.

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.

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/s2. 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 612, 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 612, 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.

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/m3, density of water is 1 t/m3, and acceleration of gravity  m/s2. The calculation is considered into the plane strain problem and 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 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 values of stress components and displacement components of the dam are plotted in Figures 1625, 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 1625, 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.

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.

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

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/m3, density of water is 1 t/m3, and acceleration of gravity  m/s2. 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.

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.

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.

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.