Numerical Modeling of Quasi-Brittle Materials Using a Phase-Field Regularized Cohesive Zone Model with Optimal Softening Law

: In this paper, we propose an approach combining optimal softening laws and a phase-ﬁeld regularized cohesive zone model (PF-CZM) for modeling the fracture and damage properties of quasi-brittle materials accurately. In this method, the optimal softening law is determined by comparing the predicted results with experimental data in the framework of the PF-CZM; three typical softening laws are considered. The PF-CZM with a length scale is used to model crack initiation and propagation without considering the mesh bias. We ﬁrst investigate the mechanical responses and crack propagations of different concrete beams based on the above approach; the predicted results are compared with the data from conventional methods and experiments. The results indicate that the mechanical properties of concrete beams with the optimal softening law are better than the data reported in the literature. Further validation indicates that once the optimal softening law is determined, it is stable for the same group of materials. Moreover, we demonstrate that the PF-CZM can naturally predict and reproduce the critical notch offset and fracture transition process of three-point bending concrete beams and the fracture features of typical double-notched concrete beams, such as the interaction between two notches objectively, together with the changes of limit load capacity.


Introduction
Quasi-brittle materials such as concrete and fiber-reinforced concrete are most commonly used in engineering structures.Predicting the limit load capacity and failure behavior is of vital importance in preventing the potential risks of these structures due to the extremely stochastic and nonlinear properties of quasi-brittle materials [1,2].The accurate analysis of failure in mechanical response, crack initiation and propagation via computational modeling is not only an indispensable tool for which systematic experiments are expensive, time-cost, difficult or even impractical, but it also provides insights into understanding the failure processes and mechanical behaviors of these materials, such as in the third numerical example provided in this work [3].
For quasi-brittle materials with strain softening and intrinsically heterogeneous properties, their fracture process and mechanical behavior cannot be analyzed simply by elastic or plastic mechanics methods.In the past decades, various approaches have been proposed to understand and model the mechanical response, crack initiation and propagation of quasi-brittle materials under different conditions [4][5][6].These approaches can be divided into two groups: damage mechanics methods and fracture mechanics methods.However, conventional damage mechanics methods often suffer from mesh bias problems due to the strain localization phenomenon [4,5] unless some additional methods are proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems exist in the smeared/diffusive and discrete methods of fracture mechanics.
In fracture mechanics, there are many methods currently proposed for crack propagation modeling.These methods may be classified into discrete and smeared/diffusive methods.Discrete methods model cracks as real discontinuities.Representative methods are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obtain mesh-independent results [9][10][11].The XFEM requires an explicit representation of the crack and additional criteria to achieve complex crack growth, and its convergence is not robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can be seen as a development of continuum damage mechanics, which treat the cracks as the result of damage accumulation.As the most popular diffuse approach, the phase field (PF) method has recently received more and more attention from researchers for its advantage of not requiring additional criteria to predict crack nucleation and propagation [1,2,14].The method started with the linear-elastic fracture variational principle proposed by Francfort and Marigo [15] Ψ(u, Γ) = Ψ b + Ψ Γ = Ω However, conventional damage mechanics methods often suffer from mesh bias problems due to the strain localization phenomenon [4,5] unless some additional methods are proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems exist in the smeared/diffusive and discrete methods of fracture mechanics.
In fracture mechanics, there are many methods currently proposed for crack propagation modeling.These methods may be classified into discrete and smeared/diffusive methods.Discrete methods model cracks as real discontinuities.Representative methods are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obtain mesh-independent results [9][10][11].The XFEM requires an explicit representation of the crack and additional criteria to achieve complex crack growth, and its convergence is not robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can be seen as a development of continuum damage mechanics, which treat the cracks as the result of damage accumulation.As the most popular diffuse approach, the phase field (PF) method has recently received more and more attention from researchers for its advantage of not requiring additional criteria to predict crack nucleation and propagation [1,2,14].The method started with the linear-elastic fracture variational principle proposed by Francfort and Marigo [15] (, ) =  +  = (, )d +  d (1) where the total potential energy functional (, ) is assumed to include the surface energy  of the crack and the bulk energy  ; more details will be discussed in Section 2. Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made further extensions to the PF model, making the phase field modeling for failure generalized.PF methods have been demonstrated to be a powerful modeling tool due to their elegant format for the prediction of failure in many engineering materials such as rock [19], asphalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] and shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality between the classical PF model and the nonlocal damage model and proposed a unified phase-field theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM) derived from the above theory has been validated to be able to overcome most of the drawbacks mentioned in the above models, such as mesh bias and length scale sensitivity, and even the size effect can be captured accurately [26,27].See [28,29] for recent progress.
Multiscale methods are also adopted by combining the above approaches for the accurate analysis of quasi-brittle materials [30,31].Bharali et al. [32] developed a two-scale phase field fracture framework and discussed computational homogenization.However, the representative volume element (RVE) concept is not completely in accordance with actual quasi-brittle materials [2,4].A more popular strategy is considering the local multiscale, i.e., only the potential cracking zone is considered with multiphase at the mesoscale [31].In this way, both indirect and direct approaches are often adopted.In the indirect approaches, the critical fracture parameters of the materials, such as tensile strength and fracture energy, are modeled as stochastic variables with random fields so that different mixtures at the mesoscale are modeled implicitly.Recently, Li et al. [33] developed a method combining the PF-CZM and random fields for quasi-brittle materials in a 2D setting; only tensile strength is modeled by random fields.Zhang et al. [34] analyzed the size effect at the mesoscale with the same method.Huang et al. [35] adopted statistical reconstruction algorithms and the Karhunen-Loève expansion technique with microscale X-ray computed tomography (CT) images in order to get the random fields of tensile strength.In addition, Lu et al. [36] generated a random field with a stochastic harmonic function (SHF).On the other hand, various phases at the mesoscale are generally taken into account explicitly in the direct approaches.Xia et al. [37] investigated the mechanical behaviors of concrete in a 2D mesoscale setting, with the aggregates generated using a random sequential addition (RSA) method.Li et al. [38] developed a 3D mesoscale method for the damage and fracture simulation of quasi-brittle materials with the PF- where the total potential energy functional Ψ(u, Γ) is assumed to include the surface energy Ψ Γ of the crack and the bulk energy Ψ b ; more details will be discussed in Section 2. Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made further extensions to the PF model, making the phase field modeling for failure generalized.PF methods have been demonstrated to be a powerful modeling tool due to their elegant format for the prediction of failure in many engineering materials such as rock [19], asphalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] and shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality between the classical PF model and the nonlocal damage model and proposed a unified phase-field theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM) derived from the above theory has been validated to be able to overcome most of the drawbacks mentioned in the above models, such as mesh bias and length scale sensitivity, and even the size effect can be captured accurately [26,27].See [28,29] for recent progress.
Multiscale methods are also adopted by combining the above approaches for the accurate analysis of quasi-brittle materials [30,31].Bharali et al. [32] developed a two-scale phase field fracture framework and discussed computational homogenization.However, the representative volume element (RVE) concept is not completely in accordance with actual quasi-brittle materials [2,4].A more popular strategy is considering the local multiscale, i.e., only the potential cracking zone is considered with multiphase at the mesoscale [31].In this way, both indirect and direct approaches are often adopted.In the indirect approaches, the critical fracture parameters of the materials, such as tensile strength and fracture energy, are modeled as stochastic variables with random fields so that different mixtures at the mesoscale are modeled implicitly.Recently, Li et al. [33] developed a method combining the PF-CZM and random fields for quasi-brittle materials in a 2D setting; only tensile strength is modeled by random fields.Zhang et al. [34] analyzed the size effect at the mesoscale with the same method.Huang et al. [35] adopted statistical reconstruction algorithms and the Karhunen-Loève expansion technique with microscale X-ray computed tomography (CT) images in order to get the random fields of tensile strength.In addition, Lu et al. [36] generated a random field with a stochastic harmonic function (SHF).On the other hand, various phases at the mesoscale are generally taken into account explicitly in the direct approaches.Xia et al. [37] investigated the mechanical behaviors of concrete in a 2D mesoscale setting, with the aggregates generated using a random sequential addition (RSA) method.Li et al. [38] developed a 3D mesoscale method for the damage and fracture simulation of quasi-brittle materials with the PF-CZM; the aggregates were built using a random generating and packing algorithm.Huang et al. [39] proposed a new mesoscale modeling method in which the PF-CZM is combined with the cell-based smoothed finite element method (CSFEM) to improve computational efficiency.However, it is not easy to extend these methods to engineering applications directly.First of all, most studies focus on small specimens, and the distribution of both aggregates and random fields is assumed to be ideal.Moreover, the extensive Monte Carlo simulations (MCS) cost a lot, even for the RVEs.Additionally, the determination of the micromechanical parameters usually needs extra discussion, and all these factors may lead to an unreal deviation.
Another promising strategy is offered by peridynamics (PD), which substitutes integrodifferential motion equations for partial differential equations (PDEs).These equations are always applicable whether the displacement field is continuous or discontinuous.PD has been drawing a lot of attention for its simple and powerful ability to capture complicated failure processes.Li et al. [40] proposed a mesoscale PD for modeling the cracking process in concrete.To properly predict cohesive crack propagation, Yang et al. [41] proposed a PD-CZM and an approach combining PD and the finite element method (FEM).See [42] for a comparative review between PD and PF models.In this work, the comparison between the PF-CZM and the modified PDs will be provided as the second numerical example.
It can be seen that most of the existing studies focus mainly on the realization of the failure simulation of quasi-brittle materials.However, there have been a few studies on the effects of softening properties in numerical simulations.One key point is the inherent drawbacks of the various models.Inspired by the work by Wang [43] and Wu [44] and the excellent properties (e.g., no mesh bias, accuracy and insensitivity to the length scale parameter, etc.) of the PF-CZM, it is thus possible to focus on the softening laws in the accurate analysis of quasi-brittle materials.Specifically, this work is aimed at proposing a simple new approach to modeling the fracture and damage properties of quasi-brittle materials more accurately and efficiently.First, the optimal softening law is determined by comparing the numerical results using several softening laws in the framework of the PF-CZM with the experimental results.After that, several examples are investigated with the above approach to validate and extend the application scope of the PF-CZM in detail.The model is conducted using Abaqus software combined with the UMAT subroutine in order to take advantage of its built-in nonlinear solver, employing the Newton-Raphson algorithm and some other modules.The rest of this paper is organized as follows: In Section 2, we introduce the unified phase field theory and the PF-CZM.The approach and a numerical example of the optimal softening law are described in Section 3.Then, in Section 4, we investigate qualitatively and quantitatively the fracture processes and the mechanical responses of concrete beams with different features under three-point bending based on the PF-CZM with optimal softening law.Finally, a summary of this work and some suggestions are given in Section 5.

Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area A arranged along the x-axis, assume that the rod is fully cracked at x = 0, as shown in Figure 1.An auxiliary field variable φ(x) ∈ [0, 1] can be introduced to describe the sharp crack topology; then, we can get Equation (2) based on the variational principle.The crack phase field φ(x) located at x = 0 can be derived by (3); see more details in [1,18].
CZM; the aggregates were built using a random generating and packing algorithm.Huang et al. [39] proposed a new mesoscale modeling method in which the PF-CZM is combined with the cell-based smoothed finite element method (CSFEM) to improve computational efficiency.However, it is not easy to extend these methods to engineering applications directly.First of all, most studies focus on small specimens, and the distribution of both aggregates and random fields is assumed to be ideal.Moreover, the extensive Monte Carlo simulations (MCS) cost a lot, even for the RVEs.Additionally, the determination of the micromechanical parameters usually needs extra discussion, and all these factors may lead to an unreal deviation.
Another promising strategy is offered by peridynamics (PD), which substitutes integrodifferential motion equations for partial differential equations (PDEs).These equations are always applicable whether the displacement field is continuous or discontinuous.PD has been drawing a lot of attention for its simple and powerful ability to capture complicated failure processes.Li et al. [40] proposed a mesoscale PD for modeling the cracking process in concrete.To properly predict cohesive crack propagation, Yang et al. [41] proposed a PD-CZM and an approach combining PD and the finite element method (FEM).See [42] for a comparative review between PD and PF models.In this work, the comparison between the PF-CZM and the modified PDs will be provided as the second numerical example.
It can be seen that most of the existing studies focus mainly on the realization of the failure simulation of quasi-brittle materials.However, there have been a few studies on the effects of softening properties in numerical simulations.One key point is the inherent drawbacks of the various models.Inspired by the work by Wang [43] and Wu [44] and the excellent properties (e.g., no mesh bias, accuracy and insensitivity to the length scale parameter, etc.) of the PF-CZM, it is thus possible to focus on the softening laws in the accurate analysis of quasi-brittle materials.Specifically, this work is aimed at proposing a simple new approach to modeling the fracture and damage properties of quasi-brittle materials more accurately and efficiently.First, the optimal softening law is determined by comparing the numerical results using several softening laws in the framework of the PF-CZM with the experimental results.After that, several examples are investigated with the above approach to validate and extend the application scope of the PF-CZM in detail.The model is conducted using Abaqus software combined with the UMAT subroutine in order to take advantage of its built-in nonlinear solver, employing the Newton-Raphson algorithm and some other modules.The rest of this paper is organized as follows: In Section 2, we introduce the unified phase field theory and the PF-CZM.The approach and a numerical example of the optimal softening law are described in Section 3.Then, in Section 4, we investigate qualitatively and quantitatively the fracture processes and the mechanical responses of concrete beams with different features under three-point bending based on the PF-CZM with optimal softening law.Finally, a summary of this work and some suggestions are given in Section 5.

Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area  arranged along the -axis, assume that the rod is fully cracked at  = 0, as shown in Figure 1.An auxiliary field variable () ∈ [0, 1] can be introduced to describe the sharp crack topology; then, we can get Equation (2) based on the variational principle.The crack phase field () located at  = 0 can be derived by (3); see more details in [1,18].
CZM; the aggregates were built using a random generating and packing algo Huang et al. [39] proposed a new mesoscale modeling method in which the PF-C combined with the cell-based smoothed finite element method (CSFEM) to improv putational efficiency.However, it is not easy to extend these methods to engineeri plications directly.First of all, most studies focus on small specimens, and the distri of both aggregates and random fields is assumed to be ideal.Moreover, the ext Monte Carlo simulations (MCS) cost a lot, even for the RVEs.Additionally, the de nation of the micromechanical parameters usually needs extra discussion, and al factors may lead to an unreal deviation.
Another promising strategy is offered by peridynamics (PD), which substitu tegrodifferential motion equations for partial differential equations (PDEs).These tions are always applicable whether the displacement field is continuous or disco ous.PD has been drawing a lot of attention for its simple and powerful ability to c complicated failure processes.Li et al. [40] proposed a mesoscale PD for modeli cracking process in concrete.To properly predict cohesive crack propagation, Yan [41] proposed a PD-CZM and an approach combining PD and the finite element m (FEM).See [42] for a comparative review between PD and PF models.In this wo comparison between the PF-CZM and the modified PDs will be provided as the s numerical example.
It can be seen that most of the existing studies focus mainly on the realization failure simulation of quasi-brittle materials.However, there have been a few stud the effects of softening properties in numerical simulations.One key point is the in drawbacks of the various models.Inspired by the work by Wang [43] and Wu [4 the excellent properties (e.g., no mesh bias, accuracy and insensitivity to the lengt parameter, etc.) of the PF-CZM, it is thus possible to focus on the softening laws accurate analysis of quasi-brittle materials.Specifically, this work is aimed at propo simple new approach to modeling the fracture and damage properties of quasi-brit terials more accurately and efficiently.First, the optimal softening law is determin comparing the numerical results using several softening laws in the framework of t CZM with the experimental results.After that, several examples are investigated w above approach to validate and extend the application scope of the PF-CZM in deta model is conducted using Abaqus software combined with the UMAT subroutine in to take advantage of its built-in nonlinear solver, employing the Newton-Raphson rithm and some other modules.The rest of this paper is organized as follows: In S 2, we introduce the unified phase field theory and the PF-CZM.The approach and merical example of the optimal softening law are described in Section 3.Then, in S 4, we investigate qualitatively and quantitatively the fracture processes and the me cal responses of concrete beams with different features under three-point bending on the PF-CZM with optimal softening law.Finally, a summary of this work and suggestions are given in Section 5.

Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area  arr along the -axis, assume that the rod is fully cracked at  = 0, as shown in Figure auxiliary field variable () ∈ [0, 1] can be introduced to describe the sharp crack ogy; then, we can get Equation (2) based on the variational principle.The crack phas () located at  = 0 can be derived by (3); see more details in [1,18].
where A denotes the fracture surface area, and the crack surface density function is M; the aggregates were built using a random generating and packing algorithm.ang et al. [39] proposed a new mesoscale modeling method in which the PF-CZM is bined with the cell-based smoothed finite element method (CSFEM) to improve comtational efficiency.However, it is not easy to extend these methods to engineering apcations directly.First of all, most studies focus on small specimens, and the distribution both aggregates and random fields is assumed to be ideal.Moreover, the extensive nte Carlo simulations (MCS) cost a lot, even for the RVEs.Additionally, the determition of the micromechanical parameters usually needs extra discussion, and all these tors may lead to an unreal deviation.
Another promising strategy is offered by peridynamics (PD), which substitutes inrodifferential motion equations for partial differential equations (PDEs).These equans are always applicable whether the displacement field is continuous or discontinus.PD has been drawing a lot of attention for its simple and powerful ability to capture plicated failure processes.Li et al. [40] proposed a mesoscale PD for modeling the cking process in concrete.To properly predict cohesive crack propagation, Yang et al. ] proposed a PD-CZM and an approach combining PD and the finite element method M).See [42] for a comparative review between PD and PF models.In this work, the parison between the PF-CZM and the modified PDs will be provided as the second merical example.
It can be seen that most of the existing studies focus mainly on the realization of the lure simulation of quasi-brittle materials.However, there have been a few studies on effects of softening properties in numerical simulations.One key point is the inherent wbacks of the various models.Inspired by the work by Wang [43] and Wu [44] and excellent properties (e.g., no mesh bias, accuracy and insensitivity to the length scale rameter, etc.) of the PF-CZM, it is thus possible to focus on the softening laws in the urate analysis of quasi-brittle materials.Specifically, this work is aimed at proposing a ple new approach to modeling the fracture and damage properties of quasi-brittle maials more accurately and efficiently.First, the optimal softening law is determined by paring the numerical results using several softening laws in the framework of the PF-M with the experimental results.After that, several examples are investigated with the ve approach to validate and extend the application scope of the PF-CZM in detail.The del is conducted using Abaqus software combined with the UMAT subroutine in order take advantage of its built-in nonlinear solver, employing the Newton-Raphson algom and some other modules.The rest of this paper is organized as follows: In Section e introduce the unified phase field theory and the PF-CZM.The approach and a nurical example of the optimal softening law are described in Section 3.Then, in Section e investigate qualitatively and quantitatively the fracture processes and the mechaniresponses of concrete beams with different features under three-point bending based the PF-CZM with optimal softening law.Finally, a summary of this work and some gestions are given in Section 5.

. Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area  arranged ng the -axis, assume that the rod is fully cracked at  = 0, as shown in Figure 1.An xiliary field variable () ∈ [0, 1] can be introduced to describe the sharp crack topoly; then, we can get Equation (2) based on the variational principle.The crack phase field ) located at  = 0 can be derived by (3); see more details in [1,18].α(φ)dφ is the scaling parameter and 0 is a length scale parameter.
where  denotes the fracture surface area, and the crack surface density function is () ≔  As demonstrated by Wu [1,25], some conventional models can be obtained with different  ∈ [0, 2] of the crack geometry function () In this paper, we intend to use the WN model, i.e.,  = 2; thus, we have () = 2 −  2 ,  0 = , () =  As demonstrated by Wu [1,25], some conventional models can be obtained with different ξ ∈ [0, 2] of the crack geometry function α(φ) In this paper, we intend to use the WN model, i.e., ξ = 2; thus, we have CZM; the aggregates were built using a random generating and packing algorithm.Huang et al. [39] proposed a new mesoscale modeling method in which the PF-CZM is combined with the cell-based smoothed finite element method (CSFEM) to improve computational efficiency.However, it is not easy to extend these methods to engineering applications directly.First of all, most studies focus on small specimens, and the distribution of both aggregates and random fields is assumed to be ideal.Moreover, the extensive Monte Carlo simulations (MCS) cost a lot, even for the RVEs.Additionally, the determination of the micromechanical parameters usually needs extra discussion, and all these factors may lead to an unreal deviation.Another promising strategy is offered by peridynamics (PD), which substitutes integrodifferential motion equations for partial differential equations (PDEs).These equations are always applicable whether the displacement field is continuous or discontinuous.PD has been drawing a lot of attention for its simple and powerful ability to capture complicated failure processes.Li et al. [40] proposed a mesoscale PD for modeling the cracking process in concrete.To properly predict cohesive crack propagation, Yang et al. [41] proposed a PD-CZM and an approach combining PD and the finite element method (FEM).See [42] for a comparative review between PD and PF models.In this work, the comparison between the PF-CZM and the modified PDs will be provided as the second numerical example.
It can be seen that most of the existing studies focus mainly on the realization of the failure simulation of quasi-brittle materials.However, there have been a few studies on the effects of softening properties in numerical simulations.One key point is the inherent drawbacks of the various models.Inspired by the work by Wang [43] and Wu [44] and the excellent properties (e.g., no mesh bias, accuracy and insensitivity to the length scale parameter, etc.) of the PF-CZM, it is thus possible to focus on the softening laws in the accurate analysis of quasi-brittle materials.Specifically, this work is aimed at proposing a simple new approach to modeling the fracture and damage properties of quasi-brittle materials more accurately and efficiently.First, the optimal softening law is determined by comparing the numerical results using several softening laws in the framework of the PF-CZM with the experimental results.After that, several examples are investigated with the above approach to validate and extend the application scope of the PF-CZM in detail.The model is conducted using Abaqus software combined with the UMAT subroutine in order to take advantage of its built-in nonlinear solver, employing the Newton-Raphson algorithm and some other modules.The rest of this paper is organized as follows: In Section 2, we introduce the unified phase field theory and the PF-CZM.The approach and a numerical example of the optimal softening law are described in Section 3.Then, in Section 4, we investigate qualitatively and quantitatively the fracture processes and the mechanical responses of concrete beams with different features under three-point bending based on the PF-CZM with optimal softening law.Finally, a summary of this work and some suggestions are given in Section 5.

Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area  arranged along the -axis, assume that the rod is fully cracked at  = 0, as shown in Figure 1.An auxiliary field variable () ∈ [0, 1] can be introduced to describe the sharp crack topology; then, we can get Equation (2) based on the variational principle.The crack phase field () located at  = 0 can be derived by (3); see more details in [1,18].
Extending the 1D crack surface density function M; the aggregates were built using a random generating and packing algorithm.ang et al. [39] proposed a new mesoscale modeling method in which the PF-CZM is bined with the cell-based smoothed finite element method (CSFEM) to improve comtational efficiency.However, it is not easy to extend these methods to engineering apcations directly.First of all, most studies focus on small specimens, and the distribution both aggregates and random fields is assumed to be ideal.Moreover, the extensive nte Carlo simulations (MCS) cost a lot, even for the RVEs.Additionally, the determition of the micromechanical parameters usually needs extra discussion, and all these tors may lead to an unreal deviation.
Another promising strategy is offered by peridynamics (PD), which substitutes inrodifferential motion equations for partial differential equations (PDEs).These equans are always applicable whether the displacement field is continuous or discontinus.PD has been drawing a lot of attention for its simple and powerful ability to capture plicated failure processes.Li et al. [40] proposed a mesoscale PD for modeling the cking process in concrete.To properly predict cohesive crack propagation, Yang et al. ] proposed a PD-CZM and an approach combining PD and the finite element method M).See [42] for a comparative review between PD and PF models.In this work, the parison between the PF-CZM and the modified PDs will be provided as the second merical example.
It can be seen that most of the existing studies focus mainly on the realization of the lure simulation of quasi-brittle materials.However, there have been a few studies on effects of softening properties in numerical simulations.One key point is the inherent wbacks of the various models.Inspired by the work by Wang [43] and Wu [44] and excellent properties (e.g., no mesh bias, accuracy and insensitivity to the length scale rameter, etc.) of the PF-CZM, it is thus possible to focus on the softening laws in the urate analysis of quasi-brittle materials.Specifically, this work is aimed at proposing a ple new approach to modeling the fracture and damage properties of quasi-brittle maials more accurately and efficiently.First, the optimal softening law is determined by paring the numerical results using several softening laws in the framework of the PF-M with the experimental results.After that, several examples are investigated with the ve approach to validate and extend the application scope of the PF-CZM in detail.The del is conducted using Abaqus software combined with the UMAT subroutine in order take advantage of its built-in nonlinear solver, employing the Newton-Raphson algom and some other modules.The rest of this paper is organized as follows: In Section e introduce the unified phase field theory and the PF-CZM.The approach and a nurical example of the optimal softening law are described in Section 3.Then, in Section e investigate qualitatively and quantitatively the fracture processes and the mechaniresponses of concrete beams with different features under three-point bending based the PF-CZM with optimal softening law.Finally, a summary of this work and some gestions are given in Section 5.

. Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area  arranged ng the -axis, assume that the rod is fully cracked at  = 0, as shown in Figure 1.An xiliary field variable () ∈ [0, 1] can be introduced to describe the sharp crack topoly; then, we can get Equation (2) based on the variational principle.The crack phase field ) located at  = 0 can be derived by (3); see more details in [1,18].

EER REVIEW 3 of 19
CZM; the aggregates were built using a random generating and packing algorithm.Huang et al. [39] proposed a new mesoscale modeling method in which the PF-CZM is combined with the cell-based smoothed finite element method (CSFEM) to improve computational efficiency.However, it is not easy to extend these methods to engineering applications directly.First of all, most studies focus on small specimens, and the distribution of both aggregates and random fields is assumed to be ideal.Moreover, the extensive Monte Carlo simulations (MCS) cost a lot, even for the RVEs.Additionally, the determination of the micromechanical parameters usually needs extra discussion, and all these factors may lead to an unreal deviation.
Another promising strategy is offered by peridynamics (PD), which substitutes integrodifferential motion equations for partial differential equations (PDEs).These equations are always applicable whether the displacement field is continuous or discontinuous.PD has been drawing a lot of attention for its simple and powerful ability to capture complicated failure processes.Li et al. [40] proposed a mesoscale PD for modeling the cracking process in concrete.To properly predict cohesive crack propagation, Yang et al. [41] proposed a PD-CZM and an approach combining PD and the finite element method (FEM).See [42] for a comparative review between PD and PF models.In this work, the comparison between the PF-CZM and the modified PDs will be provided as the second numerical example.
It can be seen that most of the existing studies focus mainly on the realization of the failure simulation of quasi-brittle materials.However, there have been a few studies on the effects of softening properties in numerical simulations.One key point is the inherent drawbacks of the various models.Inspired by the work by Wang [43] and Wu [44] and the excellent properties (e.g., no mesh bias, accuracy and insensitivity to the length scale parameter, etc.) of the PF-CZM, it is thus possible to focus on the softening laws in the accurate analysis of quasi-brittle materials.Specifically, this work is aimed at proposing a simple new approach to modeling the fracture and damage properties of quasi-brittle materials more accurately and efficiently.First, the optimal softening law is determined by comparing the numerical results using several softening laws in the framework of the PF-CZM with the experimental results.After that, several examples are investigated with the above approach to validate and extend the application scope of the PF-CZM in detail.The model is conducted using Abaqus software combined with the UMAT subroutine in order to take advantage of its built-in nonlinear solver, employing the Newton-Raphson algorithm and some other modules.The rest of this paper is organized as follows: In Section 2, we introduce the unified phase field theory and the PF-CZM.The approach and a numerical example of the optimal softening law are described in Section 3.Then, in Section 4, we investigate qualitatively and quantitatively the fracture processes and the mechanical responses of concrete beams with different features under three-point bending based on the PF-CZM with optimal softening law.Finally, a summary of this work and some suggestions are given in Section 5.

Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area  arranged along the -axis, assume that the rod is fully cracked at  = 0, as shown in Figure 1.An auxiliary field variable () ∈ [0, 1] can be introduced to describe the sharp crack topology; then, we can get Equation (2) based on the variational principle.The crack phase field () located at  = 0 can be derived by (3); see more details in [1,18].
Two-dimensional sharp and diffusive crack topology is shown in Figure 2. The phase field model diffuses sharp crack Γ into crack band B ⊆ Ω with finite scale 0 > 0. ∂B is the external boundary of the crack band B, The outward unit normal vector of ∂B is denoted as n B .The body force is denoted as b * .The regularized surface energy can be obtained by Equation (8).
suggestions are given in Section 5.

Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area  arranged along the -axis, assume that the rod is fully cracked at  = 0, as shown in Figure 1.An auxiliary field variable () ∈ [0, 1] can be introduced to describe the sharp crack topology; then, we can get Equation (2) based on the variational principle.The crack phase field () located at  = 0 can be derived by (3); see more details in [1,18].
suggestions are given in Section 5.

Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area along the -axis, assume that the rod is fully cracked at  = 0, as shown in auxiliary field variable () ∈ [0, 1] can be introduced to describe the sharp ogy; then, we can get Equation (2) based on the variational principle.The crac () located at  = 0 can be derived by (3); see more details in [1,18].
where G c is the critical energy release rate.

