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

Authors: Yijiang Peng; Qing Guo; Zhaofeng Zhang; Yanyan Shan
Journal: Mathematical Problems in Engineering
 (2015)

---

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


---

## Body

## 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 Figure1; the compliance matrix of a base force element can be obtained as [43](1)CIJ=1+νEArIJU-ν1+νrI⊗rJ,(I,J=1,2,3,4)in which E is Young’s modulus, ν is Poisson’s ratio, A is the area of an element, U is the unit tensor, and rIJ is the dot product of radius vectors rI and rJ at points I and J.Figure 1
Four-mid-node plane element.For a plane rectangular coordinate system, the radius vectorsrI and rJ of points I and J can be written as(2)rI=xIex+yIey,rJ=xJex+yJeyin which ex,  ey are the unit vectors.Further, the compliance matrix of an element can be reduced as follows:(3)CIJ=1+νEA11+νxIxJ+yIyJex⊗ex-ν1+νxIyJex⊗ey-ν1+νyIxJey⊗ex+xIxJ+11+vyIyJey⊗ey.For a plane strain problem, it is necessary to replaceE by E/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.


### 2.2. Governing Equations

The total complementary energy of the elastic system which hasn elements can be written as(4)ΠC=∑nWCe-u-I·TI,where WCe is the complementary energy of an element and TI and u-I are the resultant force vectors and the given displacement acting on the center node I of the edge I, 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:(5)ΠCe*T,λ,λθ=ΠCeT+λ∑I=14TI+λθTI×rI,where λ=λxex+λyey and λθ are the Lagrange multipliers.For the elastic system, the modified total complementary energy function of the elastic system which containsn elements should meet the following equation by means of the modified complementary energy principle:(6)δΠC*=∑nδΠCe*(T,λ,λθ)=0.Further, (6) can be expressed as (7)∂ΠC*(T,λ,λθ)∂T=0,∂ΠC*(T,λ,λθ)∂λ=0,∂ΠC*(T,λ,λθ)∂λθ=0.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 Figure1; 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(8)σ=1ATI⊗rI,where ⊗ is the dyadic symbol, TI and rI are the resultant force vectors acting on the center node I of the edge I and the radius vector of the node I, respectively, and the summation rule is implied.For a plane rectangular coordinate system, the force vectorsTIof the node I can be written as(9)TI=TxIex+TyIey,where TxI and TyI are the components of the force vector TI along coordinates x and y, respectively.Further, the stress tensors of an element can be reduced as follows:(10)σ=1A∑I=14TxIxIex⊗ex+TxIyIex⊗ey+TyIxIey⊗ex+TyIyIey⊗ey.


### 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(11)δI=CIJ·TJ+λ+λθε·rI,in which ε is the alternating tensor, λ and λθ are the Lagrange multipliers [43], and ε and λ can be expressed as(12)ε=ex⊗ey-ey⊗ex,λ=λxex+λyey.Further, the displacement vectors of an element can be reduced as follows [43]:(13)δI=CIxJxTxJ+CIxJyTyJ+λx+λθyIex+CIyJxTxJ+CIyJyTyJ+λy-λθxIey.


## 2.1. Compliance Matrix

Consider a 4-mid-node plane element as shown in Figure1; the compliance matrix of a base force element can be obtained as [43](1)CIJ=1+νEArIJU-ν1+νrI⊗rJ,(I,J=1,2,3,4)in which E is Young’s modulus, ν is Poisson’s ratio, A is the area of an element, U is the unit tensor, and rIJ is the dot product of radius vectors rI and rJ at points I and J.Figure 1
Four-mid-node plane element.For a plane rectangular coordinate system, the radius vectorsrI and rJ of points I and J can be written as(2)rI=xIex+yIey,rJ=xJex+yJeyin which ex,  ey are the unit vectors.Further, the compliance matrix of an element can be reduced as follows:(3)CIJ=1+νEA11+νxIxJ+yIyJex⊗ex-ν1+νxIyJex⊗ey-ν1+νyIxJey⊗ex+xIxJ+11+vyIyJey⊗ey.For a plane strain problem, it is necessary to replaceE by E/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.


## 2.2. Governing Equations

