**1. Introduction**

Fiber-reinforced concrete (FRC) has been widely used in the field of civil engineering for its excellent physical and mechanical properties. FRC is a kind of cement-based composite material composed of metal fiber, inorganic non-metallic fiber, synthetic fiber, or natural organic fiber as a reinforcing material. The most widely used are steel fiber-reinforced concrete (SFRC) and polypropylene fiber-reinforced concrete (PFRC). A large number of experimental studies have shown that these randomly distributed fibers can effectively hinder the expansion of micro-cracks in the concrete and the formation of macro-cracks, significantly improving the compressive, tensile, bending, impact resistance, and fatigue resistance of concrete [1–12]. Fiber materials are further used in the research of new types of concrete, such as coral concrete, geopolymer concrete, self-compacting concrete, etc. Liu et al. [7] carried out uniaxial compression tests of carbon fiber-reinforced coral concrete and established empirical piecewise constitutive model. At present, the representative theories on the mechanism of fiber reinforcement mainly include the fiber spacing theory proposed by Romualdi and Batson [13], and the reinforcement rules for composite materials proposed by Swamy [14].

The constitutive relation of concrete material is one of the key problems in the nonlinear analysis of concrete structure. Much effort has been devoted in the last decade to model the FRC by considering it as a composite material composed of concrete matrix

**Citation:** Bai, W.; Lu, X.; Guan, J.; Huang, S.; Yuan, C.; Xu, C. Stress–Strain Behavior of FRC in Uniaxial Tension Based on Mesoscopic Damage Model. *Crystals* **2021**, *11*, 689. https://doi.org/ 10.3390/cryst11060689

Academic Editors: Yifeng Ling, Chuanqing Fu, Peng Zhang, Peter Taylor and Nima Farzadnia

Received: 23 May 2021 Accepted: 12 June 2021 Published: 16 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and fibers [15–19]. Uniaxial stress–strain behavior is the most fundamental constitutive relationship of the concrete material. Due to the limitation of experimental technology and the deficiency of relevant theories, the research achievements on uniaxial compression are abundant, while the research on tension is relatively few. The research on the constitutive relation of FRC is based primarily on the relevant theoretical results of ordinary concrete, mainly focusing on the influence of fiber characteristic parameters on the constitutive relation. In recent years, based on the direct tensile test carried out on the universal test machine, the uniaxial tensile stress–strain curves for different types of FRC were obtained, and the corresponding macroscopic phenomenological constitutive models were proposed [3,8,20–22].

As a typical quasi-brittle material, the nonlinear macroscopic stress–strain behavior of concrete are closely related to the heterogeneity in microstructure. The essential nature lies not only in the initiation and propagation of individual micro-cracks, but also in the interaction and coalescence of crack populations [23,24]. Most of the existing concrete constitutive relations are macroscopic phenomenological models established by polynomial and other mathematical formulas, focusing on the fitting of test data, but lacking the clear physical meaning for the parameters. Therefore, it is difficult to reveal the meso-damage evolution law of concrete in essence by those phenomenological models.

Damage Mechanics (DM) [25–27] is a relatively new field studying the response and reliability of materials with countless randomly distributed irregular microcracks. The fundamental aspect of damage mechanics is the selection of damage parameters. The essential feature of the original Kachanov's model [28] resides in the introduction of a special internal variable defined by the state of local damage and its accumulation. Many macroscopic continuum damage models and micromechanics damage models [25,29] were proposed after Kachanov's work. Li et al. [21] suggested a continuum damage mechanicsbased model for FRC in tension, in which the quasi-brittleness of the matrix and the fiber– matrix interfacial properties were taken into consideration. Considering the stochastic characteristics of the distribution of disfigurements in the microstructure, a physically motivated damage model (i.e., the bundle parallel bar system (PBS)) was proposed by Krajcinovic and Silva [30], based on the continuous damage theory and statistical strength theory. After this work, Statistical Damage Mechanics (SDM) was gradually formed, constitutive relations of quasi-brittle material (rock and concrete) have been extensively studied, and some inspiring results were obtained [31–34]. This kind of model abstracts quasi-brittle material into a complex system composed of mesoscopic physical elements. The heterogeneity in the microstructure is introduced by assuming that the characteristic parameters of the mesoscopic unit follow Weibull or other statistical distribution forms. The progressive damage accumulation process of concrete is described by means of probability and statistics. It ignores the complicated physical details of the damage process and avoids the complicated calculation of statistical mechanics. The relationship between meso-damage mechanism and macroscopic mechanical behavior of material is effectively established. Considering there are two fundamental damage modes (fracture and yield) in the microstructure of concrete, the statistical damage models of concrete under uniaxial and multiaxial loading were proposed by Chen et al. [35] and Bai et al. [36–38], which could be used to predict the constitutive behavior of concrete under a complex loading environment and explore the damage mechanism of the concrete material.

