# Fracture Analysis of Brittle Materials Based on Nonlinear FEM and Application in Arch Dam with Fractures

Authors: Yuanwei Pan; Yaoru Liu; Zhixiong Cui; Xin Chen; Qiang Yang
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Current fracture analysis models based on fracture mechanics or continuum damage mechanics are still limited in the application to three-dimensional structure. Based on deformation reinforcement theory coming from elastoperfect plastic theory, unbalanced force is proposed to predict initiation and propagation of cracks. Unbalanced force is the driving force of time-dependent deformation according to Perzyna’s viscoplasticity theory. It is also related to the damage driving force in viscoplastic damage model. The distribution of unbalanced force indicates cracks initiation area, while its direction predicts possible cracks propagation path. Uniaxial compression test of precrack specimen is performed as verification to this method. The trend and distribution of cracks are in good agreement with numerical results, proving that unbalanced force is feasible and effective for fracture analysis. The method is applied in fracture analysis of Xiaowan high arch dam, which is subjected to some cracks in dam due to the temperature control program. The results show that the deformation and stress of cracks and the stress characteristics of dam are insensitive to grouting of cracks. The existing cracks are stable and dam heel is still the most possible cracking position.


---

## Body

## 1. Introduction

Rock and concrete are heterogeneous anisotropic materials, containing numerous microcosmic voids and flaws. Fracture is a common and significant failure mode of geotechnical structure. Fracture evaluation for structure under certain loads, including crack initiation, propagation, and penetration, is still an unsolved problem. Thus fracture analysis method for rock and concrete structure is of significant importance in the sense of cracking prevention and global stability evaluation.There is a key problem remaining in fracture analysis for brittle materials and structures, that is, a feasible fracture criterion. Common fracture criterions include stress criterions and energetic criterions [1–8]. The former state that failure occurs when the maximum principal stress in some local point exceeds the tensile strength. The latter are provided by linear elastic fracture mechanics, which covers some precise measurement, such as the stress intensity factor (SIF), energy release rate, and strain energy density. The stress criterions are acceptable for bodies without cracks, while the energetic criterions only work when a certain large crack already exists.Current numerical methods on fracture analysis include fracture mechanics [9] and continuum damage mechanics [10]. Some researchers combine both methods for numerical failure analysis [11, 12]. Crack propagation simulation could be concluded as smeared fracture model [1] and discrete fracture model [13]. Besides, numerical tools have developed rapidly, including Finite Element Method (FEM) [14], Extended FEM [15], Element Free Method [16], Boundary Element Method [17], Discrete Element Method [18], Numerical Manifold Method [19], Discontinuous Deformation Analysis [20], and Fast Lagrangian Method [21].The theories and methods mentioned above work well in planar analysis. There is still severe limitation on the applicability when extended to three-dimensional and complex structure. Besides, behavior of rock and concrete involves complex nonlinear overall deformation, which is beyond the capacity of common numerical methods.This paper presents a new approach based on deformation reinforcement theory [22] to evaluate cracks growth in three-dimensional structure. Unbalanced force is proposed as a set of equivalent nodal forces of overstress beyond the yield surface. Unbalanced force drives time-dependent deformation, as well as damage evolution. The distribution of unbalanced force indicates crack initiation area, while its direction predicts potential crack propagation path. Uniaxial compression test of precrack specimen is performed, which is in a good agreement with the numerical results. The method is also applied in fracture analysis of Xiaowan high arch dam.


## 2. Fracture Analysis Method Based on Unbalanced Force

### 2.1. Deformation Reinforcement Theory

Most geotechnical structures are under complicated configurations and working conditions. Stability analysis of these engineering structure could be summarized as a boundary value problem with some basic equations, including kinematic admissibility, equilibrium condition, and constitutive equations. The classical elastoplastic theory aims at solving the displacement and stress fields that simultaneously satisfy all the aforementioned equations. However, the existence of such solution requires that the structure is stable. Structural instability occurs when action is greater than resistance. The difference between action and resistance is overstress, which is the key concept of the Deformation Reinforcement Theory (DRT) and defined as the unbalanced force.For perfect elasto-plastic materials with associative flow rule, the constitutive equations can be stated as(1)ε˙=ε˙e+ε˙p,ε˙e=C:σ˙;Dp(σ,ε˙p)≥Dp(σyc,ε˙p),∀f(σyc)≤0,
where Dp is defined as the plastic dissipation, Dp(σyc,ε˙p)∶=σyc:ε˙p. ε˙e and ε˙p are elastic and plastic strain rates, respectively. C is the flexibility tensor. f is the yield function and σyc is an arbitrary stress state satisfying the yield criterion. Inequality (1) is known as principle of maximum plastic dissipation [23], which covers associative flow rule and the Kuhn-Tucker condition.Assuming thatt is a pseudotime, the stress and plastic strain of material are σ0 and ε0p at time t0. With a strain increment in time interval Δt, a new state σ and εp is achieved. Thus the linearized plastic strain rate is
(2)ε˙p=ΔεpΔt=Δε-C:(σ-σ0)Δt=C:(σeq-σ)Δt,
where σeq is termed the trial elastic stress while ε˙p=0. σeq is also an equilibrium stress field under given loads;
(3)σeq=σ0+C-1:Δε.Applying (2) into Inequality (1), the process of solving real stress field is turned into the following minimization problem:
(4)minE(σyc),∀f(σyc)≤0,E(σyc)=12(σeq-σyc):C:(σeq-σyc).Equation (4) is known as the closest-point projection method (CPPM) [24], as shown in Figure 1. σyc is an arbitrary stress field on the yielding surface, which represents the material resistance after previous minimization process. σeq is an equilibrium stress field, which could be regarded as certain stress under external actions.Figure 1
The closest-point projection method.The minimization objectiveE represents the difference between the plastic dissipations of the external action and the material resistance. It is defined as the volume density of the plastic complementary energy (PCE);
(5)E=Dp(σeq,Δεp)-Dp(σyc,Δεp).Thus, instability of a material point is equal to the statement that the external action is greater than the material resistance:E>0, for all f(σyc)≤0. As indicated by the minimization problem, if σeq>σyc, the real stress state minimizes the resistance deficiency in the sense of PCE. In other words, material resistance capacity is fully developed.Drucker-Prager yield criterion is adopted in the following nonlinear calculation:(6)f(σ)=αI1+J2-H,
where I1=trσ, J2=σ:σ/2+I1/6, and  α and H are material parameters. With the real stress field σ solved by (4), the following characteristic can be proved [22]:
(7)∂f∂σ|σ=σ=∂f∂σ|σ=σeq.Eventually the following analytical solution ofσ can be derived:
(8)σ=(1-n)σeq+pI,n=wGJ2,p=-mw+3nI1,m=α(3λ+2G),w=f(3αm+G),
where I1 and J2 are invariants of σeq and λ and G are the Lamé constants.


