An Anisotropic Damage Model of Quasi-Brittle Materials and Its Application to the Fracture Process Simulation

Abstract

Accurate predictions of the failure behaviors of quasi-brittle materials are of practical significance to underground engineering. In this work, a novel anisotropic damage model is proposed based on continuous damage mechanics. The anisotropic damage model includes a two-parameter parabolic-type failure criterion, a stiffness degradation model that considers anisotropic damage, and damage evolution equations for tension and shear, respectively. The advantage of this model is that the degradation of elastic stiffness only occurs in the direction parallel to the failure surface for shear damage, avoiding the interpenetration of crack surfaces. In addition, the shear damage evolution equation is established based on the equivalent shear strain on the failure face. A cyclic iterative method based on the proposed anisotropic damage model was developed to numerically simulate the fracture process of quasi-brittle materials. The developed model and method are important because the ready-made finite element software is difficult to simulate the anisotropic damage of quasi-brittle materials. The proposed anisotropic damage model was tested against a conventional damage model and validated against two benchmark experiments: uniaxial and biaxial compression tests and Brazilian splitting tests. The results demonstrate that the proposed anisotropic damage model simulates the mesoscale damage mode, macroscale fracture modes, and strength characteristics more effectively and accurately than conventional damage models.

1. Introduction

Quasi-brittle materials are widely used in civil and underground engineering. Among them, there are artificial materials, such as concrete and masonry, and natural materials, such as rock and hard clay. Quasi-brittle materials generally exhibit the following mechanical properties: inhomogeneity, anisotropy, structural discreteness, and nonlinearity. The failure of quasi-brittle materials causes engineering accidents and geological disasters, such as collapses, landslides, and tunnel collapses, creating challenges to the construction and operation of engineering structures. Therefore, exploring the failure behaviors of quasi-brittle materials is of great significance.

In the past decades, many research studies were carried out on material failure behaviors. Most of these studies are based on experimental methods, involving material splitting failure [1,2], impact failure [3], static uniaxial and triaxial compression failure [4,5,6], and dynamic uniaxial and triaxial compression failure [7,8]. This method is simple and easy to operate. However, it also has some disadvantages, such as complex sample preparation, long cycles, high cost, large dispersion of results, and unclear failure mechanism. With the development of computer technology and numerical computation, numerical simulations have been gradually used to study failures. A range of methods and models have been proposed to simulate the initiation and propagation of cracks, such as the traditional finite element method [9,10,11,12,13], discrete element method [14,15,16,17,18], extended finite element method [19,20,21], boundary element method [22,23,24], displacement discontinuity method [25,26,27,28,29], etc. The most popular ones are the finite element method and discrete element method. The discrete element method considers the material as a combination of particles and expresses cracks by bond-breaking between them. It is suitable to deal with discontinuities and large deformations in the joint rock mass. However, it encounters difficulties in simulating the deformation and failure of intact rocks. In addition, the discrete element method is relatively computationally intensive, which limits either the length of a simulation or the number of particles, and is difficult to apply to complex problems, such as multiscale and multifield problems.

The finite element method is mostly based on the damage model [30,31,32,33] or phase field model [34,35,36], and simulates the initiation and propagation of cracks through stiffness reduction. Compared with the discrete element method, this method has the advantages of a clear physical meaning of parameters, high calculation efficiency, and strong applicability. The phase field and damage modeling approaches have essential differences. The phase field method introduces an additional degree of freedom that greatly increases the dimension of the problem, and requires fine grids, thus requiring higher computing resources and being limited to small-scale problems. The damage method does not need to introduce an additional degree of freedom; thus, it is computationally efficient and suitable for large-scale engineering simulations. It assumes that materials may fail in tension or shear and introduce a damage variable to characterize the degradation of stiffness. The key to adopting the damage method to simulate material failure is to construct appropriate failure criteria, a stiffness degradation model, and the damage evolution equation.