In this paper, a statistical damage constitutive model for the FRC is presented based on the macroscopic test phenomenon and the statistical damage theory. Firstly, a deep analysis of the mesoscopic damage mechanism for concrete material under uniaxial tension is elaborated. It indicated that the macroscopic nonlinear stress–strain behavior is determined by the evolution of fracture and yield damage on a meso scale. The yield damage mode, reflecting load redistribution and adjustment of stress skeleton in the microstructure, is emphasized as the key role in the whole deformation and failure process. Evolutionary factor is introduced to reflect the development of the potential mechanical capacity of materials. Subsequently, the triangular probability distributions with four parameters are used

to simulate the cumulative evolution of fracture and yield damage. By defining the above parameters as a function of fiber content, the impact of fiber content on the nucleation and growth of microcracks and the mesoscopic damage evolution of the FRC is reflected. The determination method for the model's characteristic parameters is proposed. The validity of the proposed model is demonstrated by comparing the theoretical and experimental results. The influence rule of fiber content on the mesoscopic damage evolution of the FRC is analyzed.

### **2. Deformation and Failure of Concrete under Uniaxial Tension**

Since the failure of concrete is essentially caused by the nucleation and growth of microcracks produced by local tensile strain, the uniaxial tension could be regarded as the most fundamental failure form for concrete.

According to the macroscopic experimental phenomena, the deformation and failure of the concrete specimen under uniaxial tension seems to be divided into two stages, the uniform damage stage (US) and the local fracture stage (LS), as shown in Figure 1a. In the uniform damage stage, the whole specimen remains uniformly loaded and deformed, with the damage evolving mainly by the nucleation and growth of micro-defects, which randomly distribute in the whole specimen. In the local fracture stage, a macroscopic main-crack perpendicular to the tensile direction forms in the fracture process zone (FPZ). The macroscopic response of the specimen strongly depends on the size of the largest crack with a preferential orientation in the FPZ. Local further damage and fracture occurs in the FPZ, while the rest of the region of the specimen shows the unloading behavior.

**Figure 1.** Two-stage feature of macroscopic mechanical behavior of concrete under uniaxial tension: (**a**) deformation-failure process; (**b**) typical nominal stress–strain curve.

As shown in Figure 1b, a typical master stress–strain curve of concrete specimen under uniaxial tension is presented. The signs A and B denote, respectively, the two characteristic states, i.e., (1) the peak nominal stress state and (2) the critical state at which macro-crack starts to develop in the FPZ and the damage process switches into the local failure stage. From the shape of the curve, the stress–strain curve seems to be divided into a stress-strengthening phase and a strain-softening phase bounded by state A. Before A, the tensile stress increases gradually with the growth of tensile strain. After A, the tensile stress decreases gradually with the increase of tensile strain, until to zero. The stress value corresponding to state A, known as the so-called tensile strength, is considered to be one of the most important mechanical indexes for concrete. The descending segment of the tensile stress–strain curve for concrete was firstly measured by Rusch and Hilsdorf [39]. Hughes

and Chapman [40] confirmed that the concrete did not break at the maximum load, but the softening phenomenon occurred.

