A Multiparameter Damage Constitutive Model for Rock Based on Separation of Tension and Shear

Abstract

By analysis of the microscopic damage mechanism of rock, a multiparameter elastoplastic damage constitutive model which considers damage mechanism of tension and shear is established. A revised general form of elastoplastic damage model containing damage internal variable of tensor form is derived by considering the hypothesis that damage strain is induced by the degeneration of elastic modulus. With decomposition of plastic strain introduced, the forms of tension damage variable and shear damage variable are derived, based on which effects of tension and shear damage on material’s stiffness and strength are considered simultaneously. Through the utilizing of Zienkiewicz-Pande criterion with tension limit, the specific form of the multiparameter damage model is derived. Numerical experiments show that the established model can simulate damage behavior of rock effectively.

1. Introduction

Rock is a kind of multiphase and inhomogeneous material, with mesoscopic discontinuities randomly distributed. When subjected to loads, discontinuities emerge and develop, leading to the degeneration of material’s strength and stiffness. The exact simulation of the damage behavior of rock requires definition of damage variables based on statistic methods and determination of evolution of those damage variables according to specific physical background so as to establish models like system of microstructures [1]. Owing to the complexity of mesoscopic structure of rock, the establishment and application of mesoscopic damage model can be difficult and time consuming. Constructing a macroscopic damage model based on continuum hypothesis will be beneficial to the simplification of problem and application in engineering.

Many works have been done on damage property of concrete and geomaterials. Salari et al. [2] established a triaxial damage model of geomaterials accounting for tensile damage. Nguyen and Korsunsky [3] established an isotropic damage model for concrete which addresses the relation between local and nonlocal parameters. Li et al. [4] introduced statistical method to describe the strain softening behavior of rock. For materials like rock, the basic damage types in macroscopic view include tension, shear, and crush [5], of whom some more complicated damage forms can be viewed as superposition. To describe the damage behavior of these materials, the classical damage model with only one parameter is not enough. Damage internal variable in tensor form with several independent parameters included should be introduced to simulate the mechanical behavior of different damage types. Many works have been done along this approach. Frantziskonis and Desai [6] established a model suitable for geologic materials by introducing a tensor form damage variable to describe the structural changes in such materials. Krajcinovic and Mastilovic [7] considered scalar, second-, fourth-, and sixth-order tensor representations of damage and evaluated the accuracy with which they approximate exact, micromechanical solutions. To simulate the damage behavior of rock, it is needed to relate the damage internal variable to each damage form, thus allowing each part of damage internal variable representing different damage types to evolve according to specific laws. Resende [8] suggested a rate-independent constitutive theory for concrete which consider shear damage and hydrostatic tension damage; this model divides stress into mean value part and deviatoric part and considers the damage criteria and evolution laws, respectively. Jirásek [9, 10] studied nonlocal constitutive models for damage and fracture processes of quasibrittle materials and explored a formulation with averaging of the displacement field. Lee and Fenves [11] establish a plastic damage model for concrete subjected to cyclic loading using the concepts of fracture-energy-based damage and stiffness degradation. Li and Wu [12, 13] follow the approach of stress splitting and derived a damage constitutive model for concrete in effective stress space based on energy principal. Ju [14] established an energy-based coupled elastoplastic damage theory by considering Helmholtz free energy featuring strain-split formulation which leads to more robust algorithm.

Many works on damage model are dedicated towards the notion of effective quantities, making the establishment of damage state relatively indirect. In this paper, we follow the approach of taking damage internal variable as a two-order three-dimension tensor but adopt the viewpoint raised by Ying [15] that plastic deformation and damage of rock can be treated uniformly as a macroscopic representation of the emergence, development, and accumulation of microscopic faults within rock. In this way, a revised general form of elastoplastic damage model for rock containing damage variable of tensor form is derived. By utilizing the assumption of strain separation, the tensor form damage variable is related to different damage forms in an intuitive way and the evolution laws of strength and stiffness for rock materials in these damage states are described, respectively.

2. A Revised General Form of Elastoplastic Model with Damage Variables of Tensor Form Considered

Continuum damage models are usually established by introducing particular damage variables and evolution laws, where the damage variables show to what extent damage develops. Nowadays, what is wildly used in a large amount of literature is the concept of effective stress, in which, take the scalar damage variable , for example, takes value 0 or 1 denoting no damage or full damage, respectively, and can be seen as the damage coefficient; a model containing damage variable can be established through the principal of equivalence (stress equivalence, strain equivalence, or energy equivalence). However, when considered based on the intrinsic mechanism, damage can, together with plastic deformation, be seen as the result of emergence and development of microscopic faults that only differs from plasticity by the effects on strength and stiffness [15]. In this section, the results in reference [15] are generalized and derivation similar to plastic internal variables is made about damage internal variables.