### 2.2. Demonstration in FEM

PCE and unbalanced force are implemented in FEM analysis. Since the process is restricted to displacement method, the kinematic admissibility is naturally satisfied.In order to evaluate the global stability of structure under certain actions, two stress field sets,Seq and Syc, are defined for elastoperfect plastic problem. The two stress sets, respectively, satisfy the equilibrium condition and the yield criterion as
(9)Seq={σeq∣F=∫VBTσeqdV},Syc={σyc∣f(σyc)≤0},
where B is the displacement gradient matrix, F is equivalent nodal force vector of external loads, and V is structure volume. Seq and Syc, respectively, represent external actions and structural resistance. The real stress field must satisfy the yield criterion; σ=σyc. Equivalent nodal forces of the difference between the two stress fields, which is the driving force of structure deformation, can be defined as unbalanced force as
(10)ΔU=∫VBT(σeq-σyc)dV=F-∫VBTσdV.For a given load step, the unbalanced force is the driving force of deformation in the FE iterations. The principle of minimum PCE implies that elasto-plastic structures deform tending to the state where PCE is minimized under certain actions as(11)minE(σeq,σyc),∀σeq∈Seq,σyc∈Syc,E(σeq,σyc)=∫V12(σeq-σyc):C:(σeq-σyc)dV.Unbalanced force is the measurement between external actions and structural resistance. The plastic complementary energyE is a norm of unbalanced force. Once E reaches its minimum value Emin, the unbalanced force stays constant and the structure suffers steady plastic flow until failure occurs. Only if Emin=0, the structure is stable.In elastoplastic FEM iteration, PCE and unbalanced force can be calculated by the following steps. First, the nearest stable stress fieldσ1yc can be achieved from an arbitrary equilibrium stress field σ1eq. Then, the nearest equilibrium stress field σ2eq is also obtained from σ1yc. Finally, structure stress state tends to the two closest stress fields σ and σeq, as shown in Figure 2. The iteration converges if the plastic complementary energy is reaching a steady value ΔE0. For ideal model, ΔE0 is equal to its minimum value ΔEmin. If ΔEmin=0, the structure remains stable and the stress field σ is the real stress response. Otherwise, the structure is unstable and PCE is a magnitude estimation of its global instability. Unbalanced force reflects the structural failure behavior, including failure position and pattern.Figure 2
Demonstration of the elastoplastic iteration in FEM analysis.


### 2.3. Explanation of Unbalanced Force in Viscoplastic Damage Model