From the deformation and failure features, and also the process of microcracks propagation, the stress–strain behavior seems more suitable to be divided into an even damage phase and a local fracture phase bounded by state B, corresponding to the two-stage features in Figure 1a. During the whole process, the nucleation and growth of the random distributed microcracks typically ranges from 0.01 to 1.0 mm in width, leading to a concentration of these microcracks into a narrow zone and producing a visible macroscopic fissure wider than 1.0 mm [41,42]. Before B, the material response could be considered as statistical homogeneous. The density of micro-defects is modest, retaining a dilute degree. After this threshold, the macroscopic response of the specimen strongly depends on the macroscopic main-crack in the FPZ. Local softening and fracture (decrease of the nominal stress with increase of the strain) occurs in the FPZ, while the rest of the region of the specimen enters the unloading process. This portion of the stress–strain curve becomes heavily dependent on the measure gauge length, which could not be treated as a pure material mechanical behavior. Based on the catastrophe theory, the process from damage to fracture for quasi-brittle materials is divided into two phases, globally stable (GS) mode (relevant to the distributed damage accumulation) and evolution induced catastrophic (EIC) mode [43]. The critical state transforming from GS to EIC plays a key role during the whole process, and exhibits the critical sensitivity, with many physics showing abnormal behavior. Experiments on rocks have also shown that there exists a critical value for the fracture of the quasi-brittle material [24,44].

For the location of the critical state B, many experiments [45,46] showed that it lies in the softening region behind the peak state A in a tensile stress–strain curve. However, the accurate location of B is still uncertain, and almost all the theoretical analysis in literatures ignored the identification of the critical state, treating A and B simply as the same one.

## **3. Materials and Methods**

### *3.1. Basis of Statistical Damage Theory*

It is well known that statistic physics is a ligament that communicates continuum mechanics, damage mechanics, and material mechanics. Statistical damage theory ignores the microscopic details of damage, and the representative volume element (RVE) is abstracted as a complex system composed of *N* (*N*→∞) mesoscopic units (micro-spring, micro-bar, etc., with the same sectional area *dA* and stiffness *dk*). The heterogeneity of the material is introduced by endowing each unit with different mechanical parameters (strength, characteristic strain, etc.). The macroscopic mechanical properties of concrete depend on the statistical mechanical properties of individual mesoscopic elements. The macroscopic mechanical properties could be described by using a phenomenological model with statistical method.

According to the number of mesoscopic mechanical parameters adopted, statistical damage models could be classified into two types, i.e., the single-parameter model and the double-parameter model, which are represented by the bundle parallel bar system (PBS) [30] and the improved parallel bar system (IPBS) [35,36], respectively.

### 3.1.1. The Series-Parallel Spring Stochastic Damage Model

Based on the PBS, Li and Zhang [47] presented the series-parallel spring stochastic damage model for concrete under uniaxial tensile test. The core idea is to extend the singlelayer parallel element model to multiple layers, which can better simulate the mechanical phenomenon such as localized failure and stress drop. The tensile specimen is abstracted as a complex system composed of n typical units (a unit is composed of numerous microsprings in parallel with equal spacing) in series. *F* is the tensile force. *H* = *nh*<sup>0</sup> is the height of the model, where *h*<sup>0</sup> is the material characteristic height and three times the maximum particle size of aggregate presented by Bažant and Oh [48]. In this model, each micro-spring has the same elastic modulus and cross-sectional area. *ε*Ri is the fracture

strain for micro-spring *i*, as a random variable obeying the same probability distribution for each layer unit, where Weibull distribution and Lognormal distribution are often used.

Since the prominent feature of concrete failure is local failure, the failure of the seriesparallel spring model is caused by the fracture of a unit body, called the main crack unit body (corresponding to the FPZ), and the others as the non-main crack unit bodies. Define a characteristic strain, namely, the critical strain. Before the critical strain, all unit bodies have the same deformation. After the critical strain, the deformation of the main crack unit body continues to increase, while in the non-main crack unit bodies will occur the unloading behavior with the tensile strain no longer increased. Therefore, the average stress–strain curve after the critical state for the concrete tensile specimen shows the obvious size effect, and stress drop will occur when the specimen size is long enough.

### 3.1.2. The Improved Parallel Bar System (IPBS)

Based on the PBS, the IPBS was proposed by Chen et al. [35] to simulate the damage evolution in the FPZ of quasi-brittle material. The fundamental assumptions for the damage mechanism are as follows: (1) the essence of damage and fracture for a quasi-brittle material is the nucleation, growth, and coalescence of micro-cracks, the stress redistribution and adjustment of the force skeleton in the microstructure. (2) The influence of these irreversible micro-changes on macroscopic mechanical properties of material could be addressed by two dominant aspects: decrease of the effective cross-section area and degradation of the elastic modulus corresponding to the effective bearing position. (3) The two damage modes can be simulated, respectively, by rupture and yield of the micro-bar; and the essence of failure can be interpreted as a continuum accumulation and evolution of the two damage modes.