𝛾(𝜙
Two-dimensional sharp and diffusive crack topology is shown in Figure 2. The phase field model diffuses sharp crack  into crack band ℬ ⊆  with finite scale ℓ 0 > 0. ∂ℬ is the external boundary of the crack band ℬ, The outward unit normal vector of ∂ℬ is denoted as  ℬ .The body force is denoted as  * .The regularized surface energy can be obtained by Equation (8).
where   is the critical energy release rate.Correspondingly, the bulk energy is : ℂ 0 :  is the elastic strain energy,  and  are the stress tensor and the strain tensor, respectively. 0 and ℂ 0 are the elastic stiffness tensor and the flexibility tensor, respectively.On the energy degradation function () ∈ [0, 1], Wu gave the rational form, as follows, in [1,25,46], which can be degenerated to the energetic degradation function employed in the conventional phase field fracture model [16,18] or the gradient damage model [45], with the relevant parameters specified.More discussion can be found in the following sections.Correspondingly, the bulk energy is However, conventional damage mechanics methods often suffer from mesh bias problems due to the strain localization phenomenon [4,5] unless some additional methods are proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems exist in the smeared/diffusive and discrete methods of fracture mechanics.
In fracture mechanics, there are many methods currently proposed for crack propagation modeling.These methods may be classified into discrete and smeared/diffusive methods.Discrete methods model cracks as real discontinuities.Representative methods are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obtain mesh-independent results [9][10][11].The XFEM requires an explicit representation of the crack and additional criteria to achieve complex crack growth, and its convergence is not robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can be seen as a development of continuum damage mechanics, which treat the cracks as the result of damage accumulation.As the most popular diffuse approach, the phase field (PF) method has recently received more and more attention from researchers for its advantage of not requiring additional criteria to predict crack nucleation and propagation [1,2,14].The method started with the linear-elastic fracture variational principle proposed by Francfort and Marigo [15] (, ) =  +  = (, )d +  d (1) where the total potential energy functional (, ) is assumed to include the surface energy  of the crack and the bulk energy  ; more details will be discussed in Section 2. Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made further extensions to the PF model, making the phase field modeling for failure generalized.PF methods have been demonstrated to be a powerful modeling tool due to their elegant format for the prediction of failure in many engineering materials such as rock [19], asphalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] and shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality between the classical PF model and the nonlocal damage model and proposed a unified phase-field theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM) derived from the above theory has been validated to be able to overcome most of the drawbacks mentioned in the above models, such as mesh bias and length scale sensitivity, and even the size effect can be captured accurately [26,27].See [28,29] for recent progress.
Multiscale methods are also adopted by combining the above approaches for the accurate analysis of quasi-brittle materials [30,31].Bharali et al. [32] developed a two-scale phase field fracture framework and discussed computational homogenization.However, the representative volume element (RVE) concept is not completely in accordance with actual quasi-brittle materials [2,4].A more popular strategy is considering the local multiscale, i.e., only the potential cracking zone is considered with multiphase at the mesoscale [31].In this way, both indirect and direct approaches are often adopted.In the indirect approaches, the critical fracture parameters of the materials, such as tensile strength and fracture energy, are modeled as stochastic variables with random fields so that different mixtures at the mesoscale are modeled implicitly.Recently, Li et al. [33] developed a method combining the PF-CZM and random fields for quasi-brittle materials in a 2D setting; only tensile strength is modeled by random fields.Zhang et al. [34] analyzed the size effect at the mesoscale with the same method.Huang et al. [35] adopted statistical reconstruction algorithms and the Karhunen-Loève expansion technique with microscale X-ray computed tomography (CT) images in order to get the random fields of tensile strength.In addition, Lu et al. [36] generated a random field with a stochastic harmonic in which owever, conventional damage mechanics methods often suffer from mesh bias probems due to the strain localization phenomenon [4,5] unless some additional methods are roposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems exist n the smeared/diffusive and discrete methods of fracture mechanics.In fracture mechanics, there are many methods currently proposed for crack propaation modeling.These methods may be classified into discrete and smeared/diffusive ethods.Discrete methods model cracks as real discontinuities.Representative methods re cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obtain esh-independent results [9][10][11].The XFEM requires an explicit representation of the rack and additional criteria to achieve complex crack growth, and its convergence is not obust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can be een as a development of continuum damage mechanics, which treat the cracks as the esult of damage accumulation.As the most popular diffuse approach, the phase field PF) method has recently received more and more attention from researchers for its adantage of not requiring additional criteria to predict crack nucleation and propagation 1,2,14].The method started with the linear-elastic fracture variational principle proposed y Francfort and Marigo [15] (, ) =  +  = (, )d +  d (1) here the total potential energy functional (, ) is assumed to include the surface enrgy  of the crack and the bulk energy  ; more details will be discussed in Section 2. ubsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made further xtensions to the PF model, making the phase field modeling for failure generalized.PF ethods have been demonstrated to be a powerful modeling tool due to their elegant ormat for the prediction of failure in many engineering materials such as rock [19], ashalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] and hield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality between he classical PF model and the nonlocal damage model and proposed a unified phase-field heory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM) erived from the above theory has been validated to be able to overcome most of the rawbacks mentioned in the above models, such as mesh bias and length scale sensitivity, nd even the size effect can be captured accurately [26,27].See [28,29] for recent progress.
Multiscale methods are also adopted by combining the above approaches for the acurate analysis of quasi-brittle materials [30,31].Bharali et al. [32] developed a two-scale hase field fracture framework and discussed computational homogenization.However, he representative volume element (RVE) concept is not completely in accordance with ctual quasi-brittle materials [2,4].A more popular strategy is considering the local muliscale, i.e., only the potential cracking zone is considered with multiphase at the esoscale [31].In this way, both indirect and direct approaches are often adopted.In the ndirect approaches, the critical fracture parameters of the materials, such as tensile trength and fracture energy, are modeled as stochastic variables with random fields so hat different mixtures at the mesoscale are modeled implicitly.Recently, Li et al. [33] deeloped a method combining the PF-CZM and random fields for quasi-brittle materials in 2D setting; only tensile strength is modeled by random fields.Zhang et al. [34] analyzed he size effect at the mesoscale with the same method.Huang et al. [35] adopted statistical econstruction algorithms and the Karhunen-Loève expansion technique with microscale Two-dimensional sharp and diffusive crack topology is shown in Figure 2. The phase field model diffuses sharp crack  into crack band ℬ ⊆  with finite scale ℓ > 0. ∂ℬ is the external boundary of the crack band ℬ, The outward unit normal vector of ∂ℬ is denoted as  ℬ .The body force is denoted as  * .The regularized surface energy can be obtained by Equation (8).
where  is the critical energy release rate.Correspondingly, the bulk energy is  = (, )d (9) in which (, ) = () () is the bulk energy density,  () = :  :  = : ℂ :  is the elastic strain energy,  and  are the stress tensor and the strain tensor, respectively. and ℂ are the elastic stiffness tensor and the flexibility tensor, respectively.On the energy degradation function () ∈ [0, 1], Wu gave the rational form, as follows, in [1,25,46], which can be degenerated to the energetic degradation function employed in the conventional phase field fracture model [16,18] or the gradient damage model [45], with the relevant parameters specified.More discussion can be found in the following sections. (φ) x FOR PEER REVIEW 2 of 19 However, conventional damage mechanics methods often suffer from mesh bias problems due to the strain localization phenomenon [4,5] unless some additional methods are proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems exist in the smeared/diffusive and discrete methods of fracture mechanics.In fracture mechanics, there are many methods currently proposed for crack propagation modeling.These methods may be classified into discrete and smeared/diffusive methods.Discrete methods model cracks as real discontinuities.Representative methods are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obtain mesh-independent results [9][10][11].The XFEM requires an explicit representation of the crack and additional criteria to achieve complex crack growth, and its convergence is not robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can be seen as a development of continuum damage mechanics, which treat the cracks as the result of damage accumulation.As the most popular diffuse approach, the phase field (PF) method has recently received more and more attention from researchers for its advantage of not requiring additional criteria to predict crack nucleation and propagation [1,2,14].The method started with the linear-elastic fracture variational principle proposed by Francfort and Marigo [15] (, ) =  +  = (, )d +  d (1) where the total potential energy functional (, ) is assumed to include the surface energy  of the crack and the bulk energy  ; more details will be discussed in Section 2. Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made further extensions to the PF model, making the phase field modeling for failure generalized.PF methods have been demonstrated to be a powerful modeling tool due to their elegant format for the prediction of failure in many engineering materials such as rock [19], asphalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] and shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality between the classical PF model and the nonlocal damage model and proposed a unified phase-field theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM) derived from the above theory has been validated to be able to overcome most of the drawbacks mentioned in the above models, such as mesh bias and length scale sensitivity, and even the size effect can be captured accurately [26,27].See [28,29] for recent progress.
Multiscale methods are also adopted by combining the above approaches for the accurate analysis of quasi-brittle materials [30,31].Bharali et al. [32] developed a two-scale phase field fracture framework and discussed computational homogenization.However, the representative volume element (RVE) concept is not completely in accordance with actual quasi-brittle materials [2,4].A more popular strategy is considering the local multiscale, i.e., only the potential cracking zone is considered with multiphase at the mesoscale [31].In this way, both indirect and direct approaches are often adopted.In the indirect approaches, the critical fracture parameters of the materials, such as tensile strength and fracture energy, are modeled as stochastic variables with random fields so that different mixtures at the mesoscale are modeled implicitly.Recently, Li et al. [33] developed a method combining the PF-CZM and random fields for quasi-brittle materials in a 2D setting; only tensile strength is modeled by random fields.Zhang et al. [34] analyzed the size effect at the mesoscale with the same method.Huang et al. [35] adopted statistical reconstruction algorithms and the Karhunen-Loève expansion technique with microscale 0 ( ) is the bulk energy density, However, conventional damage mechanics methods often suffer from mesh bias pr lems due to the strain localization phenomenon [4,5] unless some additional methods proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems e in the smeared/diffusive and discrete methods of fracture mechanics.In fracture mechanics, there are many methods currently proposed for crack pro gation modeling.These methods may be classified into discrete and smeared/diffu methods.Discrete methods model cracks as real discontinuities.Representative meth are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily ob mesh-independent results [9][10][11].The XFEM requires an explicit representation of crack and additional criteria to achieve complex crack growth, and its convergence is robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can seen as a development of continuum damage mechanics, which treat the cracks as result of damage accumulation.As the most popular diffuse approach, the phase f (PF) method has recently received more and more attention from researchers for its vantage of not requiring additional criteria to predict crack nucleation and propaga [1,2,14].The method started with the linear-elastic fracture variational principle propo by Francfort and Marigo [15] (, ) =  +  = (, )d +  d where the total potential energy functional (, ) is assumed to include the surface ergy  of the crack and the bulk energy  ; more details will be discussed in Sectio Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made fur extensions to the PF model, making the phase field modeling for failure generalized methods have been demonstrated to be a powerful modeling tool due to their eleg format for the prediction of failure in many engineering materials such as rock [19], phalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality betw the classical PF model and the nonlocal damage model and proposed a unified phase-f theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZ derived from the above theory has been validated to be able to overcome most of drawbacks mentioned in the above models, such as mesh bias and length scale sensitiv and even the size effect can be captured accurately [26,27].See [28,29] for recent progr Multiscale methods are also adopted by combining the above approaches for the curate analysis of quasi-brittle materials [30,31].Bharali et al. [32] developed a two-s phase field fracture framework and discussed computational homogenization.Howe the representative volume element (RVE) concept is not completely in accordance w actual quasi-brittle materials [2,4].A more popular strategy is considering the local m tiscale, i.e., only the potential cracking zone is considered with multiphase at mesoscale [31].In this way, both indirect and direct approaches are often adopted.In indirect approaches, the critical fracture parameters of the materials, such as ten strength and fracture energy, are modeled as stochastic variables with random field that different mixtures at the mesoscale are modeled implicitly.Recently, Li et al. [33] veloped a method combining the PF-CZM and random fields for quasi-brittle material a 2D setting; only tensile strength is modeled by random fields.Zhang et al. [34] analy the size effect at the mesoscale with the same method.Huang et al. [35] adopted statist reconstruction algorithms and the Karhunen-Loève expansion technique with micros 0 ( ) = 1 2 : E 0 : = 1 2 σ : C 0 : σ is the elastic strain energy, σ and are the stress tensor and the strain tensor, respectively.E 0 and C 0 are the elastic stiffness tensor and the flexibility tensor, respectively.On the energy degradation function  Correspondingly, the bulk energy is  = (, )d in which (, ) = () () is the bulk energy density,  () = :  :  = is the elastic strain energy,  and  are the stress tensor and the strain tensor, tively. and ℂ are the elastic stiffness tensor and the flexibility tensor, resp On the energy degradation function () ∈ [0, 1], Wu gave the rational form, as in [1,25,46], which can be degenerated to the energetic degradation function emp the conventional phase field fracture model [16,18] or the gradient damage mo with the relevant parameters specified.More discussion can be found in the fo sections.
Wu gave the rational form, as follows, in [1,25,46], which can be degenerated to the energetic degradation function employed in the conventional phase field fracture model [16,18] or the gradient damage model [45], with the relevant parameters specified.More discussion can be found in the following sections.Correspondingly, the bulk energy is  = (, )d (9) in which (, ) = () () is the bulk energy density,  () = :  :  = : ℂ :  is the elastic strain energy,  and  are the stress tensor and the strain tensor, respec-tively. and ℂ are the elastic stiffness tensor and the flexibility tensor, respectively.On the energy degradation function () ∈ [0, 1], Wu gave the rational form, as follows, in [1,25,46], which can be degenerated to the energetic degradation function employed in the conventional phase field fracture model [16,18] or the gradient damage model [45], with the relevant parameters specified.More discussion can be found in the following sections.

Principle of Virtual Work and Governing Equations
By using the results from the previous section, the total potential energy functional can be written as  Correspondingly, the bulk energy is  = (, )d in which (, ) = () () is the bulk energy density,  () = :  :  = : ℂ is the elastic strain energy,  and  are the stress tensor and the strain tensor, resp tively. and ℂ are the elastic stiffness tensor and the flexibility tensor, respectiv On the energy degradation function () ∈ [0, 1], Wu gave the rational form, as follo in [1,25,46], which can be degenerated to the energetic degradation function employed the conventional phase field fracture model [16,18] or the gradient damage model [ with the relevant parameters specified.More discussion can be found in the follow sections.
However, conventional damage mechanics methods often suffer from mesh bias problems due to the strain localization phenomenon [4,5] unless some additional methods are proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems exist in the smeared/diffusive and discrete methods of fracture mechanics.In fracture mechanics, there are many methods currently proposed for crack propagation modeling.These methods may be classified into discrete and smeared/diffusive methods.Discrete methods model cracks as real discontinuities.Representative methods are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obtain mesh-independent results [9][10][11].The XFEM requires an explicit representation of the crack and additional criteria to achieve complex crack growth, and its convergence is not robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can be seen as a development of continuum damage mechanics, which treat the cracks as the result of damage accumulation.As the most popular diffuse approach, the phase field (PF) method has recently received more and more attention from researchers for its advantage of not requiring additional criteria to predict crack nucleation and propagation [1,2,14].The method started with the linear-elastic fracture variational principle proposed by Francfort and Marigo [15] (, ) =  +  = (, )d +  d where the total potential energy functional (, ) is assumed to include the surface energy  of the crack and the bulk energy  ; more details will be discussed in Section 2. Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made further extensions to the PF model, making the phase field modeling for failure generalized.PF methods have been demonstrated to be a powerful modeling tool due to their elegant format for the prediction of failure in many engineering materials such as rock [19], asphalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] and shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality between the classical PF model and the nonlocal damage model and proposed a unified phase-field theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM) Appl.Sci.2022, 12, x FOR PEER REVIEW CZM; the aggregates were built using a random generating and packing a Huang et al. [39] proposed a new mesoscale modeling method in which the P combined with the cell-based smoothed finite element method (CSFEM) to impr putational efficiency.However, it is not easy to extend these methods to engine plications directly.First of all, most studies focus on small specimens, and the di of both aggregates and random fields is assumed to be ideal.Moreover, the Monte Carlo simulations (MCS) cost a lot, even for the RVEs.Additionally, the nation of the micromechanical parameters usually needs extra discussion, and factors may lead to an unreal deviation.
Another promising strategy is offered by peridynamics (PD), which subs tegrodifferential motion equations for partial differential equations (PDEs).Th tions are always applicable whether the displacement field is continuous or di ous.PD has been drawing a lot of attention for its simple and powerful ability t complicated failure processes.Li et al. [40] proposed a mesoscale PD for mod cracking process in concrete.To properly predict cohesive crack propagation, Y [41] proposed a PD-CZM and an approach combining PD and the finite elemen (FEM).See [42] for a comparative review between PD and PF models.In this comparison between the PF-CZM and the modified PDs will be provided as th numerical example.
It can be seen that most of the existing studies focus mainly on the realizat failure simulation of quasi-brittle materials.However, there have been a few s the effects of softening properties in numerical simulations.One key point is the drawbacks of the various models.Inspired by the work by Wang [43] and Wu the excellent properties (e.g., no mesh bias, accuracy and insensitivity to the len parameter, etc.) of the PF-CZM, it is thus possible to focus on the softening la accurate analysis of quasi-brittle materials.Specifically, this work is aimed at pr simple new approach to modeling the fracture and damage properties of quasi-b terials more accurately and efficiently.First, the optimal softening law is deter comparing the numerical results using several softening laws in the framework CZM with the experimental results.After that, several examples are investigated above approach to validate and extend the application scope of the PF-CZM in d model is conducted using Abaqus software combined with the UMAT subroutin to take advantage of its built-in nonlinear solver, employing the Newton-Raph rithm and some other modules.The rest of this paper is organized as follows: I 2, we introduce the unified phase field theory and the PF-CZM.The approach merical example of the optimal softening law are described in Section 3.Then, 4, we investigate qualitatively and quantitatively the fracture processes and the cal responses of concrete beams with different features under three-point bend on the PF-CZM with optimal softening law.Finally, a summary of this work suggestions are given in Section 5.