As exposed in previous section, rock damage has several basic types, each affecting the strength and stiffness in different way. Similar to plastic internal variables, introduce in constitutive relation a damage internal variable    of tensor form, which leads to where and denote stress and strain components, respectively, and denotes plastic internal variable. Differentiate (1) with respect to time; stress rate can be expressed as where , , and each denotes elastic stress rate, plastic stress rate, and damage stress rate (Einstein’s summation convention is always assumed in this paper); here, denotes elastic stiffness. For static problem, there is no need to differentiate with respect to real time. Similarly, strain rate can be expressed as where , , , and is the elastic flexibility.

Assume and are the bound for plasticity and damage in stress and strain space. In fact, shear damage criterion for rock material is often expressed as a form similar to plastic yield function (such as functions of Drucker-Prager type) in application; this justifies the unified treatment for the judgment of plasticity and damage. For materials satisfying Ilyushin’s postulate , similar to the general theory for elastoplastic model, a generalized associated flow law can be obtained by constructing a circular path in strain space: Equation (4) is also suitable for softening material. Similarly, (4) can be extended to nonassociated materials. Still treat damage and plasticity in a unified way, and it can be assumed that or where and denote plastic and damage potential function in stress and strain space, respectively.

From the consistency condition in strain space that , it can be derived that is proportional to when loading; namely, , where is the function of strain and . More generally, to take into account situations when it is not loading, introduce Macaulay bracket satisfying it can be obtained that Combine (6) and (8); it can be shown that has direction and is proportional to . Thus, it can be assumed that and similarly for stress rate, Let and . When loading, it can be obtained that or . Suppose and satisfy the relation ; it can be derived that , where . When , we have . Let ; it can be obtained that and together with (9) and (10) we get

Equations (11) and (12) have the same looking with general elastoplastic relation but implicitly contain the effect of damage variable on rock’s strength and stiffness. To clarify it, the specific form of is derived. Introduce two hypothesizes raised in [15]:(1)Rock’s stiffness matrix degenerates as damage develops. Damage strain is induced by the degeneration of elastic stiffness matrix coefficients; namely, . Thus, damage strain rate can be defined as (2) can be expressed as a homogeneous linear form of ; that is,  where is the function of plastic strain and damage variable.Substituting (13) and (14) into (12) leads to The equation above can be regarded as a system of linear equations about plastic strain rates, which can be solved as where

Substitute (16) into the consistency condition in stress space: Also consider that plastic internal variables are usually assumed as functions of plastic strain . It can be derived that Specifically, suppose plastic internal variable satisfies ; it can be obtained that where . In conclusion, hypothesize (1) implies the effect of damage on material’s stiffness while (19) implies the effects on material’s strength (plasticity and damage bound surface).

Consider the Clausius-Duhem inequality: where denotes density, denotes spatial coordinates, denotes temperature, denotes heat flow, and and denote free energy and entropy per unit mass. Similar to stress, can be expressed as and thus we have Substituting (22) into (21) leads to Suppose the Coleman relation is satisfied; inequality (23) is valid if the inequalities are satisfied, which represent mechanical and thermal dissipation, respectively.

3. The Form of Damage Variables Based on Separation of Tension and Shear

Tension and shear are the main types of damage for rock. Simo and Ju [14, 16] separate stress (strain) tensor into two different parts and consider the damage effect, respectively, according to the assumption of separation of tension and shear. Follow this approach and omit the effect of crush; damage variable can be reduced to a simpler form relating to tension and shear independently. Let in (14) be where Here, denotes cosine of the angle between the th principal direction of plastic strain and each coordinate axis and , are dimensionless function about plastic strain invariables. Equation (25) can be seen as a separation of damage variable into the tension and shear part. In fact, the principal form of can be obtained when applied with which deviates according to the sign of its principal value. Here, and in (26) can be seen as cut-off functions that reserve only positive and negative part of strain, respectively; thus, separating damage variable into the tension and shear part, each part casts back to the component form.

Substitute (25) into (14) and let ; it can be obtained that Let To simplify derivation, omit derivatives of principal direction with respect to time; we have Damage internal variable is thus reduced to three variables   , each relating to the three plastic principal strains, respectively. It is now possible to replace the damage variable with its simplified form . By using (28), we have . The damage part in (19) can thus be expressed as

4. Multiparameter Damage Model Based on Zienkiewicz-Pande Criterion