In this model, each micro-bar is composed of spring, cementation bar, and sliding block. The micro-bar *i* has two feature strains, i.e., the fracture strain *ε*Ri and the yield strain *ε*yi, and it may have two kinds of failure modes (elastic-fracture and yield-fracture) according to the values between *ε*Ri and *ε*yi (as shown in Figure 2a). By introducing some simplifying assumptions, the two kinds of damage evolution process in the meso-scale are decoupled (further detailed discussion can be referred to from Chen et al. [35]). As shown in Figure 2b, *q*(*ε*) and *p*(*ε*) are the independent probability density functions to *ε*Ri and *ε*yi, respectively; where *ε*ymin = *ε*Rmin denote the minimums of *ε*yi and *ε*Ri; *ε*ymax < *ε*Rmax are the corresponding maximums. In the course of damage evolution, the cross-sectional area *A*<sup>0</sup> in the initial state will gradually transform into three parts, i.e., *A*1, *A*2, *A*3, denoting the cross-sectional areas corresponding to the fractured bars, the yielded bars, and the bars in elastic state, respectively. *A<sup>E</sup>* is the effective cross-section area bearing load. It satisfies *A*<sup>0</sup> = *A*<sup>1</sup> + *A*<sup>2</sup> + *A*<sup>3</sup> and *A<sup>E</sup>* = *A*<sup>2</sup> + *A*3.

The macroscopic stress–strain curves described by IPBS are drawn in Figure 2c. *σ*<sup>N</sup> is the nominal stress corresponding to *A*0; *σ<sup>E</sup>* is the effective stress corresponding to *AE*; *σ* is the elastic stress corresponding to *A*3. *ε*cr = *ε*ymax; *ε*<sup>u</sup> = *ε*Rmax.

By *ε*ymax, the whole process predicted by the IPBS can be divided into two phases, i.e., partial yield phase (0 ≤ *ε* < *εy*max) and full yield phase (*εy*max ≤ *ε* ≤ *εR*max). When *ε* = *ε*ymax, all the micro-bars in IPBS will yield, and *σ<sup>E</sup>* will reach its maximum. After this state, *σ<sup>E</sup>* will keep the maximum constant with the further increase of *ε*. Therefore, we could use the IPBS to simulate the two-stage deformation and failure characteristics of quasi-brittle materials in the FPZ, by assuming that *ε*ymax corresponds to the critical state B in Figure 1b.

**Figure 2.** Mechanical behavior of concrete under uniaxial tension described by the Improved Parallel Bar System (IPBS): (**a**) micro bar and failure mode; (**b**) damage evolution process on meso-scale; (**c**) stress–strain curves on macro-scale.

The constitutive relation can be expressed as follows:

(1) Partial yield phase (0 < *ε* < *ε*ymax)

$$
\sigma\_{\rm N} = E\_0 (1 - D\_{\rm Y}) (1 - D\_{\rm R}) \varepsilon \tag{1}
$$

$$
\sigma\_{\rm E} = \,^\circ E\_0 (1 - D\_{\rm Y}) \varepsilon \tag{2}
$$

$$D\_{\rm Y} = \int\_0^{\varepsilon} p(\varepsilon) \mathrm{d}\varepsilon - \frac{\int\_0^{\varepsilon} p(\varepsilon) \mathrm{e} \mathrm{d}\varepsilon}{\varepsilon} \tag{3}$$

$$D\_{\mathbb{R}} = \int\_0^\varepsilon q(\varepsilon) d\varepsilon \tag{4}$$

where *D*<sup>y</sup> and *D*<sup>R</sup> denote the accumulated damage variables of elastic modulus of IPBS due to the yield and fracture of the micro-bars; *D*<sup>R</sup> also represents the cumulative distribution function of *q*(*ε*).

(2) Full yield phase (*ε*ymax ≤ *ε* < *ε*Rmax)

$$
\sigma\_{\rm N} = E\_0 (1 - D\_{\rm ymax}) (1 - D\_{\rm R}) \varepsilon\_{\rm ymax} \tag{5}
$$

$$
\sigma\_{\rm E} = \,^\circ E\_0 (1 - D\_{\rm ymax}) \varepsilon\_{\rm ymax} \tag{6}
$$