The total complementary energy of the elastic system which hasn elements can be written as(4)ΠC=∑nWCe-u-I·TI,where WCe is the complementary energy of an element and TI and u-I are the resultant force vectors and the given displacement acting on the center node I of the edge I, 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:(5)ΠCe*T,λ,λθ=ΠCeT+λ∑I=14TI+λθTI×rI,where λ=λxex+λyey and λθ are the Lagrange multipliers.For the elastic system, the modified total complementary energy function of the elastic system which containsn elements should meet the following equation by means of the modified complementary energy principle:(6)δΠC*=∑nδΠCe*(T,λ,λθ)=0.Further, (6) can be expressed as (7)∂ΠC*(T,λ,λθ)∂T=0,∂ΠC*(T,λ,λθ)∂λ=0,∂ΠC*(T,λ,λθ)∂λθ=0.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 Figure1; 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(8)σ=1ATI⊗rI,where ⊗ is the dyadic symbol, TI and rI are the resultant force vectors acting on the center node I of the edge I and the radius vector of the node I, respectively, and the summation rule is implied.For a plane rectangular coordinate system, the force vectorsTIof the node I can be written as(9)TI=TxIex+TyIey,where TxI and TyI are the components of the force vector TI along coordinates x and y, respectively.Further, the stress tensors of an element can be reduced as follows:(10)σ=1A∑I=14TxIxIex⊗ex+TxIyIex⊗ey+TyIxIey⊗ex+TyIyIey⊗ey.


## 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(11)δI=CIJ·TJ+λ+λθε·rI,in which ε is the alternating tensor, λ and λθ are the Lagrange multipliers [43], and ε and λ can be expressed as(12)ε=ex⊗ey-ey⊗ex,λ=λxex+λyey.Further, the displacement vectors of an element can be reduced as follows [43]:(13)δI=CIxJxTxJ+CIxJyTyJ+λx+λθyIex+CIyJxTxJ+CIyJyTyJ+λy-λθxIey.


## 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 Figure2, and the expression can be given as(14)Qe=-ρgAt401010101Tor(15)QI=-ρgAt4ey,I=1,2,3,4,where t is the thickness of an element, ρ is the density of material, and g 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 Figure3, 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:(16)σ¯=1AT1⊗r1+T2⊗r2+T3⊗r3+T4⊗r4,where TI=TxIex+TyIey, rI=xIex+yIey,  I=1,2,3,4, xI and yI are the coordinates of node I, respectively.When the gravity of an element is taken into account and there is a free boundary, as shown in Figure4, the stress tensor of the element can be calculated as(17)σ¯=1AT1+Q1⊗r1+Q2⊗r2+T3+Q3⊗r3+T4+Q4⊗r4.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 Figure5, the stress tensor of the element can be calculated as(18)σ¯=1AT1+Q1⊗r1+F+Q2⊗r2+T3+Q3⊗r3+T4+Q4⊗r4in which F=Fxex+Fyey.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 Figure6, the stress tensor of the element can be calculated as(19)σ¯=1AT1⊗r1+F⊗r2+T3⊗r3+T4⊗r4.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.


## 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 Figure2, and the expression can be given as(14)Qe=-ρgAt401010101Tor(15)QI=-ρgAt4ey,I=1,2,3,4,where t is the thickness of an element, ρ is the density of material, and g 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 Figure3, 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:(16)σ¯=1AT1⊗r1+T2⊗r2+T3⊗r3+T4⊗r4,where TI=TxIex+TyIey, rI=xIex+yIey,  I=1,2,3,4, xI and yI are the coordinates of node I, respectively.When the gravity of an element is taken into account and there is a free boundary, as shown in Figure4, the stress tensor of the element can be calculated as(17)σ¯=1AT1+Q1⊗r1+Q2⊗r2+T3+Q3⊗r3+T4+Q4⊗r4.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 Figure5, the stress tensor of the element can be calculated as(18)σ¯=1AT1+Q1⊗r1+F+Q2⊗r2+T3+Q3⊗r3+T4+Q4⊗r4in which F=Fxex+Fyey.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 Figure6, the stress tensor of the element can be calculated as(19)σ¯=1AT1⊗r1+F⊗r2+T3⊗r3+T4⊗r4.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