Phase Field Representation of Crack Surface
Considering an infinitely long circular rod with a cross-sectional area  along the -axis, assume that the rod is fully cracked at  = 0, as shown in Fig auxiliary field variable () ∈ [0, 1] can be introduced to describe the sharp cra ogy; then, we can get Equation (2) based on the variational principle.The crack p () located at  = 0 can be derived by (3); see more details in [1,18].
Here, the compact tensor is not used.For a linear-elastic material, Appl.Sci.2022, 12, x FOR PEER REVIEW However, conventional damage mechanics methods often suffer from mesh bia lems due to the strain localization phenomenon [4,5] unless some additional meth proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problem in the smeared/diffusive and discrete methods of fracture mechanics.In fracture mechanics, there are many methods currently proposed for crack gation modeling.These methods may be classified into discrete and smeared/di methods.Discrete methods model cracks as real discontinuities.Representative m are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily mesh-independent results [9][10][11].The XFEM requires an explicit representation crack and additional criteria to achieve complex crack growth, and its convergenc robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach seen as a development of continuum damage mechanics, which treat the cracks result of damage accumulation.As the most popular diffuse approach, the phas (PF) method has recently received more and more attention from researchers for vantage of not requiring additional criteria to predict crack nucleation and propa [1,2,14].The method started with the linear-elastic fracture variational principle pr by Francfort and Marigo [15] (, ) =  +  = (, )d +  d where the total potential energy functional (, ) is assumed to include the surf ergy  of the crack and the bulk energy  ; more details will be discussed in Se Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made extensions to the PF model, making the phase field modeling for failure generaliz methods have been demonstrated to be a powerful modeling tool due to their format for the prediction of failure in many engineering materials such as rock [ phalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [2 shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality b 0 ( ) = 1 2 λε kk ε ll + µε ij ε ij , with λ and µ being the Láme constants; i, j, k, l = 1, . . ., d, with d being dimensional numbers.For small strains, i.e., ε ij = 1 2 ∂u i /∂x j + ∂u j /∂x i .The governing equations are obtained as follows: we first define the variation of the external work increment as On the other hand, the variation of the internal energy increment is given by which for the case of ( 12 lems due to the strain localization phenomenon [4,5] unless some additional methods a proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems ex in the smeared/diffusive and discrete methods of fracture mechanics. In fracture mechanics, there are many methods currently proposed for crack prop gation modeling.These methods may be classified into discrete and smeared/diffusi methods.Discrete methods model cracks as real discontinuities.Representative metho are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obta mesh-independent results [9][10][11].The XFEM requires an explicit representation of t crack and additional criteria to achieve complex crack growth, and its convergence is n robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can seen as a development of continuum damage mechanics, which treat the cracks as t result of damage accumulation.As the most popular diffuse approach, the phase fie (PF) method has recently received more and more attention from researchers for its a vantage of not requiring additional criteria to predict crack nucleation and propagatio [1,2,14].The method started with the linear-elastic fracture variational principle propose by Francfort and Marigo [15] (, ) =  +  = (, )d +  d ( where the total potential energy functional (, ) is assumed to include the surface e ergy  of the crack and the bulk energy  ; more details will be discussed in Section Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made furth extensions to the PF model, making the phase field modeling for failure generalized.P methods have been demonstrated to be a powerful modeling tool due to their elega format for the prediction of failure in many engineering materials such as rock [19], a phalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] an shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality betwe the classical PF model and the nonlocal damage model and proposed a unified phase-fie theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM derived from the above theory has been validated to be able to overcome most of t drawbacks mentioned in the above models, such as mesh bias and length scale sensitivit and even the size effect can be captured accurately [26,27].See [28,29] for recent progres Multiscale methods are also adopted by combining the above approaches for the a curate analysis of quasi-brittle materials [30,31].Bharali et al. [32] developed a two-sca phase field fracture framework and discussed computational homogenization.Howeve the representative volume element (RVE) concept is not completely in accordance wi actual quasi-brittle materials [2,4].A more popular strategy is considering the local mu tiscale, i.e., only the potential cracking zone is considered with multiphase at t mesoscale [31].In this way, both indirect and direct approaches are often adopted.In t indirect approaches, the critical fracture parameters of the materials, such as tens strength and fracture energy, are modeled as stochastic variables with random fields that different mixtures at the mesoscale are modeled implicitly.Recently, Li et al. [33] d veloped a method combining the PF-CZM and random fields for quasi-brittle materials a 2D setting; only tensile strength is modeled by random fields.Zhang et al. [34] analyze the size effect at the mesoscale with the same method.Huang et al. [35] adopted statistic reconstruction algorithms and the Karhunen-Loève expansion technique with microsca X-ray computed tomography (CT) images in order to get the random fields of tens strength.In addition, Lu et al. [36] generated a random field with a stochastic harmon function (SHF).On the other hand, various phases at the mesoscale are generally tak into account explicitly in the direct approaches.Xia et al. [37] investigated the mechanic behaviors of concrete in a 2D mesoscale setting, with the aggregates generated using random sequential addition (RSA) method.Li et al. [38] developed a 3D mesosca method for the damage and fracture simulation of quasi-brittle materials with the P By appropriate constant transformation with divergence theorem, it can be found that the above expression is equivalent to where  is the critical energy release rate.for ∂ , respectively.The outward unit normal vector of external boundary ∂ is denoted as .
Correspondingly, the bulk energy is  = (, )d (9) in which (, ) = () () is the bulk energy density,  () = :  :  = : ℂ :  is the elastic strain energy,  and  are the stress tensor and the strain tensor, respectively. and ℂ are the elastic stiffness tensor and the flexibility tensor, respectively.On the energy degradation function () ∈ [0, 1], Wu gave the rational form, as follows, in [1,25,46], which can be degenerated to the energetic degradation function employed in the conventional phase field fracture model [16,18] or the gradient damage model [45], with the relevant parameters specified.More discussion can be found in the following sections.However, conventional damage mechanics methods often suffer from mesh bias problems due to the strain localization phenomenon [4,5] unless some additional methods are proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems exist in the smeared/diffusive and discrete methods of fracture mechanics.In fracture mechanics, there are many methods currently proposed for crack propagation modeling.These methods may be classified into discrete and smeared/diffusive methods.Discrete methods model cracks as real discontinuities.Representative methods are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obtain mesh-independent results [9][10][11].The XFEM requires an explicit representation of the crack and additional criteria to achieve complex crack growth, and its convergence is not robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can be seen as a development of continuum damage mechanics, which treat the cracks as the result of damage accumulation.As the most popular diffuse approach, the phase field (PF) method has recently received more and more attention from researchers for its advantage of not requiring additional criteria to predict crack nucleation and propagation [1,2,14].The method started with the linear-elastic fracture variational principle proposed by Francfort and Marigo [15] (, ) =  +  = (, )d +  d (1) where the total potential energy functional (, ) is assumed to include the surface energy  of the crack and the bulk energy  ; more details will be discussed in Section 2. Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made further extensions to the PF model, making the phase field modeling for failure generalized.PF methods have been demonstrated to be a powerful modeling tool due to their elegant format for the prediction of failure in many engineering materials such as rock [19], asphalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] and shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality between the classical PF model and the nonlocal damage model and proposed a unified phase-field theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM) derived from the above theory has been validated to be able to overcome most of the drawbacks mentioned in the above models, such as mesh bias and length scale sensitivity, and even the size effect can be captured accurately [26,27].See [28,29] for recent progress.

𝜔(𝜙
Multiscale methods are also adopted by combining the above approaches for the accurate analysis of quasi-brittle materials [30,31].Bharali et al. [32] developed a two-scale phase field fracture framework and discussed computational homogenization.However, the representative volume element (RVE) concept is not completely in accordance with actual quasi-brittle materials [2,4].A more popular strategy is considering the local multiscale, i.e., only the potential cracking zone is considered with multiphase at the mesoscale [31].In this way, both indirect and direct approaches are often adopted.In the indirect approaches, the critical fracture parameters of the materials, such as tensile strength and fracture energy, are modeled as stochastic variables with random fields so that different mixtures at the mesoscale are modeled implicitly.Recently, Li et al. [33] developed a method combining the PF-CZM and random fields for quasi-brittle materials in a 2D setting; only tensile strength is modeled by random fields.Zhang et al. [34] analyzed the size effect at the mesoscale with the same method.Huang et al. [35] adopted statistical reconstruction algorithms and the Karhunen-Loève expansion technique with microscale X-ray computed tomography (CT) images in order to get the random fields of tensile strength.In addition, Lu et al. [36] generated a random field with a stochastic harmonic function (SHF).On the other hand, various phases at the mesoscale are generally taken into account explicitly in the direct approaches.Xia et al. [37] investigated the mechanical behaviors of concrete in a 2D mesoscale setting, with the aggregates generated using a random sequential addition (RSA) method.Li et al. [38] developed a 3D mesoscale method for the damage and fracture simulation of quasi-brittle materials with the PF- Combine the terms in ( 13) and ( 16    Wu gave the rational form, as follo in [1,25,46], which can be degenerated to the energetic degradation function employe the conventional phase field fracture model [16,18] or the gradient damage model with the relevant parameters specified.More discussion can be found in the follow sections.where  is the critical energy release rate.Correspondingly, the bulk energy is  = (, )d in which (, ) = () () is the bulk energy density,  () = is the elastic strain energy,  and  are the stress tensor and the st tively. and ℂ are the elastic stiffness tensor and the flexibility t On the energy degradation function () ∈ [0, 1], Wu gave the ration in [1,25,46], which can be degenerated to the energetic degradation fu the conventional phase field fracture model [16,18] or the gradient d with the relevant parameters specified.More discussion can be fou sections.

𝜔(𝜙
However, conventional damage mechanics methods often suffer from mesh bias problems due to the strain localization phenomenon [4,5] unless some additional methods are proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems exist in the smeared/diffusive and discrete methods of fracture mechanics.In fracture mechanics, there are many methods currently proposed for crack propagation modeling.These methods may be classified into discrete and smeared/diffusive methods.Discrete methods model cracks as real discontinuities.Representative methods are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obtain mesh-independent results [9][10][11].The XFEM requires an explicit representation of the crack and additional criteria to achieve complex crack growth, and its convergence is not robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can be seen as a development of continuum damage mechanics, which treat the cracks as the result of damage accumulation.As the most popular diffuse approach, the phase field (PF) method has recently received more and more attention from researchers for its advantage of not requiring additional criteria to predict crack nucleation and propagation [1,2,14].The method started with the linear-elastic fracture variational principle proposed by Francfort and Marigo [15] (, ) =  +  = (, )d +  d where the total potential energy functional (, ) is assumed to include the surface energy  of the crack and the bulk energy  ; more details will be discussed in Section 2. Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made further extensions to the PF model, making the phase field modeling for failure generalized.PF methods have been demonstrated to be a powerful modeling tool due to their elegant format for the prediction of failure in many engineering materials such as rock [19], asphalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] and shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality between the classical PF model and the nonlocal damage model and proposed a unified phase-field theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM) derived from the above theory has been validated to be able to overcome most of the drawbacks mentioned in the above models, such as mesh bias and length scale sensitivity, and even the size effect can be captured accurately [26,27].See [28,29] for recent progress.
Multiscale methods are also adopted by combining the above approaches for the accurate analysis of quasi-brittle materials [30,31].Bharali et al. [32] developed a two-scale phase field fracture framework and discussed computational homogenization.However, the representative volume element (RVE) concept is not completely in accordance with actual quasi-brittle materials [2,4].A more popular strategy is considering the local multiscale, i.e., only the potential cracking zone is considered with multiphase at the mesoscale [31].In this way, both indirect and direct approaches are often adopted.In the indirect approaches, the critical fracture parameters of the materials, such as tensile strength and fracture energy, are modeled as stochastic variables with random fields so that different mixtures at the mesoscale are modeled implicitly.Recently, Li et al. [33] developed a method combining the PF-CZM and random fields for quasi-brittle materials in a 2D setting; only tensile strength is modeled by random fields.Zhang et al. [34] analyzed the size effect at the mesoscale with the same method.Huang et al. [35] adopted statistical The above shows a coupled system consisting of the modified stress equilibrium (17) and the phase field evolution (20); ( 18), ( 19) and ( 21) are the boundary conditions for the modified stress equilibrium and the phase field evolution equations, respectively.
For quasi-brittle materials, there is significant tensile-compression anisotropic behavior because their compressive strength is much greater than their tensile strength.To avoid the propagation of the crack under compression, Amor [17], Miehe [18] and other scholars proposed a volumetric-deviatoric split, spectral decomposition and other methods to overcome the above problem in the classical phase field fracture model.The non-variationally consistent formula (22) introduced by Wu [47] is used in this paper to replace the elastic strain energy EVIEW 2 of 19 However, conventional damage mechanics methods often suffer from mesh bias problems due to the strain localization phenomenon [4,5] unless some additional methods are proposed, as in the nonlocal and gradient-enhanced models [6][7][8].Similar problems exist in the smeared/diffusive and discrete methods of fracture mechanics.In fracture mechanics, there are many methods currently proposed for crack propagation modeling.These methods may be classified into discrete and smeared/diffusive methods.Discrete methods model cracks as real discontinuities.Representative methods are cohesive zone models(CZMs), XFEMs, etc.However, the CZM cannot easily obtain mesh-independent results [9][10][11].The XFEM requires an explicit representation of the crack and additional criteria to achieve complex crack growth, and its convergence is not robust, especially in a 3D setting [12,13].On the other hand, the diffusive approach can be seen as a development of continuum damage mechanics, which treat the cracks as the result of damage accumulation.As the most popular diffuse approach, the phase field (PF) method has recently received more and more attention from researchers for its advantage of not requiring additional criteria to predict crack nucleation and propagation [1,2,14].The method started with the linear-elastic fracture variational principle proposed by Francfort and Marigo [15] (, ) =  +  = (, )d +  d where the total potential energy functional (, ) is assumed to include the surface energy  of the crack and the bulk energy  ; more details will be discussed in Section 2. Subsequently, Bourdin [16], Amor [17], Miehe [18] and some other scholars made further extensions to the PF model, making the phase field modeling for failure generalized.PF methods have been demonstrated to be a powerful modeling tool due to their elegant format for the prediction of failure in many engineering materials such as rock [19], asphalt concrete [20], fiber-reinforced composites [21,22], piezoelectric composites [23] and shield tunnel lining [24], to cite a few.In 2017, Wu discovered the commonality between the classical PF model and the nonlocal damage model and proposed a unified phase-field theory that is suitable for both brittle and cohesive failure [1,25].The model (e.g., PF-CZM) derived from the above theory has been validated to be able to overcome most of the drawbacks mentioned in the above models, such as mesh bias and length scale sensitivity, and even the size effect can be captured accurately [26,27].See [28,29] for recent progress.0 ( ), in order to consider the different mechanical behaviors under tensile and compressive stress.
where E 0 is Young's modulus; the Macaulay bracket • is defined as x = max(x, 0), and σ 1 is the maximum principal stress of the stress tensor σ.

Damage Irreversibility
To deal with the irreversibility of phase field evolution ( .φ ≥ 0), we adopt the history field H to prevent crack healing, which can ensure where φ t+∆t represents the phase field value in the current time increment while φ t denotes the value of the previous increment.The history field variable must satisfy the following Kuhn-Tucker conditions [48] Y − H ≤ 0, .
Therefore, the history field variable at the current moment can be taken as Then, the evolution equation of the phase field is  (φ) are the two main functions that define the above phase field theory.According to the work of Wu et al. [1,25], it is thus possible to define an equivalent cohesive zone model explicitly due to the formula g(φ) in the geometric crack function α(φ).
Two-dimensional sharp and diffusive crack topology is shown in Figure 2. The phase field model diffuses sharp crack  into crack band ℬ ⊆  with finite scale ℓ > 0. ∂ℬ is the external boundary of the crack band ℬ, The outward unit normal vector of ∂ℬ is denoted as  ℬ .The body force is denoted as  * . The regularized surface energy can be obtained by Equation (8). =  d =  (, ∇)d =  ℬ (, ∇)d (8 where  is the critical energy release rate.for ∂ , respectively.The outward unit normal vector of external boundary ∂ is denoted as . Correspondingly, the bulk energy is  = (, )d (9) in which (, ) = () () is the bulk energy density,  () = :  :  = : ℂ :  is the elastic strain energy,  and  are the stress tensor and the strain tensor, respectively. and ℂ are the elastic stiffness tensor and the flexibility tensor, respectively.On the energy degradation function () ∈ [0, 1], Wu gave the rational form, as follows, in [1,25,46], which can be degenerated to the energetic degradation function employed in the conventional phase field fracture model [16,18] or the gradient damage model [45], with the relevant parameters specified.More discussion can be found in the following sections.

𝜔(𝜙
in which the exponent p ≥ 2, coefficients a 1 > 0, a 2 , a 3 can be determined by the failure strength f t (uniaxial tensile strength) and the initial slope k 0 and the ultimate crack opening w c of the target softening law σ(w), respectively.
ξ 2 (28) where P(1) = 1 + a 2 + a 2 a 3 .It can be found since the ultimate crack opening displacement is zero when p < 2. This is not in accordance with the features of quasi-brittle failure.Therefore, the exponential parameter in the energy degradation function  Correspondingly, the bulk energy is  = in which (, ) = () () is the bulk e is the elastic strain energy,  and  are the tively. and ℂ are the elastic stiffness te On the energy degradation function () ∈ in [1,25,46], which can be degenerated to the the conventional phase field fracture model with the relevant parameters specified.Mor sections.
for Griffith's characteristic length l ch := E 0 G c / f 2 t .The general softening law used in quasi-brittle failure models of concrete can be obtained by the parameterizing characteristic function (27).In particular, we give the equivalent parameters of Petersson's [49] softening curve for the first time to the best knowledge of the authors:

Computer Implementation
The governing equations of the PF-CZM are often solved with the finite element method [47,48].In this work, a monolithic solution scheme was first adopted, and the model was implemented through the Abaqus user material (UMAT) subroutine.Heat transfer analogy and a subdomain implementation scheme were considered for convenience and efficiency.See more details in [31,48].

Optimal Softening Law
A typical softening law was often assumed for convenience in many studies [39,43,50].However, it is well known that quasi-brittle materials are usually heterogeneous, especially concrete-like materials.It should be noted that we have confined our investigation mainly to concrete in this work due to the complexity of other quasi-brittle materials, such as fiberreinforced concrete, which deserve additional analyses.Due to the excellent properties of the PF-CZM, such as no mesh bias, accuracy and insensitivity to the length scale parameter, it is possible to focus only on the effects of softening laws on the mechanical responses of different concretes [1,25,26].In this section, we propose and validate a simple new approach to determining the optimal softening law in the framework of the PF-CZM.It could be found in the following section that the geometry of the specimen would not affect the softening properties significantly, which means the approach would be suitable for engineering analysis directly once the optimal softening law is determined.It should be noted that complex, large-scale materials and structures are not involved in this work.
A set of notched concrete beams under three-point bending is modeled, as tested by [43] and shown in Figure 3.Only two groups of the test are considered in this work, i.e., C60 and C100.Additionally, there are three specimens in each group.The material properties are shown in Table 1.The maximum tensile stress failure criterion is adopted.It can be easily found in Figure 4 that the predicted curves using the concrete damage plasticity (CDP) model with Petersson's [49] softening law are lower than the experimental results, and the error increases when the strength of the specimen increases.The same tendency could be observed in the P-CMOD curves using the PF-CZM.It is well known that the softening law can reflect the stochastic distribution of mesoscale mixtures in FPZ implicitly.For instance, we only compared the results with three typical softening curves in this paper to determine the optimal softening law.It should be noted that the consideration of other softening laws is naturally based on this approach [26,51].In addition, the length scale parameter is 0 = 2.5 mm, and the plane stress assumption is adopted for 2D problems in this work.Displacement-controlled is used for the acquisition of mechanical response curves, and the prescribed displacement u = 0.8 mm is applied in 80 incremental steps, with each incremental step ∆u fixed at 0.01 mm.
Figure 4 depicts the P-CMOD curves predicted using the PF-CZM with linear, exponential and Petersson's softening curves, together with the experimental results.Since the crack path was not given, here we mainly analyzed the P-CMOD curves.It can easily be found in Figure 4 for both C60 and C100 concrete samples that the P-CMOD curves obtained with the linear softening law are higher than with other softening laws, and the results using the exponential softening law are the minimum.In Figure 4a, all predicted curves are higher than the experimental results.The numerical results are consistent with the experimental data on the whole, except when using the linear softening law.The specific predicted peak loads are listed in Table 2.By comparing the experimental mean peak load and the predicted peak load, we can assume that the exponential softening curve is the optimal softening law for C60 concrete, while for C100 concrete, it is the Petersson softening curve.Furthermore, it can also be found that for C60 concrete, the predicted results are higher than in the experiments in general; this deviation can be modified by adjusting the softening law.A more accurate method of obtaining an accurate softening curve, such as inverse analysis [51,52], will be considered further in the framework of the PF-CZM for analysis.
Appl.Sci.2022, 12, x FOR PEER REVIEW 9 of 19 other softening laws is naturally based on this approach [26,51].In addition, the length scale parameter is ℓ 0 = 2.5 mm, and the plane stress assumption is adopted for 2D problems in this work.Displacement-controlled is used for the acquisition of mechanical response curves, and the prescribed displacement  = 0.8 mm is applied in 80 incremental steps, with each incremental step Δ fixed at 0.01 mm.other softening laws is naturally based on this approach [26,51].In addition, the length scale parameter is ℓ 0 = 2.5 mm, and the plane stress assumption is adopted for 2D problems in this work.Displacement-controlled is used for the acquisition of mechanical response curves, and the prescribed displacement  = 0.8 mm is applied in 80 incremental steps, with each incremental step Δ fixed at 0.01 mm.   Figure 4 depicts the P-CMOD curves predicted using the PF-CZM with linear, exponential and Petersson's softening curves, together with the experimental results.Since the crack path was not given, here we mainly analyzed the P-CMOD curves.It can easily be found in Figure 4 for both C60 and C100 concrete samples that the P-CMOD curves obtained with the linear softening law are higher than with other softening laws, and the results using the exponential softening law are the minimum.In Figure 4a, all predicted curves are higher than the experimental results.The numerical results are consistent with the experimental data on the whole, except when using the linear softening law.The specific predicted peak loads are listed in Table 2.By comparing the experimental mean peak  The fracture mode would be changed with a different notch, and the corresponding mechanical response would also be varied.Wu [44] modeled the fracture mode transitions of notched concrete beams with different offsets under three-point bending based on the modified PDs.The author focused mainly on the fracture mode transitions of the concrete beams with a large offset and critical value when the transition happened.Since the PF-CZM does not need to consider complex crack tracking, it can be presumed that it can capture the mode transitions in concrete beams as the notch offset increases.Consider a set of quasi-static tests carried out by John [53]; these tests contained six different notch offsets of the specimens.The geometry and loading boundary conditions of the specimen are shown in Figure 5. Since the material parameters are not provided in the literature and the material parameters taken in different studies vary, the material parameters we adopted here are: Young's modulus E 0 = 31.5GPa, Poisson's ratio v 0 = 0.17, failure strength f t = 3.19 MPa, fracture energy G c = 99.67N/m.The maximum tensile stress failure criterion and the Petersson [49] softening curve are adopted with the method proposed in Section 3. To improve computational efficiency, we used a subdomain implementation scheme.The mesh size of the element used to discretize the damage subdomain is taken as h = 0 /5, while a larger mesh size is used for the elastic domain, as shown in Figure 6.In addition, the length scale parameter 0 is 2.5 mm, the plane stress assumption is considered, and the CPS4T and CPS4R elements are adopted for the damage subdomain and the elastic domain, respectively, in Abaqus.Displacement-controlled is used for the acquisition of response curves, and the prescribed displacement u = 0.6 mm is applied in 60 incremental steps, with each incremental step ∆u fixed at 0.01 mm.CZM does not need to consider complex crack tracking, it can be presumed that it can capture the mode transitions in concrete beams as the notch offset increases.Consider a set of quasi-static tests carried out by John [53]; these tests contained six different notch offsets of the specimens.The geometry and loading boundary conditions of the specimen are shown in Figure 5. Since the material parameters are not provided in the literature and the material parameters taken in different studies vary, the material parameters we adopted here are: Young's modulus  0 = 31.5GPa, Poisson's ratio  0 = 0.17, failure strength   = 3.19 MPa, fracture energy   = 99.67N/m.The maximum tensile stress failure criterion and the Petersson [49] softening curve are adopted with the method proposed in Section 3. To improve computational efficiency, we used a subdomain implementation scheme.The mesh size of the element used to discretize the damage subdomain is taken as ℎ = ℓ 0 5, while a larger mesh size is used for the elastic domain, as shown in Figure 6.In addition, the length scale parameter ℓ 0 is 2.5 mm, the plane stress assumption is considered, and the CPS4T and CPS4R elements are adopted for the damage subdomain and the elastic domain, respectively, in Abaqus.Displacement-controlled is used for the acquisition of response curves, and the prescribed displacement  = 0.6 mm is applied in 60 incremental steps, with each incremental step Δ fixed at 0.01 mm.  Figure 7 gives a comparison of the crack growth with the experimental results using the PF-CZM and the peridynamic model for six different notch offset specimens, with γ = 0, 0.5, 0.7, 0.72, 0.77 and 0.875 from top to bottom.The two groups of specimens with γ = 0.5 and 0.7 in Figure 7a are directly calculated using the span length L. As can be seen from Figure 7, PF-CZM (Figure 7a) and IH-PD models (Figure 7c) are in good agreement with the experimental results, including the crack growth angle and fracture mode, except for a slight difference with the experimental results at γ = 0.7.However, the IH-PD model has taken into account the factor of mesoscale inhomogeneity and the predicted results of the FH-PD model (Figure 7d) without this factor differing greatly from the experimental results; the crack initiates and grows from the notch tip, no matter where the notch is located.For the PF-CZM, without considering the mesoscale inhomogeneity, it is still possible to capture the twice transitions of the fracture mode of the specimen, and the damage zone is clearer and more explicit compared to the peridynamic modeling.However, it should be noted that PF-CZM still has the problem of computational efficiency.When the notch is located in the middle of the span (γ = 0), the specimen has a typical mode-I failure.Within a certain offset range, the crack growths will gradually change from mode-I to a mixed mode fracture as the notch offset increases.When the offset is greater than the critical value of offset distance (γ = 0.7), the fracture will occur anywhere in the mid-span of the beam instead of extending along the cut tip; this feature can be well reproduced by the PF-CZM.It can also be found that when the offset is 50 mm (γ = 0.5), the results obtained using the PF-CZM exhibit a tendency to change to mode-I in the late stage of the fracture.
Figure 6.In addition, the length scale parameter ℓ 0 is 2.5 mm, the plane stress assumption is considered, and the CPS4T and CPS4R elements are adopted for the damage subdomain and the elastic domain, respectively, in Abaqus.Displacement-controlled is used for the acquisition of response curves, and the prescribed displacement  = 0.6 mm is applied in 60 incremental steps, with each incremental step Δ fixed at 0.01 mm.  Figure 7 gives a comparison of the crack growth with the experimental results using the PF-CZM and the peridynamic model for six different notch offset specimens, with γ = 0, 0.5, 0.7, 0.72, 0.77 and 0.875 from top to bottom.The two groups of specimens with γ = 0.5 and 0.7 in Figure 7a are directly calculated using the span length L. As can be seen from Figure 7, the PF-CZM (Figure 7a) and IH-PD models (Figure 7c) are in good agreement with the experimental results, including the crack growth angle and fracture mode, except for a slight difference with the experimental results at γ = 0.7.However, the IH-PD model has taken into account the factor of mesoscale inhomogeneity and the predicted results of the FH-PD model (Figure 7d) without this factor differing greatly from the experimental results; the crack initiates and grows from the notch tip, no matter where the notch is located.For the PF-CZM, without considering the mesoscale inhomogeneity, it is still possible to capture the twice transitions of the fracture mode of the specimen, and the damage zone is clearer and more explicit compared to the peridynamic modeling.However, it should be noted that PF-CZM still has the problem of computational efficiency.
When the notch is located in the middle of the span γ = 0), the specimen has a typical mode-I failure.Within a certain offset range, the crack growths will gradually change from mode-I to a mixed mode fracture as the notch offset increases.When the offset is greater than the critical value of offset distance γ = 0.7), the fracture will occur anywhere in the mid-span of the beam instead of extending along the cut tip; this feature can be well reproduced by the PF-CZM.It can also be found that when the offset is 50 mm γ = 0.5 , the results obtained using the PF-CZM exhibit a tendency to change to mode-I in the late stage of the fracture.In the experimental results, a critical state for the second transition of the fracture mode exists at γ = 0.7.The PF-CZM can naturally capture this process of fracture mode transition states.Figure 8 gives the damage profile of the three groups of specimens at γ = 0.70, 0.72 and 0.77, when the crack initiates in the span for the first time.It can be seen that when γ gradually approaches 0.7, the damage zone at the tip of the notch gradually expands.Compared with γ = 0.50 in Figure 7a, the damage continues to extend directly along the notch tip and does not initiate in the span, which can prove that, at this time, there is indeed a critical state.In the experimental results, a critical state for the second transition of the fracture mode exists at γ = 0.7.The PF-CZM can naturally capture this process of fracture mode transition states.Figure 8 gives the damage profile of the three groups of specimens at γ = 0.70, 0.72 and 0.77, when the crack initiates in the span for the first time.It can be seen that when γ gradually approaches 0.7, the damage zone at the tip of the notch gradually expands.Compared with γ = 0.50 in Figure 7a, the damage continues to extend directly along the notch tip and does not initiate in the span, which can prove that, at this time, there is indeed a critical state.To further validate the model, we quantitatively compared the peak loads of concrete beams at different offsets using different methods, and the results are shown in Figure 9.The peak loads obtained from the experimental and theoretical analyses in [53], the PF-CZM, and the peridynamic model [44] are shown in Figure 9.It can be noted that the peak loads of the PF-CZM do not differ much from the experimental results as well as the analytical solution, and the evolution tendency with γ is consistent with the test, i.e., the effect of the notch offset on the limit peak capacity of the specimen gradually decreases with an increase in the notch offset, which effectively verifies the accuracy of the model.To further validate the model, we quantitatively compared the peak loads of concrete beams at different offsets using different methods, and the results are shown in Figure 9.The peak loads obtained from the experimental and theoretical analyses in [53], the PF-CZM, and the peridynamic model [44] are shown in Figure 9.It can be noted that the peak loads of the PF-CZM do not differ much from the experimental results as well as the analytical solution, and the evolution tendency with γ is consistent with the test, i.e., the effect of the notch offset on the limit peak capacity of the specimen gradually decreases with an increase in the notch offset, which effectively verifies the accuracy of the model.To further validate the model, we quantitatively compared the peak loads of concrete beams at different offsets using different methods, and the results are shown in Figure 9.The peak loads obtained from the experimental and theoretical analyses in [53], the PF-CZM, and the peridynamic model [44] are shown in Figure 9.It can be noted that the peak loads of the PF-CZM do not differ much from the experimental results as well as the analytical solution, and the evolution tendency with γ is consistent with the test, i.e., the effect of the notch offset on the limit peak capacity of the specimen gradually decreases with an increase in the notch offset, which effectively verifies the accuracy of the model.Where δ = |P 1 −P 2 | P 2 × 100%, P 1 is the predicted peak load using the PF-CZM and IH-PD, and P 2 is the analytical results from [53].It can be found in Figure 10 that the relative error of the PF-CZM is smaller than IH-PD for the same group materials with different notch features.Thus, we could assume that once the optimal softening law is determined, it is suitable for the same group of materials.However, it should be noted that whether the method is applicable to large-scale structures depends on further research.[44], the experimental results and the analytical solution [53].
Where  = | 1 − 2 |  2 × 100%,  1 is the predicted peak load using the PF-CZM and IH-PD, and  2 is the analytical results from [53].It can be found in Figure 10 that the relative error of the PF-CZM is smaller than IH-PD for the same group materials with different notch features.Thus, we could assume that once the optimal softening law is determined, it is suitable for the same group of materials.However, it should be noted that whether the method is applicable to large-scale structures depends on further research.

Analysis of the Mechanical Behaviors of Double-Notched Concrete Beams under Three-Point Bending
In actual concrete structures, there is usually more than one crack, and its distribution is often stochastic.It is costly to conduct systematic experimental studies.Even for simple double-notched concrete beams, the location of the notch is diverse.In this section, the mechanical response and fracture features of a typical double-notched concrete beam under three-point bending are analyzed, and the results are compared qualitatively with the experimental results.
The geometry and boundary conditions are shown in Figure 11, where the primary notch depth  0 is located in the middle of the concrete beam, the secondary notch depth is noted as  1 , and the distance between the primary and secondary notches is noted as .The following material parameters are considered:  0 = 33.5 GPa,  0 = 0.15,   = 2.96 MPa, with the linear softening curve considered.A prescribed displacement,  = 1 mm, is applied in 100 incremental steps, with each incremental step Δ fixed at 0.01 mm.For convenience, the primary notch depth  0 is fixed at 60 mm with the corresponding fracture energy   = 64 N/m.Only two groups of notch distance  = 30 and 60 mm are considered in this work; for each notch distance, three different secondary notch depths are selected (i.e.,  1 = 30, 45 and 60 mm).The secondary notch is not higher than the primary notch.

Analysis of the Mechanical Behaviors of Double-Notched Concrete Beams under Three-Point Bending
In actual concrete structures, there is usually more than one crack, and its distribution is often stochastic.It is costly to conduct systematic experimental studies.Even for simple double-notched concrete beams, the location of the notch is diverse.In this section, the mechanical response and fracture features of a typical double-notched concrete beam under three-point bending are analyzed, and the results are compared qualitatively with the experimental results.
The geometry and boundary conditions are shown in Figure 11, where the primary notch depth a 0 is located in the middle of the concrete beam, the secondary notch depth is noted as a 1 , and the distance between the primary and secondary notches is noted as d.The following material parameters are considered: E 0 = 33.5 GPa, v 0 = 0.15, f t = 2.96 MPa, with the linear softening curve considered.A prescribed displacement, u = 1 mm, is applied in 100 incremental steps, with each incremental step ∆u fixed at 0.01 mm.For convenience, the primary notch depth a 0 is fixed at 60 mm with the corresponding fracture energy G c = 64 N/m.Only two groups of notch distance d = 30 and 60 mm are considered in this work; for each notch distance, three different secondary notch depths are selected (i.e., a 1 = 30, 45 and 60 mm).The secondary notch is not higher than the primary notch.Figure 12 gives the numerical results for two typical specimens, and it is evident that the crack initiates and propagates at the primary notch of the specimens and that damage exists at the secondary notch as well.The closer the secondary notch is to the primary notch, the higher the notch depth and the more obvious it is.This is qualitatively in agreement with the experimental results of Zhang [54].In particular, the specimen exhibits a typical mode-I failure in Figure 12a; the effect of the secondary notch is negligible compared to Figure 12b.It can be clearly observed in Figure 12b that as the secondary notch depth increases and the notch distance decreases, the primary notch crack path exhibits a pronounced reverse perturbation by the secondary notch, even without considering the nonuniformity at the mesoscale.The phenomenon is consistent with the experimental observation in Figure 13, and it cannot be objectively explained using the CZM.It can be assumed that the crack path perturbation caused by local stress concentrations would be more obvious as the distance between the primary and secondary notches becomes   12 gives the numerical results for two typical specimens, and it is evident that the crack initiates and propagates at the primary notch of the specimens and that damage exists at the secondary notch as well.The closer the secondary notch is to the primary notch, the higher the notch depth and the more obvious it is.This is qualitatively in agreement with the experimental results of Zhang [54].In particular, the specimen exhibits a typical mode-I failure in Figure 12a; the effect of the secondary notch is negligible compared to Figure 12b.It can be clearly observed in Figure 12b that as the secondary notch depth increases and the notch distance decreases, the primary notch crack path exhibits a pronounced reverse perturbation by the secondary notch, even without considering the nonuniformity at the mesoscale.The phenomenon is consistent with the experimental observation in Figure 13, and it cannot be objectively explained using the CZM.It can be assumed that the crack path perturbation caused by local stress concentrations would be more obvious as the distance between the primary and secondary notches becomes smaller.However, the crack would turn back to mode-I failure (without considering the heterogeneity at the mesoscale explicitly), as shown in Figure 12b.Figure 12 gives the numerical results for two typical specimens, and it is evident that the crack initiates and propagates at the primary notch of the specimens and that damage exists at the secondary notch as well.The closer the secondary notch is to the primary notch, the higher the notch depth and the more obvious it is.This is qualitatively in agreement with the experimental results of Zhang [54].In particular, the specimen exhibits a typical mode-I failure in Figure 12a; the effect of the secondary notch is negligible compared to Figure 12b.It can be clearly observed in Figure 12b that as the secondary notch depth increases and the notch distance decreases, the primary notch crack path exhibits a pronounced reverse perturbation by the secondary notch, even without considering the nonuniformity at the mesoscale.The phenomenon is consistent with the experimental observation in Figure 13, and it cannot be objectively explained using the CZM.It can be assumed that the crack path perturbation caused by local stress concentrations would be more obvious as the distance between the primary and secondary notches becomes smaller.However, the crack would turn back to mode-I failure (without considering the heterogeneity at the mesoscale explicitly), as shown in Figure 12b.Furthermore, the predicted results of the peak load are shown in Table 3, and it can be found that when the secondary notch depth is the same, the limit load capacity of the beam is increased accordingly as the notch distance decreases.When the notch distance is fixed, the higher the secondary notch depth, the higher the limit load capacity of the beam.This phenomenon would be more obvious with an increase in the secondary notch depth.However, it should be noted that the increments of peak load are not obvious.On the other hand, this also shows the prospect of the applicability of the PF-CZM on the accurate failure modeling of multi-cracked structures.Furthermore, the predicted results of the peak load are shown in Table 3, and it can be found that when the secondary notch depth is the same, the limit load capacity of the beam is increased accordingly as the notch distance decreases.When the notch distance is fixed, the higher the secondary notch depth, the higher the limit load capacity of the beam.This phenomenon would be more obvious with an increase in the secondary notch depth.However, it should be noted that the increments of peak load are not obvious.On the other hand, this also shows the prospect of the applicability of the PF-CZM on the accurate failure modeling of multi-cracked structures.

Conclusions
In this study, we have proposed an approach combining optimal softening laws and the PF-CZM to model the fracture behaviors and mechanical properties of quasi-brittle materials.The main conclusions are: (1) A PF-CZM with Petersson's softening law has been realized and adopted to simulate the damage and fracture properties of quasi-brittle materials.
(2) The method has been validated by two tests of notched concrete beams under threepoint bending.The results indicate that the mechanical properties of concrete beams with the optimal softening law are better than the data reported using different methods.Further validation indicates that once the optimal softening law is determined, it is suitable for the same group of materials.(3) The modeling of concrete beams with different notch offsets indicates that the PF-CZM can naturally predict and reproduce the critical notch offset and the fracture transition process with the rising notch offset without considering mesoscale inhomogeneity.(4) The modeling of typical double-notched concrete beams indicates that the PF-CZM can predict the interaction between two notches objectively and the changes in the limit load capacity.

Figure 1 .
Figure 1.(a) A one-dimensional (1D) rod with a crack located at  = 0 and a cross-section .(b) Damage representation with a shape crack at  = 0. (c) Diffusive phase field approximation with WN (black), AT1 (blue) and AT2 (red) model scaling by the length scale parameter ℓ 0 .

Figure 1 .
Figure 1.(a) A one-dimensional (1D) rod with a crack located at x = 0 and a cross-section A. (b) Damage representation with a shape crack at x = 0. (c) Diffusive phase field approximation with WN (black), AT1 (blue) and AT2 (red) model scaling by the length scale parameter 0 .

Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a 2D solid  with its external boundary ∂.The external boundary ∂ can be divided into two parts: ∂  and ∂  , with ∂  ∩ ∂  = ∅ and ∂  ∪ ∂  = ∂, in which it is given displacements  * for ∂  and tractions  * for ∂  , respectively.The outward unit normal vector of external boundary ∂ is denoted as .

Figure 2 .
Figure 2. Schematic illustration of sharp (Γ, left) and diffusive (B, right) crack topology in a 2D solid Ω with its external boundary ∂Ω.The external boundary ∂Ω can be divided into two parts: ∂Ω t and ∂Ω u , with ∂Ω t ∩ ∂Ω u = ∅ and ∂Ω t ∪ ∂Ω u = ∂Ω, in which it is given displacements u * for ∂Ω u and tractions t * for ∂Ω t , respectively.The outward unit normal vector of external boundary ∂Ω is denoted as n.

Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a 2D solid  with its external boundary ∂.The external boundary ∂ can be divided into two parts: ∂ and ∂ , with ∂ ∩ ∂ = ∅ and ∂ ∪ ∂ = ∂, in which it is given displacements  * for ∂ and tractions  * for ∂ , respectively.The outward unit normal vector of external boundary ∂ is denoted as .
Appl.Sci.2022, 12, x FOR PEER REVIEW 2 o Appl.Sci.2022, 12, x FOR PEER REVIEW (, ∇) sharp and diffusive crack topology is shown in Figure 2. Th field model diffuses sharp crack  into crack band ℬ ⊆  with finite scale ℓ is the external boundary of the crack band ℬ, The outward unit normal vector o denoted as  ℬ .The body force is denoted as  * .The regularized surface energy obtained by Equation (8). =  d =  (, ∇)d =  ℬ (, ∇)d where  is the critical energy release rate.

Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a  with its external boundary ∂.The external boundary ∂ can be divided into two p and ∂ , with ∂ ∩ ∂ = ∅ and ∂ ∪ ∂ = ∂, in which it is given displacemen ∂ and tractions  * for ∂ , respectively.The outward unit normal vector of external b ∂ is denoted as .
sharp and diffusive crack topology is shown in Figure 2. The phase field model diffuses sharp crack  into crack band ℬ ⊆  with finite scale ℓ > 0. ∂ℬ is the external boundary of the crack band ℬ, The outward unit normal vector of ∂ℬ is denoted as  ℬ .The body force is denoted as  * .The regularized surface energy can be obtained by Equation (8). =  d =  (, ∇)d =  ℬ (, ∇)d (8)where  is the critical energy release rate.

Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a 2D solid  with its external boundary ∂.The external boundary ∂ can be divided into two parts: ∂ and ∂ , with ∂ ∩ ∂ = ∅ and ∂ ∪ ∂ = ∂, in which it is given displacements  * for ∂ and tractions  * for ∂ , respectively.The outward unit normal vector of external boundary ∂ is denoted as .

Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a 2D s  with its external boundary ∂.The external boundary ∂ can be divided into two parts: and ∂ , with ∂ ∩ ∂ = ∅ and ∂ ∪ ∂ = ∂, in which it is given displacements  * ∂ and tractions  * for ∂ , respectively.The outward unit normal vector of external bound ∂ is denoted as .

Figure 2 .Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a 2D so  with its external boundary ∂.The external boundary ∂ can be divided into two parts: ∂ and ∂ , with ∂ ∩ ∂ = ∅ and ∂ ∪ ∂ = ∂, in which it is given displacements  * ∂ and tractions  *for ∂ , respectively.The outward unit normal vector of external bound ∂ is denoted as .

Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a 2D solid  with its external boundary ∂.The external boundary ∂ can be divided into two parts: ∂ and ∂ , with ∂ ∩ ∂ = ∅ and ∂ ∪ ∂ = ∂, in which it is given displacements  * for ∂ and tractions  *for ∂ , respectively.The outward unit normal vector of external boundary ∂ is denoted as .
) =   ⋅ (), () = 1 +   +    (11) (φ)δφ i. 2022, 12, x FOR PEER REVIEW 2 of 19 ), imposing that δW int = δW ext should hold for the arbitrary values of δu and δφ.This leads to the strong form of the governing equations: Two-dimensional sharp and diffusive crack topology is shown in Figure 2. The ph field model diffuses sharp crack  into crack band ℬ ⊆  with finite scale ℓ > 0. is the external boundary of the crack band ℬ, The outward unit normal vector of ∂ℬ denoted as  ℬ .The body force is denoted as  * .The regularized surface energy can obtained by Equation (8). =  d =  (, ∇)d =  ℬ (, ∇)dwhere  is the critical energy release rate.

Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a 2D so  with its external boundary ∂.The external boundary ∂ can be divided into two parts: and ∂ , with ∂ ∩ ∂ = ∅ and ∂ ∪ ∂ = ∂, in which it is given displacements  * ∂ and tractions  *for ∂ , respectively.The outward unit normal vector of external bound ∂ is denoted as .
sharp and diffusive crack topology is shown in Figure 2. The ph field model diffuses sharp crack  into crack band ℬ ⊆  with finite scale ℓ > 0. is the external boundary of the crack band ℬ, The outward unit normal vector of ∂ denoted as  ℬ .The body force is denoted as  * .The regularized surface energy can obtained by Equation (8). =  d =  (, ∇)d =  ℬ (, ∇)dwhere  is the critical energy release rate.

Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a 2D s  with its external boundary ∂.The external boundary ∂ can be divided into two parts: and ∂ , with ∂ ∩ ∂ = ∅ and ∂ ∪ ∂ = ∂, in which it is given displacements  * ∂ and tractions  *for ∂ , respectively.The outward unit normal vector of external boun ∂ is denoted as .
sharp and diffusive crack topology is shown in field model diffuses sharp crack  into crack band ℬ ⊆  with fin is the external boundary of the crack band ℬ, The outward unit nor denoted as  ℬ .The body force is denoted as  * .The regularized su obtained by Equation (8). =  d =  (, ∇)d =  ℬ (, ∇)