$$D\_{\text{ymax}} = 1 - \frac{\int\_0^{\varepsilon\_{\text{ymax}}} p(\varepsilon) \varepsilon \text{d}\varepsilon}{\varepsilon\_{\text{ymax}}} = \text{constant} \tag{7}$$

$$D\_{\mathbf{R}} = \int\_0^{\varepsilon\_{\text{ymax}}} q(\varepsilon) \mathbf{d}\varepsilon + \int\_{\varepsilon\_{\text{ymax}}}^{\varepsilon} q(\varepsilon) \mathbf{d}\varepsilon \tag{8}$$

where *D*ymax is the value of *D*<sup>y</sup> corresponding to *ε*ymax.

IPBS also could simulate the unloading process of the rest of the region of the specimen in the local fracture phase. *σ*<sup>N</sup> and *σ*<sup>E</sup> can be expressed by Equations (1) and (2). *D*<sup>y</sup> and *D*<sup>R</sup> can be expressed by Equations (9) and (10).

$$D\_{\rm Y} = \frac{\int\_0^{\frac{\varepsilon\_{\rm{f}} \text{max} - \varepsilon}{2}} p(\varepsilon) \varepsilon \text{d}\varepsilon}{\varepsilon} + \int\_0^{\frac{\varepsilon\_{\rm{f}} \text{max} - \varepsilon}{2}} p(\varepsilon) \text{d}\varepsilon + \frac{\varepsilon\_{\rm{f}\text{max}} \int\_{\frac{\varepsilon\_{\rm{f}} \text{max} - \varepsilon}{2}}^{\varepsilon\_{\rm{f}} \text{max}} p(\varepsilon) \text{d}\varepsilon}{\varepsilon} - \frac{\int\_{\frac{\varepsilon\_{\rm{f}} \text{max} - \varepsilon}{2}}^{\varepsilon\_{\rm{f}} \text{max} - \varepsilon} p(\varepsilon) \text{e}\varepsilon}{\varepsilon} \tag{9}$$

$$D\_{\mathbb{R}} = D\_{\mathbb{R}}(\varepsilon\_{\text{ymax}}) = \int\_0^{\varepsilon\_{\text{ymax}}} q(\varepsilon) d\varepsilon = \text{constant} \tag{10}$$

If the length of specimen (*H*) and the length of the FPZ(*h0*) are given, the average stress–strain curve of the specimen during the local fracture phase could be determined by Equations (5) to (10), which would show the obvious size effect [35].

### 3.1.3. Description of the Damage Evolution Process by IPBS

A typical nominal stress–strain curve and a predicted effective stress–strain curve under uniaxial tension are shown in Figure 3a,b. Signs D and E denote the peak nominal stress state and the critical state, respectively. A denotes the limit elastic state, B and C denote two states in strengthen segment, F and F 0 denote two states in soften segment. The uniaxial tensile damage evolution process of concrete is simulated by using the microscopic spring beam model and the cylinder model, respectively.

In Figure 3c, the micro-bar has two kinds of failure modes, brittle-fracture and yieldfracture. A corresponds to the limit elastic state, where all the micro-bars remain in the elastic state. B, C, D correspond to the two states in strengthen segment and the peak nominal stress state, where some micro-bars are out of work due to fracture, some are yield, and the others still in the elastic state. E corresponds to the critical state, where all the residual bars are yield. F corresponds to the local breach phase in the FPZ, where the yielded bars continue to fracture, exhibiting the softening phenomenon. F 0 corresponds to the local breach phase in the other region, exhibiting the unloading phenomenon.

In Figure 3d, the tensile specimen is divided into three regions, the elastic region, the yield-damage region, and the rupture-damage region, during the whole damage process. With the growth of damage, the elastic region will gradually transform into the other two kinds of regions. When it reaches the critical state E, the elastic region will disappear, the tensile traction will be fully borne by the yield-damage region, and then the damage of specimen will change into the local breach phase.

In Figure 2c, it shows the stress–strain curves predicted by the IPBS, where *σ* is the ideal elastic stress corresponding to the elastic region, *σ*<sup>E</sup> is the effective stress corresponding to *A*<sup>E</sup> (the cross-sectional area relevant to the elastic region and the yield-damage region), *σ*<sup>N</sup> is the nominal stress corresponding to *A*<sup>N</sup> (the initial total cross-sectional area).