WhenTn<0 and Ts<-f·Tn+c·l or 0<Tn<t0 and Ts<c·l, the contact interface is not cracked and in the elastic state. Here, Tn and Ts are the normal interface force and tangential interface force at contact interface between two elements, respectively. c and f are the cohesion and the internal friction coefficient of the contact interface, respectively. l is the length of an element. t0 is the tensile strength of an element, and it is positive in tension.


#### 4.1.2. Condition of Positive Slip at Contact Interface

WhenTn<0 and Ts≥-f·Tn+c·l, 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:(20)Ts0=f·Ts.When0<Tn<t0 and Ts≥c·l, 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

WhenTn<0 and Ts≤f·Tn-c·l, 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(21)Ts0=-f·Ts.When0<Tn<t0 and Ts≤-c·l, 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

WhenTn≥t0, 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 Figure7.Figure 7
Flow chart of the nonlinear BEFM for crack propagation problems.


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

WhenTn<0 and Ts<-f·Tn+c·l or 0<Tn<t0 and Ts<c·l, the contact interface is not cracked and in the elastic state. Here, Tn and Ts are the normal interface force and tangential interface force at contact interface between two elements, respectively. c and f are the cohesion and the internal friction coefficient of the contact interface, respectively. l is the length of an element. t0 is the tensile strength of an element, and it is positive in tension.


### 4.1.2. Condition of Positive Slip at Contact Interface

WhenTn<0 and Ts≥-f·Tn+c·l, 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:(20)Ts0=f·Ts.When0<Tn<t0 and Ts≥c·l, 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

WhenTn<0 and Ts≤f·Tn-c·l, 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(21)Ts0=-f·Ts.When0<Tn<t0 and Ts≤-c·l, 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

WhenTn≥t0, 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.1.1. Condition of Elastic State at Contact Interface

WhenTn<0 and Ts<-f·Tn+c·l or 0<Tn<t0 and Ts<c·l, the contact interface is not cracked and in the elastic state. Here, Tn and Ts are the normal interface force and tangential interface force at contact interface between two elements, respectively. c and f are the cohesion and the internal friction coefficient of the contact interface, respectively. l is the length of an element. t0 is the tensile strength of an element, and it is positive in tension.


## 4.1.2. Condition of Positive Slip at Contact Interface

WhenTn<0 and Ts≥-f·Tn+c·l, 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:(20)Ts0=f·Ts.When0<Tn<t0 and Ts≥c·l, 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

WhenTn<0 and Ts≤f·Tn-c·l, 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(21)Ts0=-f·Ts.When0<Tn<t0 and Ts≤-c·l, 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