Figure 2 .
Figure 2. Schematic illustration of sharp (, left) and diffusive (ℬ, right) crack topology in a 2D solid  with its external boundary ∂.The external boundary ∂ can be divided into two parts: ∂ and ∂ , with ∂ ∩ ∂ = ∅ and ∂ ∪ ∂ = ∂, in which it is given displacements  * for ∂ and tractions  *for ∂ , respectively.The outward unit normal vector of external boundary ∂ is denoted as .

Figure 3 .
Figure 3. Geometry of the notched concrete beam under three-point bending.

Figure 3 .
Figure 3. Geometry of the notched concrete beam under three-point bending.

Figure 4 .
Figure 4. P-CMOD curves obtained from exponential, linear and Petersson's softening laws using the PF-CZM, the concrete damage plasticity (CDP) model with Petersson's softening law and the experimental results [43] of (a) C60 concrete and (b) C100 concrete.

Figure 4 .
Figure 4. P-CMOD curves obtained from exponential, linear and Petersson's softening laws using the PF-CZM, the concrete damage plasticity (CDP) model with Petersson's softening law and the experimental results [43] of (a) C60 concrete and (b) C100 concrete.

19 Figure 8 .
Figure 8. Damage profile when crack initiates at the mid-span for the first time: (a γ = 0.70; (b γ = 0.72; (c γ = 0.77.(The phase field values are rounded to one decimal place). b

19 Figure 8 .
Figure 8. Damage profile when crack initiates at the mid-span for the first time: (a γ = 0.70; (b γ = 0.72; (c γ = 0.77.(The phase field values are rounded to one decimal place).

Figure 10 .
Figure 10.Relative error  of predicted peak loads with different γ values from PF-CZM and IH-PD.

Figure 10 .
Figure 10.Relative error δ of predicted peak loads with different γ values from PF-CZM and IH-PD.

Figure 11 .
Figure 11.A concrete specimen: geometry, loading and boundary conditions.

Figure
Figure 12  gives the numerical results for two typical specimens, and it is evident that the crack initiates and propagates at the primary notch of the specimens and that damage exists at the secondary notch as well.The closer the secondary notch is to the primary notch, the higher the notch depth and the more obvious it is.This is qualitatively in agreement with the experimental results of Zhang[54].In particular, the specimen exhibits a typical mode-I failure in Figure12a; the effect of the secondary notch is negligible compared to Figure12b.It can be clearly observed in Figure12bthat as the secondary notch depth increases and the notch distance decreases, the primary notch crack path exhibits a pronounced reverse perturbation by the secondary notch, even without considering the nonuniformity at the mesoscale.The phenomenon is consistent with the experimental observation in Figure13, and it cannot be objectively explained using the CZM.It can be assumed that the crack path perturbation caused by local stress concentrations would be more obvious as the distance between the primary and secondary notches becomes smaller.However, the crack would turn back to mode-I failure (without considering the heterogeneity at the mesoscale explicitly), as shown in Figure12b.

Figure 11 .
Figure 11.A concrete specimen: geometry, loading and boundary conditions.

Table 1 .
Material parameters for notched concrete beams.Geometry of the notched concrete beam under three-point bending.

Table 1 .
Material parameters for notched concrete beams.

Table 1 .
Material parameters for notched concrete beams.

Table 2 .
Predicted peak loads using the PF-CZM with different softening laws.