The most-used tension failure criterion is the maximum tensile stress criterion. The commonly used shear failure criteria include the Mohr–Coulomb strength criterion [31,32,37,38] and the Drucker–Prager criterion [33,39]. These criteria approximate the envelope of the Mohr circles as a straight line. The error of this approximation is small under low-stress conditions but could be large under high-stress conditions. In addition, the simulation results based on these criteria in the tension–shear composite zone are far from the experimental results. Notably, it is impossible to describe the continuous change of failure mode with stress state because the connection between tension and shear failure envelopes is not smooth. In terms of the stiffness degradation model, most studies use the isotropic damage model. That is, stiffness degradation is assumed to be the same for all directions. This can be reasonable for tension damage, but not for compression shear damage. When compression shear damage occurs, stiffness weakens mainly along the crack plane direction. The isotropic damage assumption will cause the interpenetration of crack faces. With respect to the damage evolution equation, conventional damage models generally adopt the tension damage evolution equation based on tensile strain and the shear damage evolution equation based on compressive strain. As the tension cracks are opened and the two crack faces are separated, the physical meaning of the tensile strain as the primary indicator to measure the degree of tension damage is clear. The shear cracks are mostly closed and only a relative slip occurs between the two crack faces. Therefore, the adoption of compressive strain as a basic indicator for measuring the degree of shear damage is unreasonable.

Most studies use biaxial or triaxial compression tests of rock specimens to validate the damage models. A comparative analysis found that although these models simulate the general failure process, the following deficiencies still exist: (1) Tension–shear failure or compression–shear failure under low-stress conditions cannot be reflected. (2) The fracture morphology of the rock specimen exhibited a discontinuous mutation with confining pressure. (3) The fracture angle between the failure plane and the first principal stress plane did not change with the confining pressure. (4) There was no apparent regularity between the residual strength and confining pressure, which was inconsistent with the experimental results. (5) The shear crack faces interpenetrate each other and the relative sliding direction is unreasonable. The above problems existing in existing damage models cannot be ignored. The aim of this study is to propose a novel anisotropic damage model to more accurately simulate the meso- and macro-failure behaviors and strength characteristics of quasi-brittle materials.

Compared to conventional damage models, the anisotropic damage model has the following three improvements. First, a two-parameter parabolic-type yield criterion was used as the failure criterion. It can reflect the nonlinear strength characteristics of quasi-brittle materials. Moreover, it describes the pure tension, tension–shear, and compression–shear damages using only one continuous criterion, enabling the continuous transition of different failure modes. Second, an anisotropic stiffness degradation model was used instead of an isotropic degradation model to reflect the anisotropic mechanical characteristics of quasi-brittle materials after shear damage. Third, the shear strain on the failure plane was selected as the independent variable rather than compression strain to establish the shear damage evolution equation, for which the physical meaning is clearer. The proposed anisotropic damage model was implemented using COMSOL Multiphysics and MATLAB and validated against two benchmark simulations: uniaxial and biaxial compression tests and Brazilian splitting tests. The effectiveness and performance of the proposed anisotropic damage model were verified by comparing the simulation results with experimental results and other numerical simulation results.

The remainder of this paper is organized as follows: Section 2 presents a novel anisotropic damage model. Section 3 presents the implementation strategy of the numerical simulation. Section 4 evaluated the effectiveness of the proposed damage model. Section 5 presents the validation against the uniaxial and biaxial compression tests. Section 6 presents the validation against the Brazilian splitting test. Section 7 summarizes the applications of the proposed anisotropic damage model and discusses the scope for future research.

2. Anisotropic Damage Model

2.1. Failure Criterion

Two damage modes are typically assumed for quasi-brittle materials at the mesoscale: tension and shear. To unify the two types of damages with one criterion, the two-parameter parabolic-type yield criterion [40] was adopted as the failure criterion. It is assumed that the material will fail once the compressive stress σ normal to the failure plane and shear stress τ on the failure plane satisfy the following parabolic relationship:

$${\tau }^2=n(\sigma +{\sigma }_{\mathrm{t}0}),$$

(1)

where $n={\sigma }_{\mathrm{t}0}{(\sqrt{m+1}-1)}^2$, $m={\sigma }_{\mathrm{c}0}/{\sigma }_{\mathrm{t}0}$, σ_t0 is the uniaxial tensile strength, and σ_c0 is the uniaxial compressive strength. In the absence of compressive stress, this failure criterion is reduced to $\tau =\sqrt{n{\sigma }_{\mathrm{t}0}}$, where $\sqrt{n{\sigma }_{\mathrm{t}0}}$ is the shear strength of the failure plane without normal stress.