To expose the model further, in this section, Zienkiewicz-Pande criterion with tension limit is taken as an example and the specific form of the multiparameter damage model is derived. The criteria currently used very often in research and engineering include Drucker-Prager, Mohr-Coulomb, Hoek-Brown, and unified strength theory suggested by Yu et al. [17], of which the Mohr-Coulomb criterion is widely adopted. The existence of vertex singularity in bound surface will bring some difficulties in iteration. Zienkiewicz-Pande criterion is a smooth version of Mohr-Coulomb criterion [18], which uses a quadratic curve to approximate straight line in -plane so as to eliminate singular vertex and benefit numerical implementation; also it takes intermediate principal stress into account to some extent. The Zienkiewicz-Pande criterion has been used in many engineering projects for the simulation like excavation of underground caverns and shows a good agreement with monitoring results [19].

4.1. The Specific Form of Multiparameter Damage Model

In damage judgment, plastic yield surface is often used as the bound for shear damage and stress or strain limit as the bound for tension. Thus, in this section, Zienkiewicz-Pande criterion with tension limit is assumed as the bound for plasticity and damage; namely, Here, the hyperbolic type of approximation in [20] is adopted; that is, where is the cohesion, is the angle of internal friction, is the tension limit, is the mean stress , is lode angle, is lode radius, namely, , is the revising coefficient introduced in [20], and controls to what extent approximation can be made to the Mohr-Coulomb yield surface. Here, the damage type is determined by the first bound the state reaches when loading. Let potential function be where Here, denotes the angle of dilatancy. Consider that , , and , , and ; it can be derived that where Consider that ; it is easy to get accordingly. Similarly, and can be obtained.

Let elastic flexibility be an isotropic tensor: where denotes Young’s modulus and denotes Poisson’s ratio. Suppose only Young’s modulus is affected by damage internal variables; then, similar to (30), we can get where Substitute (38a) and (38b) into (17); we can get , which together with (20) and (35a) and (35b) (also , , and which can be obtained in a similar approach) leads to a constitutive relation in the form as (11) and (12).

4.2. The Damage Evolution of Material Strength and Stiffness

Material strength and stiffness degenerate as damage evolves. Here, those strength and stiffness parameters are treated as functions of the reduced form of damage variables and their specific forms are examined. Consider that damage variable has a degenerating effect on cohesion and angle of internal friction; it can be derived that Similarly, consider the degenerating effect damage variable has on tension limit; we have Since Zienkiewicz-Pande yield surface is a smooth approximation to the Mohr-Coulomb yield surface and the Mohr-Coulomb criterion has the form during uniaxial compression, it can be assumed that the cohesion in Zienkiewicz-Pande criterion is proportional to uniaxial compressive strength . Instead of plastic hardening effect, damage makes material’s strength degenerate. To simplify implementation, omit the effect damage variable has on angle of internal friction, and also consider the function type adopted in [21] for the degeneration of material’s parameters; let and be linear functions for damage variables: where and denote the original compression and tension strength and and are damage modulus for shear and tension, respectively.

As internal microscopic cracks develop very fast after reaching damage [22], exponential function is adopted: where is the original Young’s modulus, is material’s damage coefficient, and is taken as the linear combination of damage variables . Specifically, let Given the irreversibility of damage during the damage evolution of material’s parameters, the minima of material parameters along the evolution path are used.

5. Numerical Example

The algorithm for solving group of nonlinear equations includes Newton-Raphson method, modified Newton-Raphson method, and quasi-Newton-Raphson method, of which Newton-Raphson method and modified Newton-Raphson method can effectively utilize the sparsity of stiffness matrix. Consider that the model raised in this paper has many factors interacting with each other and thus highly nonlinearized; a framework for nonlinear finite element analysis by using Newton-Raphson method is implemented, which calls specific constitutive model to handle iterative solving.

Typical implementation of material model can be reduced to the following problems:(1)The establishment of tangent stiffness matrix, which can be implemented using (11).(2)State update, namely, update of stress, strain, and material parameters of each quadrature point when strain increment is known. Typical implementation includes explicit method (e.g., Runge-Kutta method) and implicit method (e.g., return map method [23, 24]). To simplify implementation, explicit method is adopted which takes the incremental form of (11) and perform update and revision for state and material parameters at the end of each increment step.

To validate the effectiveness of the model established, consider the cylindrical sample in [25]. Let  GPa, ,  MPa, and . Apply uniaxial compression to the sample model. The stress-strain curve is depicted in Figure 1. Comparison with the experiment result of uniaxial compression in [25] shows that the established model can simulate this process well.

Figure 1 

Stress-strain curve under uniaxial compression (MPa).

Consider a cubic finite element model that each side has length 1 m. Let  GPa, ,  MPa, , and  MPa. Apply cyclic compression and tension uniaxial load. The stress-strain curves are shown in Figure 2. When strain reaches 0.0018 and 0.0028 for compression or when strain reaches 0.00016 and 0.0022 for tension, unload stress to zero and then apply again. The tangent stiffness during unloading is smaller than elastic stiffness as the introduction of damage. The model presented in this paper can simulate damage behavior under cyclic compression and tension effectively.