Stress state beyond the yielding surface is unacceptable for perfect elastoplasticity. It exists only in the iteration process of FEM calculation. However, in viscoplasticity models, it is of explicit physical significance as the driving force of visco-plasticity deformation. The Perzyna associative visco-plasticity strain rate could be stated as [25, 26]
(12)ε˙ijvp={0,f≤0,ΓvpΦ(f)∂Q∂σ,f>0,
where ε˙ijvp is the visco-plasticity strain rate. Γvp is the viscosity parameter. f=f(σ,κ)=0 is the yield function. Q is plastic potential function, and for associative flow rule, Q=f(σ). Φ(f) is overstress function characterized by the yield function, commonly expressed as Φ(f)=(f/f0)n. f0 is a reference constant of the same dimension with yield function f.According to Perzyna’s visco-plasticity theory, visco-plasticity strain rate is in proportion to overstress function.∂f/∂σ represents the direction of ε˙ijvp, while Φ(f) reflects the magnitude of ε˙ijvp. If the overstress function is restricted to yield function, plastic flow occurs only on the condition that the yield function value is positive. As the DRT  states, unbalanced force is a set of the equivalent nodal forces of stress exceeding the yield function, which can be termed as the driving force of time-dependent deformation.Unbalanced force is also involved in structural crack growth analysis based on damage mechanics methods. In the concept of continuum damage mechanics, the damage variable is usually related to equivalent plastic strain, that is, the accumulation of time-dependent deformation. A classical definition of damage variable isϕ=1-A-/A, where A is total area and A- is intact area [27, 28]. Effective stress tensor in the undamaged configuration is determined by the damage variable and the nominal Cauchy stress tensor [29, 30] as
(13)σ-ij=σij(1-ϕ)2,
where σ-ij is the effective stress tensor and σij is the nominal Cauchy stress tensor.Darabi et al. [31] presented a thermoviscodamage model to describe the damage evolution law, which could be stated as
(14)ϕ˙=Γ0φ(Y-(1-ϕ)2Y0)qexp(kε-effTot)G(T),
where ε-effTot=ε-ijε-ij is the effective total strain in the effective configuration. ε-ij includes both viscoelastic and viscoplastic components. q is the stress dependency parameter and k is a material parameter. Γ0φ and Y0 are the reference damage viscosity parameter and the reference damage force, respectively. G(T) is a damage temperature function. Y- is the damage driving force in the effective configuration, which can be assumed to have a modified Drucker-Prager-type form as
(15)Y-=αI-1+τ-,τ-=J-22[1+1d+(1-1d)J-3J-23],
where I-1 and J-2 are the invariants of the effective stress tensor. d gives the distinction of material behavior in compression and extension loading conditions, where d=1 implies that τ-=J-2. Thus the damage driving force is expressed as
(16)Y-=αI-1+J-2.Damage evolution requires thatY->0; that is, the equilibrium stress field exceeds the yield function. As stated above, the equivalent nodal force of stress exceeding the yield function is unbalanced force, which is the driving force of damage evolution. Unbalanced force is more explicit than equivalent plastic strain in the sense of physical significance.


### 2.4. Fracture Analysis Method

Unbalanced force is a set of equivalent nodal forces of plastic stress in the elements around a node, which reflects the difference between external actions and structural resistance. The process of FEM iteration is to find a set of additional forces to prevent failure from occurring, which is minimized in the sense of PCE. If the additional force, that is, unbalanced force, tends to be zero in the iteration, the cracks will not initiate or stay in a limit state. Otherwise, cracks will initiate and propagate, and the direction and distribution of unbalanced force indicate the potential failure area.Thus, occurrence of unbalanced force can be the identification of local cracks initiation. Meanwhile, as unbalanced force is related to the damage driving force in the thermo-viscodamage model, the distribution of unbalanced force will expand gradually as material damage accumulates. This process corresponds to the cracks growth and propagation. In general, its direction predicts the possible path of crack’s propagation.Current fracture analysis methods are based on planar problems, which provide good results dealing with a single crack or two cracks. When extended to 3D structures, the fracture criterion and nonlinear calculation efficiency remain to be solved. Unbalanced force is presented in this paper, which is of clear physical concept and could be effectively applied to elato-plastic FEM calculation.


## 2.1. Deformation Reinforcement Theory

Most geotechnical structures are under complicated configurations and working conditions. Stability analysis of these engineering structure could be summarized as a boundary value problem with some basic equations, including kinematic admissibility, equilibrium condition, and constitutive equations. The classical elastoplastic theory aims at solving the displacement and stress fields that simultaneously satisfy all the aforementioned equations. However, the existence of such solution requires that the structure is stable. Structural instability occurs when action is greater than resistance. The difference between action and resistance is overstress, which is the key concept of the Deformation Reinforcement Theory (DRT) and defined as the unbalanced force.For perfect elasto-plastic materials with associative flow rule, the constitutive equations can be stated as(1)ε˙=ε˙e+ε˙p,ε˙e=C:σ˙;Dp(σ,ε˙p)≥Dp(σyc,ε˙p),∀f(σyc)≤0,
where Dp is defined as the plastic dissipation, Dp(σyc,ε˙p)∶=σyc:ε˙p. ε˙e and ε˙p are elastic and plastic strain rates, respectively. C is the flexibility tensor. f is the yield function and σyc is an arbitrary stress state satisfying the yield criterion. Inequality (1) is known as principle of maximum plastic dissipation [23], which covers associative flow rule and the Kuhn-Tucker condition.Assuming thatt is a pseudotime, the stress and plastic strain of material are σ0 and ε0p at time t0. With a strain increment in time interval Δt, a new state σ and εp is achieved. Thus the linearized plastic strain rate is
(2)ε˙p=ΔεpΔt=Δε-C:(σ-σ0)Δt=C:(σeq-σ)Δt,
where σeq is termed the trial elastic stress while ε˙p=0. σeq is also an equilibrium stress field under given loads;
(3)σeq=σ0+C-1:Δε.Applying (2) into Inequality (1), the process of solving real stress field is turned into the following minimization problem:
(4)minE(σyc),∀f(σyc)≤0,E(σyc)=12(σeq-σyc):C:(σeq-σyc).Equation (4) is known as the closest-point projection method (CPPM) [24], as shown in Figure 1. σyc is an arbitrary stress field on the yielding surface, which represents the material resistance after previous minimization process. σeq is an equilibrium stress field, which could be regarded as certain stress under external actions.Figure 1
The closest-point projection method.The minimization objectiveE represents the difference between the plastic dissipations of the external action and the material resistance. It is defined as the volume density of the plastic complementary energy (PCE);
(5)E=Dp(σeq,Δεp)-Dp(σyc,Δεp).Thus, instability of a material point is equal to the statement that the external action is greater than the material resistance:E>0, for all f(σyc)≤0. As indicated by the minimization problem, if σeq>σyc, the real stress state minimizes the resistance deficiency in the sense of PCE. In other words, material resistance capacity is fully developed.Drucker-Prager yield criterion is adopted in the following nonlinear calculation:(6)f(σ)=αI1+J2-H,
where I1=trσ, J2=σ:σ/2+I1/6, and  α and H are material parameters. With the real stress field σ solved by (4), the following characteristic can be proved [22]:
(7)∂f∂σ|σ=σ=∂f∂σ|σ=σeq.Eventually the following analytical solution ofσ can be derived:
(8)σ=(1-n)σeq+pI,n=wGJ2,p=-mw+3nI1,m=α(3λ+2G),w=f(3αm+G),
where I1 and J2 are invariants of σeq and λ and G are the Lamé constants.


## 2.2. Demonstration in FEM

PCE and unbalanced force are implemented in FEM analysis. Since the process is restricted to displacement method, the kinematic admissibility is naturally satisfied.In order to evaluate the global stability of structure under certain actions, two stress field sets,Seq and Syc, are defined for elastoperfect plastic problem. The two stress sets, respectively, satisfy the equilibrium condition and the yield criterion as
(9)Seq={σeq∣F=∫VBTσeqdV},Syc={σyc∣f(σyc)≤0},
where B is the displacement gradient matrix, F is equivalent nodal force vector of external loads, and V is structure volume. Seq and Syc, respectively, represent external actions and structural resistance. The real stress field must satisfy the yield criterion; σ=σyc. Equivalent nodal forces of the difference between the two stress fields, which is the driving force of structure deformation, can be defined as unbalanced force as
(10)ΔU=∫VBT(σeq-σyc)dV=F-∫VBTσdV.For a given load step, the unbalanced force is the driving force of deformation in the FE iterations. The principle of minimum PCE implies that elasto-plastic structures deform tending to the state where PCE is minimized under certain actions as(11)minE(σeq,σyc),∀σeq∈Seq,σyc∈Syc,E(σeq,σyc)=∫V12(σeq-σyc):C:(σeq-σyc)dV.Unbalanced force is the measurement between external actions and structural resistance. The plastic complementary energyE is a norm of unbalanced force. Once E reaches its minimum value Emin, the unbalanced force stays constant and the structure suffers steady plastic flow until failure occurs. Only if Emin=0, the structure is stable.In elastoplastic FEM iteration, PCE and unbalanced force can be calculated by the following steps. First, the nearest stable stress fieldσ1yc can be achieved from an arbitrary equilibrium stress field σ1eq. Then, the nearest equilibrium stress field σ2eq is also obtained from σ1yc. Finally, structure stress state tends to the two closest stress fields σ and σeq, as shown in Figure 2. The iteration converges if the plastic complementary energy is reaching a steady value ΔE0. For ideal model, ΔE0 is equal to its minimum value ΔEmin. If ΔEmin=0, the structure remains stable and the stress field σ is the real stress response. Otherwise, the structure is unstable and PCE is a magnitude estimation of its global instability. Unbalanced force reflects the structural failure behavior, including failure position and pattern.Figure 2
Demonstration of the elastoplastic iteration in FEM analysis.


## 2.3. Explanation of Unbalanced Force in Viscoplastic Damage Model

Stress state beyond the yielding surface is unacceptable for perfect elastoplasticity. It exists only in the iteration process of FEM calculation. However, in viscoplasticity models, it is of explicit physical significance as the driving force of visco-plasticity deformation. The Perzyna associative visco-plasticity strain rate could be stated as [25, 26]
(12)ε˙ijvp={0,f≤0,ΓvpΦ(f)∂Q∂σ,f>0,
where ε˙ijvp is the visco-plasticity strain rate. Γvp is the viscosity parameter. f=f(σ,κ)=0 is the yield function. Q is plastic potential function, and for associative flow rule, Q=f(σ). Φ(f) is overstress function characterized by the yield function, commonly expressed as Φ(f)=(f/f0)n. f0 is a reference constant of the same dimension with yield function f.According to Perzyna’s visco-plasticity theory, visco-plasticity strain rate is in proportion to overstress function.∂f/∂σ represents the direction of ε˙ijvp, while Φ(f) reflects the magnitude of ε˙ijvp. If the overstress function is restricted to yield function, plastic flow occurs only on the condition that the yield function value is positive. As the DRT  states, unbalanced force is a set of the equivalent nodal forces of stress exceeding the yield function, which can be termed as the driving force of time-dependent deformation.Unbalanced force is also involved in structural crack growth analysis based on damage mechanics methods. In the concept of continuum damage mechanics, the damage variable is usually related to equivalent plastic strain, that is, the accumulation of time-dependent deformation. A classical definition of damage variable isϕ=1-A-/A, where A is total area and A- is intact area [27, 28]. Effective stress tensor in the undamaged configuration is determined by the damage variable and the nominal Cauchy stress tensor [29, 30] as
(13)σ-ij=σij(1-ϕ)2,
where σ-ij is the effective stress tensor and σij is the nominal Cauchy stress tensor.Darabi et al. [31] presented a thermoviscodamage model to describe the damage evolution law, which could be stated as
(14)ϕ˙=Γ0φ(Y-(1-ϕ)2Y0)qexp(kε-effTot)G(T),
where ε-effTot=ε-ijε-ij is the effective total strain in the effective configuration. ε-ij includes both viscoelastic and viscoplastic components. q is the stress dependency parameter and k is a material parameter. Γ0φ and Y0 are the reference damage viscosity parameter and the reference damage force, respectively. G(T) is a damage temperature function. Y- is the damage driving force in the effective configuration, which can be assumed to have a modified Drucker-Prager-type form as
(15)Y-=αI-1+τ-,τ-=J-22[1+1d+(1-1d)J-3J-23],
where I-1 and J-2 are the invariants of the effective stress tensor. d gives the distinction of material behavior in compression and extension loading conditions, where d=1 implies that τ-=J-2. Thus the damage driving force is expressed as
(16)Y-=αI-1+J-2.Damage evolution requires thatY->0; that is, the equilibrium stress field exceeds the yield function. As stated above, the equivalent nodal force of stress exceeding the yield function is unbalanced force, which is the driving force of damage evolution. Unbalanced force is more explicit than equivalent plastic strain in the sense of physical significance.


## 2.4. Fracture Analysis Method

Unbalanced force is a set of equivalent nodal forces of plastic stress in the elements around a node, which reflects the difference between external actions and structural resistance. The process of FEM iteration is to find a set of additional forces to prevent failure from occurring, which is minimized in the sense of PCE. If the additional force, that is, unbalanced force, tends to be zero in the iteration, the cracks will not initiate or stay in a limit state. Otherwise, cracks will initiate and propagate, and the direction and distribution of unbalanced force indicate the potential failure area.Thus, occurrence of unbalanced force can be the identification of local cracks initiation. Meanwhile, as unbalanced force is related to the damage driving force in the thermo-viscodamage model, the distribution of unbalanced force will expand gradually as material damage accumulates. This process corresponds to the cracks growth and propagation. In general, its direction predicts the possible path of crack’s propagation.Current fracture analysis methods are based on planar problems, which provide good results dealing with a single crack or two cracks. When extended to 3D structures, the fracture criterion and nonlinear calculation efficiency remain to be solved. Unbalanced force is presented in this paper, which is of clear physical concept and could be effectively applied to elato-plastic FEM calculation.


## 3. Fracture Test and Numerical Analysis of Precrack Specimen

The strong correlation is found between the distribution of unbalanced force and the initial point of structural cracking, which is verified by the following model test and numerical results.


### 3.1. Specimen and Numerical Model

The precrack cubic specimen is made with gypsum, which is 15 cm wide, 15 cm high, and 5 cm thick, as shown in Figure3. Steel slices are inserted during the specimen modeling process and pulled out once the initial set begins. By this means the precracks are manufactured. The angle between precrack and horizontal direction is 30°, and precrack length is 12 mm.Figure 3
Sketch map of the precrack specimen.Numerical model is exactly established according to the physical test. Precracks are simulated where the corresponding elements are set null. The precrack specimen and numerical model are illustrated in Figure4.The precrack specimen and numerical model.
(a)
(b)


### 3.2. Test and Calculation Results

Uniaxial compression test of the precrack specimen is performed with its bottom restrained. The pressure on the top is gradually increased from 0 MPa to 2.0 MPa by the increment of 0.1 MPa.The precracks in the left part of specimen initiate first during the uniaxial compression test, when the pressure achieves 1.6 MPa. The cracks initiation direction is almost perpendicular to the precracks trend. With the pressure being gradually increased, cracks propagation path tends to be vertical, and precracks at the top right corner begin to initiate as well. Cracks propagate and eventually penetrate both sides of the specimen, where failure occurs with a compressive strength 2.0 MPa. The cracks growth process and final failure mode are shown in Figure5.Cracks growth and failure mode of specimen test.
(a)
(b)
(c)
(d)Nonlinear numerical calculation is applied with the same constraint conditions and loading procedure. As clearly shown in (14), the damage evolution law is influenced by various factors, including stress, strain, strain rate, temperature, and damage history. The loading process duration of this test is relatively short, so the temperature effect is ignored. Since gypsum is quasi-brittle material, the time-dependent deformation of damage evolution process is simplified in the numerical calculation. For the region where the damage driving force Y- is greater than the reference damage force Y0, gypsum is regarded as damaged and elastic modulus E,  and shear strengths f and c are reduced by 50%. In this case, Y0=10MPa. Table 1 shows the materials parameters.Table 1
Materials parameters.
E/GPa
μ
r/ kN·m−3
f
c/MPa
Gypsum
5
0.35
2
0.1
0.2
Damaged gypsum
2.5
0.35
2
0.05
0.1Unbalanced force distribution and evolution process at crack tips are shown in Figure6. When the specimen stays in a stress state of elasticity, there is no unbalanced force. With the pressure being increased, unbalanced force occurs at the precrack tips around the middle left part of the specimen. As the distribution area expands, the direction of unbalanced force growth tends to be perpendicular to the precracks. The unbalanced force eventually goes through the nearest crack tips, which indicates the most possible propagation path of specimen fracture.Unbalanced force evolution process at crack tips.
(a)
(b)
(c)
(d)Damaged area expanding process is shown in Figure7. Damaged area appears around the tip of crack and gradually extends to neighbor elements in the perpendicular direction to the precracks. At last, the damaged area unites between the tips of the two parallel precracks.Damaged area expansion process at crack tips.
(a)
(b)
(c)
(d)The existence of unbalanced force indicates that the model is unable to balance the loads and fracture occurs, which agrees with the test results. As cracks propagate, a new structure is achieved and the stress field is redistributed. Since the crack propagation is a local and quasi-static failure process, the new structure retains its bearing capacity. Thus the model is able to sustain the loads and reaches a new equilibrium state. This process continues until the load reaches the compressive strength and structural failure occurs. The distribution of unbalanced force indicates cracks initiation area, while its direction predicts the possible cracks propagation path.


## 3.1. Specimen and Numerical Model

The precrack cubic specimen is made with gypsum, which is 15 cm wide, 15 cm high, and 5 cm thick, as shown in Figure3. Steel slices are inserted during the specimen modeling process and pulled out once the initial set begins. By this means the precracks are manufactured. The angle between precrack and horizontal direction is 30°, and precrack length is 12 mm.Figure 3
Sketch map of the precrack specimen.Numerical model is exactly established according to the physical test. Precracks are simulated where the corresponding elements are set null. The precrack specimen and numerical model are illustrated in Figure4.The precrack specimen and numerical model.
(a)
(b)


## 3.2. Test and Calculation Results

Uniaxial compression test of the precrack specimen is performed with its bottom restrained. The pressure on the top is gradually increased from 0 MPa to 2.0 MPa by the increment of 0.1 MPa.The precracks in the left part of specimen initiate first during the uniaxial compression test, when the pressure achieves 1.6 MPa. The cracks initiation direction is almost perpendicular to the precracks trend. With the pressure being gradually increased, cracks propagation path tends to be vertical, and precracks at the top right corner begin to initiate as well. Cracks propagate and eventually penetrate both sides of the specimen, where failure occurs with a compressive strength 2.0 MPa. The cracks growth process and final failure mode are shown in Figure5.Cracks growth and failure mode of specimen test.
(a)
(b)
(c)
(d)Nonlinear numerical calculation is applied with the same constraint conditions and loading procedure. As clearly shown in (14), the damage evolution law is influenced by various factors, including stress, strain, strain rate, temperature, and damage history. The loading process duration of this test is relatively short, so the temperature effect is ignored. Since gypsum is quasi-brittle material, the time-dependent deformation of damage evolution process is simplified in the numerical calculation. For the region where the damage driving force Y- is greater than the reference damage force Y0, gypsum is regarded as damaged and elastic modulus E,  and shear strengths f and c are reduced by 50%. In this case, Y0=10MPa. Table 1 shows the materials parameters.Table 1
Materials parameters.
E/GPa
μ
r/ kN·m−3
f
c/MPa
Gypsum
5
0.35
2
0.1
0.2
Damaged gypsum
2.5
0.35
2
0.05
0.1Unbalanced force distribution and evolution process at crack tips are shown in Figure6. When the specimen stays in a stress state of elasticity, there is no unbalanced force. With the pressure being increased, unbalanced force occurs at the precrack tips around the middle left part of the specimen. As the distribution area expands, the direction of unbalanced force growth tends to be perpendicular to the precracks. The unbalanced force eventually goes through the nearest crack tips, which indicates the most possible propagation path of specimen fracture.Unbalanced force evolution process at crack tips.
(a)
(b)
(c)
(d)Damaged area expanding process is shown in Figure7. Damaged area appears around the tip of crack and gradually extends to neighbor elements in the perpendicular direction to the precracks. At last, the damaged area unites between the tips of the two parallel precracks.Damaged area expansion process at crack tips.
(a)
(b)
(c)
(d)The existence of unbalanced force indicates that the model is unable to balance the loads and fracture occurs, which agrees with the test results. As cracks propagate, a new structure is achieved and the stress field is redistributed. Since the crack propagation is a local and quasi-static failure process, the new structure retains its bearing capacity. Thus the model is able to sustain the loads and reaches a new equilibrium state. This process continues until the load reaches the compressive strength and structural failure occurs. The distribution of unbalanced force indicates cracks initiation area, while its direction predicts the possible cracks propagation path.


## 4. Fracture Analysis of Arch Dam with Cracks

### 4.1. Numerical Model

Xiaowan arch dam, with a height of 292 m, is subjected to some cracks in dam due to the temperature control program. Most of the cracks distribute on the middle and bottom elevations of the dam and tend tangentially to the arch. The average width of these cracks is 1 mm. Distribution of cracks on typical elevation and on typical dam section is shown in Figures8 and 9, respectively.Figure 8
Distribution of cracks on typical elevation.Figure 9
Distribution of cracks in typical dam section.The size of numerical model is as follows: upstream 300 m, downstream 900 m, and left and right banks 700 m each. Elevations of dam crest and model bottom are∇1245 m and ∇595 m, and thus the model height is 650 m. The model includes 58989 nodes and 53850 elements, as shown in Figure 10. There are 11 major cracks that are numbered and simulated, as shown in Figure 11 and Table 2.Table 2
Numbers of cracks in the dam.
Number
Crack
1
13-1
2
13-2
3
20-1
4
20-2
5
22-1
6
22-3
7
22-4
8
25-1
9
28-1
10
28-2
11
30-1
Crack 13-1 in the table means the 1st crack in dam section number 13, and so on.Figure 10
Numerical model of Xiaowan arch dam.Figure 11
Distribution of cracks in the dam.The dam cracks are induced by the improper temperature control program during construction. Since concrete grouting in the cracks is operated, the dam cracks are simulated with thin layer elements. The material parameters are shown in Table3.Table 3
Material parameters.
E/GPa
μ
r/kN·m−3
f
c/MPa
Foundation
18
0.28
28
1.52
1.33
Dam
23.1
0.2
24
1.4
1.6
Crack
5
0.3
27.5
1.12
0.9


### 4.2. Fracture Analysis

Structural fracture and parameters sensitivity analysis are presented. Three schemes are performed, including normal parameters (scheme 1), crack’s shear strengthsfand c reduced by 50% (scheme 2), and crack’s shear strengths fand c reduced by 75% (scheme 3). Both the hydrostatic pressure P0 and overload conditions are included in each scheme.


#### 4.2.1. Dam Heel Cracking Analysis

The distribution of unbalanced force at dam heel under different loads in scheme 1 is shown in Figure12. The other two schemes are of similar distribution. Under the hydrostatic pressure P0, unbalanced force rarely occurs and dam heel stays as a global steady state. As the hydrostatic pressure is overloaded, unbalanced force area expands and its magnitude increases rapidly, which is concentrated in the dam heel.Distribution of dam heel unbalanced force of scheme 1.
(a)
(b)
(c)Unbalanced forces of dam heel in three schemes are summarized in Table4. In the overload conditions, unbalanced force increases rapidly in dam heel, which is the most possible location of cracks initiation. As the variation of unbalanced forces between the three schemes is relatively small, the dam heel cracking behavior is insensitive to the change of crack parameters.Table 4
Unbalanced forces of dam heel (104 N).
Scheme
P
0
1.5P0
2.0P0
1
40249.56
257289.9
415601.5
2
40216.39
257269.2
415560.7
3
40126.59
256869.2
414973.0


#### 4.2.2. Dam Body Cracking Analysis

The distribution of unbalanced force in dam under different loads in scheme 1 is shown in Figure13. Under the hydrostatic pressure P0, there is a little unbalanced force existing in the right abutment of the arch dam. As the pressure is overloaded, unbalanced force occurs in both the dam abutments. Distribution of unbalanced force in dam is close to the foundation, where cracks initiate and failure may occur.Distribution of unbalanced force in dam of scheme 1.
(a)
(b)
(c)Distribution of dam body unbalanced force under overload2.0P0 in schemes 2 and  3 is shown in Figures 14 and 15, respectively. Dam unbalanced force is insensitive to the change of crack parameters. Crack number 10, that is, the second crack in dam section number 28, suffers from some unbalanced force in scheme 3, which indicates the most sensitive and dangerous crack in the dam. Meanwhile, the unbalanced force of cracks is much less than that of dam body, and the latter is still the emphasis of dam cracking prevention.Figure 14
Distribution of dam body unbalanced force under overload2.0P0 in scheme 2.Figure 15
Distribution of dam body unbalanced force under overload2.0P0 in scheme 3.


#### 4.2.3. Analysis for Existing Cracks

The length of unbalanced force vectors in the cracks is enlarged 20 times since it is relatively less than that of dam body. Distribution of cracks unbalanced force under overload2.0P0 in the three schemes is shown in Figures 16, 17, and 18, respectively. Cracks initiation and propagation are very sensitive to the crack parameters. In other words, the stability of the existing cracks relies on the quality of concrete grouting in the cracks.Figure 16
Distribution of cracks unbalanced force under overload2.0P0 in scheme 1.Figure 17
Distribution of cracks unbalanced force under overload2.0P0 in scheme 2.Figure 18
Distribution of cracks unbalanced force under overload2.0P0 in scheme 3.Unbalanced force of cracks under overload2.0P0 in all three schemes is summarized in Table 5. Unbalanced force in dam heel is also listed in Table 5 as comparison. Unbalanced forces in dam heel increase earlier than the cracks in the dam. Dam heel contributes most of the unbalanced forces in all conditions. Namely, dam heel cracking occurs before any existing crack propagates. Among all existing cracks, 20-2 and 28-2 are the dominating cracks in the process of fracture propagation, while 30-1 and 13-2 are also possible initiation area, as shown in Figure 19. The cracks near the crown cantilever are stable.Table 5
Unbalanced force of cracks under overload2P0 (104 N).
Number
Crack
Scheme   1
Scheme   2
Scheme   3
1
13-1
0.00
0.00
12.54
2
13-2
0.02
0.87
221.69
3
20-1
0.00
0.00
2.83
4
20-2
0.03
0.70
7614.90
*
5
22-1
0.00
0.00
0.00
6
22-3
0.00
0.00
0.00
7
22-4
0.00
0.00
0.00
8
25-1
0.00
0.00
8.00
9
28-1
0.00
0.94
65.95
10
28-2
1.77
345.54
15810.91
*
11
30-1
0.02
6.04
687.85
Dam heel
415601.5
415560.7
414973
*Dominating cracks in the propagation process.Figure 19
Local distribution of cracks unbalanced force (scheme 1, under overload2.0P0).


### 4.3. Comparison with Geomechanical Model Test

Results and conclusions of geomechanical model test are presented as comparison with numerical method. The final distribution of upstream cracks is shown in Figure20. Dam heel cracking occurs as the overload increases to 1.7~3.0P0. There is no sign of existing cracks growth in the dam during the test. Instead, cracks that occur on the dam surface begin to extend after the overload reaches 4.0P0. Experimental results indicate that dam heel cracking, compared with existing cracks in dam, is the dominating problem of Xiaowan arch dam, which is corresponding to numerical results.Figure 20
The final distribution of upstream cracks.


## 4.1. Numerical Model

Xiaowan arch dam, with a height of 292 m, is subjected to some cracks in dam due to the temperature control program. Most of the cracks distribute on the middle and bottom elevations of the dam and tend tangentially to the arch. The average width of these cracks is 1 mm. Distribution of cracks on typical elevation and on typical dam section is shown in Figures8 and 9, respectively.Figure 8
Distribution of cracks on typical elevation.Figure 9
Distribution of cracks in typical dam section.The size of numerical model is as follows: upstream 300 m, downstream 900 m, and left and right banks 700 m each. Elevations of dam crest and model bottom are∇1245 m and ∇595 m, and thus the model height is 650 m. The model includes 58989 nodes and 53850 elements, as shown in Figure 10. There are 11 major cracks that are numbered and simulated, as shown in Figure 11 and Table 2.Table 2
Numbers of cracks in the dam.
Number
Crack
1
13-1
2
13-2
3
20-1
4
20-2
5
22-1
6
22-3
7
22-4
8
25-1
9
28-1
10
28-2
11
30-1
Crack 13-1 in the table means the 1st crack in dam section number 13, and so on.Figure 10
Numerical model of Xiaowan arch dam.Figure 11
Distribution of cracks in the dam.The dam cracks are induced by the improper temperature control program during construction. Since concrete grouting in the cracks is operated, the dam cracks are simulated with thin layer elements. The material parameters are shown in Table3.Table 3
Material parameters.
E/GPa
μ
r/kN·m−3
f
c/MPa
Foundation
18
0.28
28
1.52
1.33
Dam
23.1
0.2
24
1.4
1.6
Crack
5
0.3
27.5
1.12
0.9


## 4.2. Fracture Analysis

Structural fracture and parameters sensitivity analysis are presented. Three schemes are performed, including normal parameters (scheme 1), crack’s shear strengthsfand c reduced by 50% (scheme 2), and crack’s shear strengths fand c reduced by 75% (scheme 3). Both the hydrostatic pressure P0 and overload conditions are included in each scheme.


### 4.2.1. Dam Heel Cracking Analysis

The distribution of unbalanced force at dam heel under different loads in scheme 1 is shown in Figure12. The other two schemes are of similar distribution. Under the hydrostatic pressure P0, unbalanced force rarely occurs and dam heel stays as a global steady state. As the hydrostatic pressure is overloaded, unbalanced force area expands and its magnitude increases rapidly, which is concentrated in the dam heel.Distribution of dam heel unbalanced force of scheme 1.
(a)
(b)
(c)Unbalanced forces of dam heel in three schemes are summarized in Table4. In the overload conditions, unbalanced force increases rapidly in dam heel, which is the most possible location of cracks initiation. As the variation of unbalanced forces between the three schemes is relatively small, the dam heel cracking behavior is insensitive to the change of crack parameters.Table 4
Unbalanced forces of dam heel (104 N).
Scheme
P
0
1.5P0
2.0P0
1
40249.56
257289.9
415601.5
2
40216.39
257269.2
415560.7
3
40126.59
256869.2
414973.0


### 4.2.2. Dam Body Cracking Analysis

The distribution of unbalanced force in dam under different loads in scheme 1 is shown in Figure13. Under the hydrostatic pressure P0, there is a little unbalanced force existing in the right abutment of the arch dam. As the pressure is overloaded, unbalanced force occurs in both the dam abutments. Distribution of unbalanced force in dam is close to the foundation, where cracks initiate and failure may occur.Distribution of unbalanced force in dam of scheme 1.
(a)
(b)
(c)Distribution of dam body unbalanced force under overload2.0P0 in schemes 2 and  3 is shown in Figures 14 and 15, respectively. Dam unbalanced force is insensitive to the change of crack parameters. Crack number 10, that is, the second crack in dam section number 28, suffers from some unbalanced force in scheme 3, which indicates the most sensitive and dangerous crack in the dam. Meanwhile, the unbalanced force of cracks is much less than that of dam body, and the latter is still the emphasis of dam cracking prevention.Figure 14
Distribution of dam body unbalanced force under overload2.0P0 in scheme 2.Figure 15
Distribution of dam body unbalanced force under overload2.0P0 in scheme 3.


### 4.2.3. Analysis for Existing Cracks

The length of unbalanced force vectors in the cracks is enlarged 20 times since it is relatively less than that of dam body. Distribution of cracks unbalanced force under overload2.0P0 in the three schemes is shown in Figures 16, 17, and 18, respectively. Cracks initiation and propagation are very sensitive to the crack parameters. In other words, the stability of the existing cracks relies on the quality of concrete grouting in the cracks.Figure 16
Distribution of cracks unbalanced force under overload2.0P0 in scheme 1.Figure 17
Distribution of cracks unbalanced force under overload2.0P0 in scheme 2.Figure 18
Distribution of cracks unbalanced force under overload2.0P0 in scheme 3.Unbalanced force of cracks under overload2.0P0 in all three schemes is summarized in Table 5. Unbalanced force in dam heel is also listed in Table 5 as comparison. Unbalanced forces in dam heel increase earlier than the cracks in the dam. Dam heel contributes most of the unbalanced forces in all conditions. Namely, dam heel cracking occurs before any existing crack propagates. Among all existing cracks, 20-2 and 28-2 are the dominating cracks in the process of fracture propagation, while 30-1 and 13-2 are also possible initiation area, as shown in Figure 19. The cracks near the crown cantilever are stable.Table 5
Unbalanced force of cracks under overload2P0 (104 N).
Number
Crack
Scheme   1
Scheme   2
Scheme   3
1
13-1
0.00
0.00
12.54
2
13-2
0.02
0.87
221.69
3
20-1
0.00
0.00
2.83
4
20-2
0.03
0.70
7614.90
*
5
22-1
0.00
0.00
0.00
6
22-3
0.00
0.00
0.00
7
22-4
0.00
0.00
0.00
8
25-1
0.00
0.00
8.00
9
28-1
0.00
0.94
65.95
10
28-2
1.77
345.54
15810.91
*
11
30-1
0.02
6.04
687.85
Dam heel
415601.5
415560.7
414973
*Dominating cracks in the propagation process.Figure 19
Local distribution of cracks unbalanced force (scheme 1, under overload2.0P0).


## 4.2.1. Dam Heel Cracking Analysis

The distribution of unbalanced force at dam heel under different loads in scheme 1 is shown in Figure12. The other two schemes are of similar distribution. Under the hydrostatic pressure P0, unbalanced force rarely occurs and dam heel stays as a global steady state. As the hydrostatic pressure is overloaded, unbalanced force area expands and its magnitude increases rapidly, which is concentrated in the dam heel.Distribution of dam heel unbalanced force of scheme 1.
(a)
(b)
(c)Unbalanced forces of dam heel in three schemes are summarized in Table4. In the overload conditions, unbalanced force increases rapidly in dam heel, which is the most possible location of cracks initiation. As the variation of unbalanced forces between the three schemes is relatively small, the dam heel cracking behavior is insensitive to the change of crack parameters.Table 4
Unbalanced forces of dam heel (104 N).
Scheme
P
0
1.5P0
2.0P0
1
40249.56
257289.9
415601.5
2
40216.39
257269.2
415560.7
3
40126.59
256869.2
414973.0


## 4.2.2. Dam Body Cracking Analysis

The distribution of unbalanced force in dam under different loads in scheme 1 is shown in Figure13. Under the hydrostatic pressure P0, there is a little unbalanced force existing in the right abutment of the arch dam. As the pressure is overloaded, unbalanced force occurs in both the dam abutments. Distribution of unbalanced force in dam is close to the foundation, where cracks initiate and failure may occur.Distribution of unbalanced force in dam of scheme 1.
(a)
(b)
(c)Distribution of dam body unbalanced force under overload2.0P0 in schemes 2 and  3 is shown in Figures 14 and 15, respectively. Dam unbalanced force is insensitive to the change of crack parameters. Crack number 10, that is, the second crack in dam section number 28, suffers from some unbalanced force in scheme 3, which indicates the most sensitive and dangerous crack in the dam. Meanwhile, the unbalanced force of cracks is much less than that of dam body, and the latter is still the emphasis of dam cracking prevention.Figure 14
Distribution of dam body unbalanced force under overload2.0P0 in scheme 2.Figure 15
Distribution of dam body unbalanced force under overload2.0P0 in scheme 3.


## 4.2.3. Analysis for Existing Cracks

The length of unbalanced force vectors in the cracks is enlarged 20 times since it is relatively less than that of dam body. Distribution of cracks unbalanced force under overload2.0P0 in the three schemes is shown in Figures 16, 17, and 18, respectively. Cracks initiation and propagation are very sensitive to the crack parameters. In other words, the stability of the existing cracks relies on the quality of concrete grouting in the cracks.Figure 16
Distribution of cracks unbalanced force under overload2.0P0 in scheme 1.Figure 17
Distribution of cracks unbalanced force under overload2.0P0 in scheme 2.Figure 18
Distribution of cracks unbalanced force under overload2.0P0 in scheme 3.Unbalanced force of cracks under overload2.0P0 in all three schemes is summarized in Table 5. Unbalanced force in dam heel is also listed in Table 5 as comparison. Unbalanced forces in dam heel increase earlier than the cracks in the dam. Dam heel contributes most of the unbalanced forces in all conditions. Namely, dam heel cracking occurs before any existing crack propagates. Among all existing cracks, 20-2 and 28-2 are the dominating cracks in the process of fracture propagation, while 30-1 and 13-2 are also possible initiation area, as shown in Figure 19. The cracks near the crown cantilever are stable.Table 5
Unbalanced force of cracks under overload2P0 (104 N).
Number
Crack
Scheme   1
Scheme   2
Scheme   3
1
13-1
0.00
0.00
12.54
2
13-2
0.02
0.87
221.69
3
20-1
0.00
0.00
2.83
4
20-2
0.03
0.70
7614.90
*
5
22-1
0.00
0.00
0.00
6
22-3
0.00
0.00
0.00
7
22-4
0.00
0.00
0.00
8
25-1
0.00
0.00
8.00
9
28-1
0.00
0.94
65.95
10
28-2
1.77
345.54
15810.91
*
11
30-1
0.02
6.04
687.85
Dam heel
415601.5
415560.7
414973
*Dominating cracks in the propagation process.Figure 19
Local distribution of cracks unbalanced force (scheme 1, under overload2.0P0).


## 4.3. Comparison with Geomechanical Model Test

Results and conclusions of geomechanical model test are presented as comparison with numerical method. The final distribution of upstream cracks is shown in Figure20. Dam heel cracking occurs as the overload increases to 1.7~3.0P0. There is no sign of existing cracks growth in the dam during the test. Instead, cracks that occur on the dam surface begin to extend after the overload reaches 4.0P0. Experimental results indicate that dam heel cracking, compared with existing cracks in dam, is the dominating problem of Xiaowan arch dam, which is corresponding to numerical results.Figure 20
The final distribution of upstream cracks.


## 5. Conclusion

Unbalanced force is proposed based on deformation reinforcement theory to analyze fracture behavior, including initiation and propagation of cracks in 3D structures. Unbalanced force is a set of the equivalent nodal forces of stress exceeding the yield function, which can be termed the driving force of time-dependent deformation, as well as damage evolution.The unbalanced force and damaged area are in good agreement with precrack specimen test results. The distribution of unbalanced force indicates cracks initiation area, while its direction predicts the possible cracks propagation path.The method is applied in fracture analysis of Xiaowan high arch dam. Dam heel cracking occurs before any existing crack propagates, which is the most possible failure mode. Among all existing cracks, 20-2 and 28-2 are the dominating cracks in the process of fracture propagation.



---

*Source: 658160-2013-10-28.xml*