Figure 1 shows the two-parameter parabolic-type yield envelope and the Mohr circles. It can be observed that the parabolic envelope shows strong nonlinear characteristics in both tensile and compressive stress regions, making up for the deficiency of the linear yield criterion. The left-end point of the envelope corresponds to the uniaxial tensile strength of the material. According to the relationship between the Mohr circles and strength envelope in the critical state, the shear and compressive stresses on the failure plane are as follows:

$$\tau =\frac{{\sigma }_1-{\sigma }_3}{2}\cos \phi ,$$

(2)

$$\sigma =\frac{{\sigma }_1+{\sigma }_3}{2}-\tau \tan \phi ,$$

(3)

where φ is the internal friction angle; σ₁ and σ₃ are the first and third principal stresses, respectively. Moreover, the compressive stresses are positive.

According to the failure condition Equation (1), tan φ = ∂τ/∂σ = n/2τ. Therefore, the stresses on the failure plane are as follows:

$$\sigma =\frac{{\sigma }_1+{\sigma }_3}{2}-\frac{n}{2},$$

(4)

$$\tau =\frac{{\sigma }_1-{\sigma }_3}{2}\cdot \frac{\sqrt{2({\sigma }_1+{\sigma }_3)-2n+4{\sigma }_t}}{\sqrt{2({\sigma }_1+{\sigma }_3)-n+4{\sigma }_t}}.$$

(5)

Hence, the failure condition can be expressed in terms of the principal stresses as follows:

$$\{\begin{array}{r}\begin{array}{cc}-{\sigma }_3\ge {\sigma }_{t0},& {\sigma }_1\le n-{\sigma }_{t0}\end{array}\\ \begin{array}{cc}{({\sigma }_1-{\sigma }_3)}^2-2n({\sigma }_1+{\sigma }_3)-4n{\sigma }_{t0}+n^2\ge 0,& {\sigma }_1>n-{\sigma }_{t0}\end{array}\end{array}.$$

(6)

When σ₁ ≤ n − σ_t0, the material fails in tension. When σ₁ > n − σ_t0, it fails in shear. Therefore, σ₁ = n − σ_t0 is the critical principal stress for switching the failure mode. In addition, for shear failure, $\frac{{\sigma }_1+{\sigma }_3}{2}-\frac{n}{2}\ge 0$ corresponds to compression–shear failure, whereas $\frac{{\sigma }_1+{\sigma }_3}{2}-\frac{n}{2}<0$ corresponds to the tension–shear failure.

It is necessary to further establish the failure conditions expressed by strains for shear failure. According to the shear strength of the failure plane without normal stress and Hooke’s law for shear, the failure initiation shear strain on the failure plane is as follows:

$${\gamma }_0=\frac{2(1+\mu )}{E}\sqrt{n{\sigma }_{\mathrm{t}}}.$$

(7)

Hence, the shear failure condition can be expressed as follows:

$${\gamma }_{\mathrm{eq}}={\gamma }_0,$$

(8)

where γ_eq is the equivalent shear strain on the failure plane, deduced using Equation (6), as follows:

$${\gamma }_{\mathrm{eq}}=\sqrt{{({\epsilon }_1-{\epsilon }_3)}^2+\alpha ({\epsilon }_1+{\epsilon }_3)+\beta },$$

(9)

where $\alpha =-\frac{2n{(1+\mu )}^2}{(1-\mu )E}$ and $\beta =\frac{n^2{(1+\mu )}^2}{E^2}$.

As shown in Equation (8), the equivalent shear strain based on the two-parameter parabolic-type yield criterion is a nonlinear function of the strain, differing from the existing equivalent strain based on the linear yield criterion [41].

2.2. Stiffness Degradation Model

Based on the theory of continuous damage mechanics [42], the stress–strain relationship of the damaged zone can be established by introducing the degenerated elastic stiffness coefficient C_ijkl as follows:

$${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}.$$

(10)

Under the assumption of isotropic damage, the degenerated stiffness coefficient C_ijkl is expressed as

$$C_{ijkl}=(1-D)C_{ijkl}^0,$$

(11)

where $C_{ijkl}^0=\lambda {\delta }_{ij}{\delta }_{kl}+\mu \left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)$ is the initial elastic stiffness coefficient for isotropic materials, $\lambda =\frac{\nu E}{\left(1+\nu \right)\left(1-2\nu \right)}$ and $\mu =\frac{E}{2\left(1+\nu \right)}$ are Lamé’s constants, and E and ν are Young’s modulus and Poisson’s ratio, respectively. D is the damage variable reflecting the damage state of the material. D = 0 corresponds to the undamaged state, D = 1 corresponds to the complete damage (fracture or failure) state, and 0 < D < 1 corresponds to the different degrees of damage.