(a) Compression

(b) Tension

Figure 2 

Stress-strain curve under cyclic uniaxial compression and tension (MPa).

6. Engineering Application

In this section, the established model is used for the excavation simulation of Zhen’an underground power station.

6.1. Project Overview and Numerical Modeling

Zhen’an pumped storage power station is located in Shanxi Province in China with installed capacity of 1400 MW (4 × 350 MW). The underground powerhouse is made up of the main powerhouse, main transform house, tailrace surge chamber, and several other tunnels connecting each cavern. Mechanical parameters for rock are shown in Table 1. Finite element model shown in Figure 3 contains 480608 isoparametric 8-node elements, with 109011 elements to be excavated divided into 10 stages. -axis of the model is perpendicular to longitudinal axis of main powerhouse with a range of 318.64 m; -axis is along the longitudinal axis with a range of 389.6 m; -axis is along the vertical axis with a range of 568.6 m. Initial stress is shown in Figure 4. It is obtained by inverse analysis based on 4 gauging points near the carven and the result shows that it is generated mainly by gravity. The established model is used for excavation simulation of the underground powerhouse.

Table 1 

Mechanical parameters for rock.

(a)

(b)

Figure 3 

Computation model (a) and excavation elements (b).

(a)

(b)

Figure 4 

Initial first (a) and third (b) principal stress (MPa).

6.2. Result Analysis

6.2.1. Damage Zone Distribution

To simulate the process of excavation, divide the computation into 10 steps, each with a layer shown in Figure 3 removed, and let the iteration reach convergence from the unbalanced state. Figure 5 gives the development of damage zone where the third and sixth stage are the excavation of lateral caverns, resulting in some finite elements returning from damage bound. When excavation is complete, damage zone distribution of 2# unit section’s surrounding rocks is shown in Figure 6. Depth of damage zone reaches a maximum of 2.1 m near the arch roof. Depth of damage zone near side wall ranges from 18.4 m to 20.1 m and ranges from 7.0 m to 9.0 m near main transform house. Tension damage zone is mostly distributed near the intersection of tunnels and caverns, where release of initial stress along multiple direction makes it easy to crack.

Figure 5 

Volume of damage zone (10⁵ m³).

Figure 6 

Damage zone of 2# unit section.

6.2.2. Stress Distribution

After excavation stress, distributions of surrounding rock of each unit section are roughly similar. The first and third principal stress of 2# unit section’s surrounding rocks are shown in Figure 7; the direction of the third principal stress is almost vertical, which is the same as the direction of gravity, and first principal stress is along the radial direction. Near the intersection of each tunnel with the main powerhouse, the third principal stress concentrates and reaches a maximum of −20 MPa. First principal stress in the middle part of 2# unit section’s side wall reaches a peak of −23.7 MPa and third principal stress is around −2.0 MPa.

(a)

(b)

Figure 7 

First and third principal stress of 2# unit section (MPa).

6.2.3. Displacement Distribution

What is shown in Figure 8 is the displacement of 2# unit section. Release of initial stress makes the displacement directing towards the caverns. Displacement near the intersection of each cavern is relatively larger than that of other areas. In other areas, displacement field is distributed relatively uniformly. After excavation, 2# unit section’s arch roof’s displacement resiles to 3.0 mm. Maximum displacement of main powerhouse’s upstream and downstream side wall is 42.1 mm and 35.5 mm, respectively. The values are 13.8 mm and 19.6 mm for main transform house.

Figure 8 

Displacement of 2# unit section (mm).

From the analysis above, it is easy to know that the damage zone, stress, and displacement distribution obtained using the model established in this paper are similar to those using Zienkiewicz-Pande or Mohr-Coulomb criterion because of the form of shear damage bound adopted. The introduction of damage mechanism makes it more natural to figure out the shear and tension damage area.

7. Conclusion

This paper approaches based on the microscopic mechanism of damage for rock material and establishes in an intuitive way a multiparameter elastoplastic damage model for rock that is applicable to engineering. With the assumption that damage comes into existence as the material’s strength and stiffness degenerate and that damage is interconnected with plastic deformation, a revised general form for elastoplastic damage model containing damage variable of tensor form is established. By considering plastic strain separation, the expression of damage variable reflecting the damage mechanism for shear and tension simultaneously is derived. By adopting Zienkiewicz-Pande criterion with tension limit as the bound for plasticity and damage, the specific form for the damage model is derived and implemented.

The model established in this paper is physically intuitive and has relatively well-based theoretical background. Numerical experiments and engineering application show that this model can reflect the damage behavior of rock effectively.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.