**Figure 3.** Uniaxial tensile damage evolution process: (**a**) nominal stress–strain curve; (**b**) effective stress–strain curve; (**c**) microscopic spring beam model; (**d**) cylinder model.

### 3.1.4. Dialectical Unification between Degeneration and Evolution

Natural dialectics [49] holds that, the process of evolution in nature is the unity of opposites between evolution and degeneration, which exist and occur simultaneously. Evolution-oriented processes often inherently involve degradation, and vice versa. Evolution and degradation are two opposite trends in nature. They are closely combined and inseparable. Each side is the condition for the other side to occur. The combination of the two forms a circular spiral propulsion mode for the evolution of nature, which makes the evolution process of nature appear periodic.

A novel fundamental assumption of mesoscopic damage for concrete material was proposed by Chen et al. [35]. As shown in Figure 4a, it indicated that the damage evolution of concrete materials could be summarized as two kinds of mechanisms on a mesoscopic scale. On the one hand, with the increase of deformation, the initiation, propagation, and coalescence of micro-cracks and micro-defects, as well as the acoustic emission phenomenon, will occur randomly in the microstructure, which is the so-called degradation effect. On the other hand, due to the initiation and propagation of microcracks, stress redistribution and adjustment of force skeleton will take place in the microstructure. We may be able to understand the second mechanism as the active adjustment of the material system itself. In this way, the force skeleton of the microstructure is further adjusted and optimized, and the potential mechanical ability of the material is further developed to

bear more external loads (effective stress). Therefore, the second type of mechanism is called the strengthening effect here. Berthier et al. [24] indicated the quasi-brittle failure emerges from the interaction between the elements constituting the material. They also highlighted the central role played by the mechanism of load redistribution to control the failure behavior of quasi-brittle solids.

**Figure 4.** Two types of mesoscopic mechanism of quasi-brittle material: (**a**) damage evolution mechanism on a mesoscopic scale; (**b**) degradation factor; (**c**) evolutionary factor.

According to the assumption, the deformation and failure of concrete material is not only the process of "deterioration" of macroscopic and microscopic mechanical properties, but also the process of "evolution" in which the material goes through their own active-adjustment to adapt to the external load environment. When the inner adjustment capabilities of the material exert to its limit, the damage evolution of the material will switch from the statistical uniform damage to the local breach. This local failure process generally takes place in the weakest region (such as the fracture process zone in concrete).

The above viewpoints are consistent with the catastrophe theory [43], which holds that the failure process of solid materials could be divided into two phases, GS and EIC. The critical state transforming from GS to EIC plays a key role during the whole process, and exhibits the critical sensitivity. The two-stages embody the process from quantitative change to qualitative change.

Figure 4b,c shows the two mechanisms of the deformation and failure of concrete on a mesoscopic scale, i.e., the "degradation" and "strengthening", respectively.

As shown in Figure 4c, *E*<sup>v</sup> is introduced as the evolutionary factor to describe the "strengthening" process, corresponding to the yield damage mode. The expression is as follows:

$$E\_{\rm V} = \int\_{0}^{\varepsilon} p(\varepsilon) d\varepsilon \quad \text{(} 0 \le \varepsilon \le \varepsilon\_{\rm ymax} \text{)}\tag{11}$$

*E*<sup>v</sup> could be used to assess the extent to which the potential mechanical capabilities (adjustment capabilities of force skeleton in microstructure) of materials are developed, ranging from 0 to 1. When *E*<sup>v</sup> = 0, it corresponds to the initial undamaged state. When *E*<sup>v</sup> = 1, it corresponds to the critical state, at which the potential adjustment capabilities of materials reach their limits, *σ*<sup>E</sup> reaches its maximum value, and then the materials enter into the local catastrophic stage. The whole process embodies the characteristics of "quantum" to "qualitative", in which the yield damage mode plays a key role.

As shown in Figure 4b, *D*<sup>R</sup> is as the fracture damage variable, and could be used to describe the "degradation" processes. It has the similar physical meanings with *D* defined by Dougill [50], which characterizes the initiation and propagation of microcracks and leads to the reduction of the effective stress area. At the uniform damage stage, it ranges from 0 to H (a certain smaller value), and it seems not to be playing a key role in the whole process.