The above model applies to tension damage. For compression shear damage, the damage parallel to the failure surface is usually larger, and the damage perpendicular to the failure surface is usually smaller and can be ignored. This will lead to obvious weakening of the shear stress parallel to the failure surface, while the normal stress vertical to the failure surface is basically unaffected. Therefore, it is not necessary to weaken the stiffness in all directions, but only in the direction parallel to the failure surface.

Assuming the unit normal tensor and tangent tensor of the failure surface are n and t respectively, then the shear stress on the failure surface is

$${\tau }_{\alpha }=\frac{1}{2}{\beta }_{ij}{\sigma }_{ij}.$$

(12)

where ${\beta }_{ij}=n_i{t}_j+t_i{n}_j$

Substituting Equation (10) into Equation (12) and letting D = 0, yields

$${\tau }_{\alpha }=\mu {\beta }_{kl}{\epsilon }_{kl}.$$

(13)

Hence, the contribution of shear stress on the failure surface to the stress tensor is

$${\sigma }_{ij}^{\tau }={\tau }_{\alpha }{\beta }_{ij}=C_{ijkl}^{\tau }{\epsilon }_{kl}.$$

(14)

where $C_{ijkl}^{\tau }=\mu {\beta }_{ij}{\beta }_{kl}$ is the elastic stiffness coefficient contributing to the shear stress of the failure surface. Assuming the angle between the normal vector of the failure plane and the positive X-axis is α, $\mathrm{then}{C}_{ijkl}^{\tau }$ is as follows:

$$C_{ijkl}^{\tau }=\left[\begin{array}{cccccc}\mu {(\sin 2\alpha )}^2& -\mu {(\sin 2\alpha )}^2& 0& -\mu \sin 2\alpha \cos 2\alpha & 0& 0\\ -\mu {(\sin 2\alpha )}^2& \mu {(\sin 2\alpha )}^2& 0& \mu \sin 2\alpha \cos 2\alpha & 0& 0\\ 0& 0& 0& 0& 0& 0\\ -\mu \sin 2\alpha \cos 2\alpha & \mu \sin 2\alpha \cos 2\alpha & 0& \mu {(\cos 2\alpha )}^2& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right],$$

(15)

Hence, the stiffness degradation model considering anisotropic damage can be written as follows:

$$C_{ijkl}=C_{ijkl}^0-D{C}_{ijkl}^{\tau }.$$

(16)

2.3. Damage Evolution Equations

In this section, tension and shear damage evolution equation were established based on simplified constitutive relationships respectively. It should be noted that the shear damage evolution equation was proposed based on the equivalent shear strain on the failure face.

2.3.1. Tension Damage

An elastic–brittle strain-softening model was used for tension damage [43], as shown in the stress–strain relationship in Figure 2.

Hence, the damage variable can be obtained as follows:

$$D=\{\begin{array}{l}0\begin{array}{cc}& \overline{\epsilon }<{\epsilon }_{t0}\end{array}\\ 1-\frac{{\lambda }_{\mathrm{t}}\epsilon }{\overline{\epsilon }}\begin{array}{cc}& {\epsilon }_{t0}\le \overline{\epsilon }<{\epsilon }_{tu}\end{array}\\ 1\begin{array}{cc}& \overline{\epsilon }\ge {\epsilon }_{tu}\end{array}\end{array},$$

(17)

where λ_t = σ_tr/σ_t0 is the residual tensile strength coefficient; σ_tr is the residual tensile strength; ε_t0 = σ_t0/E₀ is the elastic ultimate tensile strain; ε_tu is the ultimate tensile strain before complete failure; $\overline{\epsilon }=\sqrt{\langle {\epsilon }_1\rangle {}^2+\langle {\epsilon }_2\rangle {}^2+\langle {\epsilon }_3\rangle {}^2}$ is the equivalent tensile strain and ε₁, ε₂, ε₃ are the first, second, third principal strains, respectively.

2.3.2. Shear Damage

An elastic–brittle–plastic strain-softening model was used for shear damage, as shown in the shear stress–shear strain relationship in Figure 3.

Hence, the damage variable can be obtained as follows:

$$D=\{\begin{array}{l}0\begin{array}{cc}& {\gamma }_{\mathrm{eq}}<{\gamma }_0\end{array}\\ 1-\frac{{\lambda }_{\mathrm{s}}{\gamma }_0}{{\gamma }_{\mathrm{eq}}}\begin{array}{cc}& {\gamma }_{\mathrm{eq}}\ge {\gamma }_0\end{array}\end{array},$$

(18)

where λ_s = τ_r /τ₀ is the residual shear strength coefficient, τ_r is the shear residual strength, and τ₀ is the shear strength.

Owing to the nonlinear relationship between the equivalent shear strain γ_eq and principal strain component, the nonlinearity of the damage evolution equations is further increased compared with the conventional damage mode. Moreover, the complexity and uncertainty of the damage problem are also significantly increased.

3. Simulation of the Fracture Process

During the fracture process of quasi-brittle materials, some elements with a large elastic modulus and low strength fail first as the load increases. This reduces the bearing capacity, redistributing the stress and resulting in the successive damage and failure of other elements. In the mathematical model, the element damage is measured using the damage variable D, which is affected by the element’s stress and strain. The element’s stress and strain are also affected by its elastic modulus, which is related to the damage variable D. Therefore, the damage problem is a strongly nonlinear problem that cannot be solved using only ready-made finite element software. This study used a cyclic iterative method based on a joint simulation of MATLAB and COMSOL to solve this problem. Figure 4 shows a detailed calculation flowchart in which the COMSOL Multiphysics or MATLAB icon in each step indicates the software used in the simulation. First, set the initial damage variable D₀ to 0 before loading. Then, in each loading step determine the elastic tensor according to the initial damage variable D₀, and calculate the stress, strain, and new damage variable D₁. Judge whether the difference between D₁ and D₀ meets the specified tolerance. If yes, go to the next loading step. If not, assign D₁ to D₀, update the elastic tensor, and recalculate the stress, strain, and new damage variable D₁. Repeat the above process until the last loading step.

4. Test of the Anisotropic Damage Model

This section evaluates the effectiveness of the proposed anisotropic damage model. A rectangular body, 25 mm × 50 mm, subjected to an increasing downward displacement with confining pressure pc, was analyzed, as shown in Figure 5. The rectangle was divided into 12,806 triangular grids with a maximum length of 0.5 mm. The vertical displacement of the bottom boundary and the horizontal displacement of the lower left corner are fixed. The material parameters were considered as follows: E = 60 GPa, ν = 0.3, σ_c0 = 200 MPa, and λ_s = 0.1. The damage variable D was zero everywhere except within the crack band, where D was computed using Equation (12). Figure 6 shows the distribution of the damage variable and elastic modulus after failure. The proposed anisotropic damage model with two cases of confining pressure, p_c = 0 MPa and p_c = 10 MPa, was compared with the conventional isotropic damage model in the literature [44,45]. In the conventional isotropic damage model, the failure criterion is the Mohr–Coulomb strength criterion, and the shear damage evolution equation is based on compressive strain.

Figure 6 shows the simulation results of failure modes corresponding to the two damage models. For the anisotropic damage model, only tangential slip occurs on the failure surface. For the isotropic damage model, the tangential slip and vertical compression occur on the failure surface. It is obviously unreasonable that the two fracture surfaces interpenetrate each other.

5. Uniaxial and Biaxial Compression Tests

Influenced by the test cases with the conventional damage model in reference [44,45], the failure processes of a rectangular rock specimen under uniaxial and biaxial compression were simulated for direct comparison. As shown in Figure 7a, the rectangular specimen had the same size and boundary conditions as the case in Section 4, except that the material was assumed to be mesoscale heterogeneous. Figure 7b shows the mesh grids and the initial elastic modulus of the rectangular specimen. The rectangle was divided into 25 × 50 = 1250 square grids with side length 1 mm.

The mechanical parameters of each grid element were considered as random variables with a Weibull distribution, whose probability density function can be described as follows:

$$f(u)=\frac{m}{u_0}{\left(\frac{u}{u_0}\right)}^{m-1}\exp \left[-{\left(\frac{u}{u_0}\right)}^m\right],$$

(19)

where u is a random variable; that is, the mechanical parameters of an individual grid element, such as elastic modulus, Poisson’s ratio, and uniaxial tensile strength. u₀ is the scale parameter related to the average value of variable u. m is the shape parameter, the homogeneity index, which reflects the degree of material homogeneity.

In this test case, the following mechanical parameters of the rock specimen were adopted: m = 4, E₀ = 60 GPa, σ_t0 = 20 MPa, n₀ = 107, and λ_t = λ_s = 0.1. Moreover, four confining pressures were studied: pc = 0 MPa, 10 MPa, 20 MPa, and 40 MPa. Confining pressure was imposed on the specimen in the first loading step, and displacement increment was applied to the top boundary in the subsequent loading steps. A numerical sensitivity analysis of the loading step length on the uniaxial compression test indicated that the loading step length in a reasonable quasi-static range did not affect the macroscopic strength and the main crack propagation significantly [31]. Considering calculation accuracy and efficiency, a constant displacement increment of 0.005 mm/step was used here.

5.1. Damage Evolution

Figure 8, Figure 9, Figure 10 and Figure 11 show the failure processes of rectangular specimens under different confining pressures. The figures show the five steps of damage: primary damage, secondary damage, pre-peak damage, post-peak damage, and final damage.

The mesoscale failure process can be analyzed through evolutionary diagrams of the damage variable, in which red represents the damaged elements. It can be observed that the loading step of the primary damage varied with different confining pressures. For low confining pressures, such as 0 MPa and 10 MPa, the primary damage occurred after a certain number of loading steps when the stresses of some elements exceeded their strength owing to an increase in the axial load. For confining pressures greater than or equal to 20 MPa, primary damage occurred in the first step when the axial load was not applied, and the damage was only caused by the confining pressure. The greater the confining pressure, the greater the number of primary damaged elements. The confining pressure also affected the loading step of the secondary damage. The greater the confining pressure, the later the secondary damage occurred. Owing to the randomness of the material properties of the specimens, the early damaged elements were arbitrarily scattered. In general, the damage developed slowly and accumulated insignificantly before the peak load. The damage spread rapidly and accumulated significantly once the peak load was attained, finally forming a macrocrack penetrating the rock specimen. Subsequently, the damage stopped because the specimen deformation mainly occurred in the completely damaged area. Moreover, the damage became more concentrated as the confining pressure increased.

The failure mechanism can be analyzed from the evolutionary diagrams of the damage mode. For low confining pressures, early damage before the peak load was dominated by compression–shear damage with minimal pure tensile damage and sporadic tensile–shear damage. Late damage after the peak load was dominated by pure tensile damage. With an increase in the confining pressures, the proportion of compression–shear damage in early damage increased gradually. Moreover, late damage also gradually changed to compression–shear-dominated damage.

5.2. Fracture Modes

The fracture modes of rock specimens under different confining pressures are shown in Figure 12, where the simulation results show the final distribution of the elastic modulus. It can be seen that the fracture surfaces were rough and complex under low confining pressures. In contrast, the fracture surface became smoother and more regular with increased confining pressure. The fracture angles between the fracture plane and longitudinal axis were small for low confining pressures. In contrast, the fracture angle increased with an increase in confining pressure. The fracture angle was predicted to stabilize at approximately 45°. The fracture morphology and orientation of the rock specimen exhibited the same strong graduality and continuity with confining pressure as the experimental results [46,47,48], which cannot be achieved by the conventional damage model [45].

5.3. Stress–Strain Curves

Figure 13 shows the stress–strain curves under different confining pressures. These stress–strain curves can be roughly divided into four stages: elastic deformation, stable damage, rapid damage, and complete fracture. In the elastic deformation stage, axial loading did not cause damage, and the stress–strain curve was linear. In the stable damage stage, known as the plastic deformation stage, the slope of the stress–strain curve decreased gradually because the local stress exceeded the material strength, and more damage was initiated in the specimen. The stress increased slowly at this stage until the peak stress was attained. The stress at the starting point of this stage was approximately two-thirds of the peak stress. In the rapid damage stage, the stress–strain curve dropped sharply from the peak and then remained stable in the complete fracture stage. It can be seen that for confining pressures smaller than 20 MPa, the post-peak behavior in compression showed a typical brittle behavior; that is, the axial stress dropped immediately once it reached its peak. With increased confining pressure, plastic deformation became apparent, strain softening occurred, and the failure gradually changed from brittle to ductile.

In general, the tendency of the simulated stress–strain curves matched well with the experimental curves [47], the data of which are scaled for comparison, as shown in Figure 13. The main deviation difference is that the elastic modulus obtained by numerical simulation does not increase with confining pressure as in the experiment. The main reason for this deviation is that the microcracking closure and specimen compaction were not considered in the numerical simulation.

Figure 14 shows the peak and residual strengths under different confining pressures. The peak and residual strengths increased gradually with the increase in the confining pressure, showing a better regularity than the numerical simulation results with conventional damage models [44].

The numerical simulation results were also well fitted by the generalized Hoek–Brown failure criterion [49], expressed as follows:

$${\sigma }_1={\sigma }_3+{\sigma }_{\mathrm{c}}\sqrt{m{\sigma }_3/{\sigma }_{\mathrm{c}}+s},$$

(20)

where σ_c is 74, the fitting parameter m is 0.9, and s is 1 or 0.

6. Brazilian Splitting Test

A numerical simulation of the Brazilian splitting test was also performed to evaluate the effectiveness of the developed damage model and simulation procedure. A flattened Brazilian disc specimen placed between two loading plates with a diameter of 50 mm and a loading angle of 20° was considered in the simulation, as shown in Figure 15a. Boundary conditions were imposed on the top and bottom loading plates. The mechanical parameters of the disc specimen were the same as those of the rectangular specimen in Section 5. The elastic modulus and Poisson’s ratio of the loading plates were 5000 GPa and 0.3, respectively. The mesh grids and the initial elastic modulus of the disc specimen are shown in Figure 15b. A constant displacement increment of 0.005 mm/step was used in this simulation.

Figure 16 and Figure 17 show the failure processes and load–displacement curves of the disc specimen, respectively. The early damage before the peak load was mainly pure tensile and distributed near the center of the specimen, which is consistent with the simulation result of another method [50]. The load–displacement curve at this moment remained linear because only a few elements were damaged. As the microcracking closure and specimen compaction were not considered, the simulated load–displacement curve cannot show the concave shape as in the experiment.

Once the peak load was reached, the damaged elements rapidly propagated, connected, and penetrated the specimen, resulting in a splitting crack along the loading direction. Moreover, the external load decreased sharply to a certain value. The entire failure process exhibited evident brittle characteristics consistent with the experiment. Pure tensile damage was also dominant in the final fracture, while shear damage, especially tension–shear damage, was rare and negligible. Figure 18 shows the final fracture mode along with the experimental results. The simulation results agree with the experimental results [51,52].

7. Conclusions

• A novel anisotropic damage model is proposed to improve the conventional isotropic damage model. It avoids the interpenetration of crack surfaces and unifies three forms of damage, which include pure tension, tension–shear, and compression–shear, as one smooth criterion. In addition, the shear damage evolution equation considers the equivalent shear strain on the failure face as the independent variable, which has a clearer physical meaning.
• The proposed anisotropic damage model is able to simulate the failure behavior of quasi-brittle materials more effectively. It can effectively simulate not only the mesoscale pure tension, tension shear, and compression shear damages but also the macroscale fracture mode and its evolution. It can also simulate the strength characteristics, including the peak and residual strengths, approximately.
• The simplified two-dimensional geometries were used to verify the proposed anisotropic damage model by comparison with a conventional damage model. However, the proposed anisotropic damage model and finite element implementation strategy proposed in this study also apply to three-dimensional, multi-scale, and multi-field coupling problems.
• Numerical tests of the failure behavior of quasi-brittle materials have the advantages of strong universality, convenience, flexibility, and repeatability compared with laboratory tests. In addition, they can provide mesoscale mechanical information that is difficult to observe in laboratory tests. However, owing to the limitations of the nonlinear solution method, numerical tests will take longer if there are more loading steps. The combination of numerical simulation and machine learning algorithms will help improve the efficiency of numerical tests better to serve the scientific research and teaching of quasi-brittle materials, which is also a follow-up research goal of the authors.