WhenTn≥t0, 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 Figure7.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:(22)K=∑i=1n-Tnifi+cili∑i=1nTsi,where Tni and Tsi are the normal interface force and tangential interface force at contact interface between two elements along the sliding path, respectively. ci and fi are the cohesion and the internal friction coefficient of the contact interface between two elements along the sliding path, respectively. li 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 Figure8. And its width is 5 m, its height is 10 m, modulus of elasticity E=1×108Pa, Poisson ratio ν=0.3, density of rock ρ=2.45t/m3, and acceleration of gravity g=9.8 m/s2. 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 Figure9.Figure 9
Three kinds of meshes for a rock pillar.
(a)
100 elements, 220 nodes
(b)
200 elements, 430 nodes
(c)
300 elements, 640 nodesThe values of stress components and displacement components of the rock pillar are listed in Tables1–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
Displacementuy at the top of the pillar.
Meshes
BFEM model (×10−2 m)
Q4 model (×10−2 m)
Q4R model (×10−2 m)
10 × 10
−1.1939
−1.1917
−1.2048
10 × 20
−1.1940
−1.1931
−1.1964
10 × 30
−1.1940
−1.1934
−1.1940Table 2
Displacementuy at h=5 m of the pillar.
Meshes
BFEM (×10−3 m)
Q4 model (×10−3 m)
Q4R model (×10−3 m)
10 × 10
−8.7592
−8.6975
−8.8045
10 × 20
−8.7630
−8.7118
−8.6726
10 × 30
−8.7637
−8.7149
−8.7363Table 3
Stressσx 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.002Table 4
Stressσy 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 Figure10. And its modulus of elasticity E1=109, E2=108, E3=107, and E4=106, 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.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 Figure11.Figure 11
Mesh with 200 elements, 430 nodes.The values of stress components of the rock pillar are listed in Table5, 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.
σ
x
τ
x
y
σ
y
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 Figure12. And its modulus of elasticity E=108Pa, Poisson’ ratio ν=0.3, density of rock ρ1=2.4t/m3, density of water ρ2=1.0t/m3, and acceleration of gravity g=9.8 m/s2. 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 Figure13.Figure 13
Four kinds of meshes for a rock pillar.
(a)
800 elements, 1660 nodes
(b)
400 elements, 850 nodes
(c)
320 elements, 688 nodes
(d)
240 elements, 526 nodesThe values of stress components and displacement components of the rock pillar are listed in Tables6–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 σx of the element in lower right corner of rock pillar between the BFEM and the Q4R model, as shown in Table 10.Table 6
Displacementux at x=2.5 m of the pillar top.
Meshes
BFEM (×10−2 m)
Q4 model (×10−2 m)
Q4R model (×10−2 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.0863Table 7
Displacementux at h=5 m on the right of the pillar.
Meshes
BFEM (×10−2 m)
Q4 model (×10−2 m)
Q4R model (×10−2 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.0862Table 8
Displacementuy at x=2.5 m of the pillar top.
Meshes
BFEM (×10−3 m)
Q4 model (×10−3 m)
Q4R model (×10−3 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.8926Table 9
Displacementuy at h=5 m on the right of the pillar.
Meshes
BFEM model (×10−2 m)
Q4 model (×10−2 m)
Q4R model (×10−2 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.4975Table 10
Stressσx 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.941Table 11
Stressτxy 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.288Table 12
Stressσy 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 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 E1=15GPa, Poisson ratio of concrete ν1=0.2, the elastic modulus of rock E2=30GPa, Poisson ratio of rock ν2=0.3, density of concrete is 2.45 t/m3, density of water is 1 t/m3, and acceleration of gravity g=9.8 m/s2. 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 Figure15. 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 Figures16–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
Theh-σx curves at the upstream face of the dam.Figure 17
Theh-τxy curves at the upstream face of the dam.Figure 18
Theh-σy curves at the upstream face of the dam.Figure 19
Theh-ux curves at the upstream face of the dam.Figure 20
Theh-uy curves at the upstream face of the dam.Figure 21
Theh-σx curves at the half height of the dam.Figure 22
Theh-σy curves at the half height of the dam.Figure 23
Theh-τxy curves at the half height of the dam.Figure 24
Theh-ux curves at the half height of the dam.Figure 25
Theh-uy 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 Figure26. For the convenience of study, we do not consider the weight. And we use the dimensionless numerical analysis. Assuming elastic modulus E=1, Poisson ratio ν=0.3, tensile strength of a large number, and the uniform load p1=1 and p2=1. 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 Figure27.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 interfacef=0.5 and change the value of interface cohesion c. 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)
Whenc=10, there is no cracks.(2)
Whenc=2.8, there is one element interface crack.(3)
Whenc=2, there are three elements’ interfaces cracks.(4)
Whenc=1.5, there are five elements’ interfaces cracks.(5)
Whenc=1.1, there are seven elements’ interfaces cracks.(6)
Whenc=1.0, the cracks are too long, and too little structural constraints have been insufficient to solve the equations.When the interface cohesionc is 1.5, the case of crack propagation of rock block is shown in Figure 28, and the safety factor of stability is K=2 which is consistent with the result of the theoretical analysis using the rigid limit equilibrium method.Figure 28
Crack propagation path of the rock block.
(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


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

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 E=15GPa, Poisson ratio of concrete ν=0.2, density of concrete is 2.45 t/m3, density of water is 1 t/m3, and acceleration of gravity g=9.8 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 Figure29.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 Section4 and the failure criteria in Section 4.1 are as follows.(1)
Whenc=106Pa and f=0.95, there is no cracks.(2)
Whenc=0.5×106Pa and f=0.5, there are two elements’ interface cracks.(3)
Whenc=0.1×106Pa and f=0.4, there are eight elements’ interfaces cracks.The case of crack propagation of dam is shown in Figure30. When c=106Pa and f=0.5, the safety factor of stability K=4.0 which is consistent with the result of the theoretical analysis. When c=0.5×106Pa and f=0.5, the safety factor of stability K=2.55 which is consistent with the results of the theoretical analysis using the rigid limit equilibrium method.Figure 30
Crack propagation path of the gravity dam.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


#### 6.6.2. Existing an Initial Crack in Dam Foundation

There is an initial crack in the dam heel as shown in Figure31. 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)
Whenc=106Pa and f=0.95, there is no cracks.(2)
Whenc=0.5×106Pa and f=0.5, there is one element interface cracks.(3)
Whenc=0.1×106Pa and f=0.5, there are three element’ interface cracks.(4)
Whenc=0.1×106Pa and f=0.4, there are seven elements’ interfaces cracks.Figure 31
A gravity dam with a crack at the dam heel.


## 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 Figure8. And its width is 5 m, its height is 10 m, modulus of elasticity E=1×108Pa, Poisson ratio ν=0.3, density of rock ρ=2.45t/m3, and acceleration of gravity g=9.8 m/s2. 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 Figure9.Figure 9
Three kinds of meshes for a rock pillar.
(a)
100 elements, 220 nodes
(b)
200 elements, 430 nodes
(c)
300 elements, 640 nodesThe values of stress components and displacement components of the rock pillar are listed in Tables1–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
Displacementuy at the top of the pillar.
Meshes
BFEM model (×10−2 m)
Q4 model (×10−2 m)
Q4R model (×10−2 m)
10 × 10
−1.1939
−1.1917
−1.2048
10 × 20
−1.1940
−1.1931
−1.1964
10 × 30
−1.1940
−1.1934
−1.1940Table 2
Displacementuy at h=5 m of the pillar.
Meshes
BFEM (×10−3 m)
Q4 model (×10−3 m)
Q4R model (×10−3 m)
10 × 10
−8.7592
−8.6975
−8.8045
10 × 20
−8.7630
−8.7118
−8.6726
10 × 30
−8.7637
−8.7149
−8.7363Table 3
Stressσx 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.002Table 4
Stressσy 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 Figure10. And its modulus of elasticity E1=109, E2=108, E3=107, and E4=106, 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.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 Figure11.Figure 11
Mesh with 200 elements, 430 nodes.The values of stress components of the rock pillar are listed in Table5, 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.
σ
x
τ
x
y
σ
y
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 Figure12. And its modulus of elasticity E=108Pa, Poisson’ ratio ν=0.3, density of rock ρ1=2.4t/m3, density of water ρ2=1.0t/m3, and acceleration of gravity g=9.8 m/s2. 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 Figure13.Figure 13
Four kinds of meshes for a rock pillar.
(a)
800 elements, 1660 nodes
(b)
400 elements, 850 nodes
(c)
320 elements, 688 nodes
(d)
240 elements, 526 nodesThe values of stress components and displacement components of the rock pillar are listed in Tables6–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 σx of the element in lower right corner of rock pillar between the BFEM and the Q4R model, as shown in Table 10.Table 6
Displacementux at x=2.5 m of the pillar top.
Meshes
BFEM (×10−2 m)
Q4 model (×10−2 m)
Q4R model (×10−2 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.0863Table 7
Displacementux at h=5 m on the right of the pillar.
Meshes
BFEM (×10−2 m)
Q4 model (×10−2 m)
Q4R model (×10−2 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.0862Table 8
Displacementuy at x=2.5 m of the pillar top.
Meshes
BFEM (×10−3 m)
Q4 model (×10−3 m)
Q4R model (×10−3 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.8926Table 9
Displacementuy at h=5 m on the right of the pillar.
Meshes
BFEM model (×10−2 m)
Q4 model (×10−2 m)
Q4R model (×10−2 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.4975Table 10
Stressσx 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.941Table 11
Stressτxy 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.288Table 12
Stressσy 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 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 E1=15GPa, Poisson ratio of concrete ν1=0.2, the elastic modulus of rock E2=30GPa, Poisson ratio of rock ν2=0.3, density of concrete is 2.45 t/m3, density of water is 1 t/m3, and acceleration of gravity g=9.8 m/s2. 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 Figure15. 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 Figures16–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
Theh-σx curves at the upstream face of the dam.Figure 17
Theh-τxy curves at the upstream face of the dam.Figure 18
Theh-σy curves at the upstream face of the dam.Figure 19
Theh-ux curves at the upstream face of the dam.Figure 20
Theh-uy curves at the upstream face of the dam.Figure 21
Theh-σx curves at the half height of the dam.Figure 22
Theh-σy curves at the half height of the dam.Figure 23
Theh-τxy curves at the half height of the dam.Figure 24
Theh-ux curves at the half height of the dam.Figure 25
Theh-uy 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 Figure26. For the convenience of study, we do not consider the weight. And we use the dimensionless numerical analysis. Assuming elastic modulus E=1, Poisson ratio ν=0.3, tensile strength of a large number, and the uniform load p1=1 and p2=1. 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 Figure27.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 interfacef=0.5 and change the value of interface cohesion c. 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)
Whenc=10, there is no cracks.(2)
Whenc=2.8, there is one element interface crack.(3)
Whenc=2, there are three elements’ interfaces cracks.(4)
Whenc=1.5, there are five elements’ interfaces cracks.(5)
Whenc=1.1, there are seven elements’ interfaces cracks.(6)
Whenc=1.0, the cracks are too long, and too little structural constraints have been insufficient to solve the equations.When the interface cohesionc is 1.5, the case of crack propagation of rock block is shown in Figure 28, and the safety factor of stability is K=2 which is consistent with the result of the theoretical analysis using the rigid limit equilibrium method.Figure 28
Crack propagation path of the rock block.
(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


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

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 E=15GPa, Poisson ratio of concrete ν=0.2, density of concrete is 2.45 t/m3, density of water is 1 t/m3, and acceleration of gravity g=9.8 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 Figure29.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 Section4 and the failure criteria in Section 4.1 are as follows.(1)
Whenc=106Pa and f=0.95, there is no cracks.(2)
Whenc=0.5×106Pa and f=0.5, there are two elements’ interface cracks.(3)
Whenc=0.1×106Pa and f=0.4, there are eight elements’ interfaces cracks.The case of crack propagation of dam is shown in Figure30. When c=106Pa and f=0.5, the safety factor of stability K=4.0 which is consistent with the result of the theoretical analysis. When c=0.5×106Pa and f=0.5, the safety factor of stability K=2.55 which is consistent with the results of the theoretical analysis using the rigid limit equilibrium method.Figure 30
Crack propagation path of the gravity dam.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


### 6.6.2. Existing an Initial Crack in Dam Foundation

There is an initial crack in the dam heel as shown in Figure31. 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)
Whenc=106Pa and f=0.95, there is no cracks.(2)
Whenc=0.5×106Pa and f=0.5, there is one element interface cracks.(3)
Whenc=0.1×106Pa and f=0.5, there are three element’ interface cracks.(4)
Whenc=0.1×106Pa and f=0.4, there are seven elements’ interfaces cracks.Figure 31
A gravity dam with a crack at the dam heel.


## 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 Section4 and the failure criteria in Section 4.1 are as follows.(1)
Whenc=106Pa and f=0.95, there is no cracks.(2)
Whenc=0.5×106Pa and f=0.5, there are two elements’ interface cracks.(3)
Whenc=0.1×106Pa and f=0.4, there are eight elements’ interfaces cracks.The case of crack propagation of dam is shown in Figure30. When c=106Pa and f=0.5, the safety factor of stability K=4.0 which is consistent with the result of the theoretical analysis. When c=0.5×106Pa and f=0.5, the safety factor of stability K=2.55 which is consistent with the results of the theoretical analysis using the rigid limit equilibrium method.Figure 30
Crack propagation path of the gravity dam.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


## 6.6.2. Existing an Initial Crack in Dam Foundation

There is an initial crack in the dam heel as shown in Figure31. 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)
Whenc=106Pa and f=0.95, there is no cracks.(2)
Whenc=0.5×106Pa and f=0.5, there is one element interface cracks.(3)
Whenc=0.1×106Pa and f=0.5, there are three element’ interface cracks.(4)
Whenc=0.1×106Pa and f=0.4, 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.



---

*Source: 292809-2015-05-20.xml*