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

Abstract

In this paper, we propose an approach combining optimal softening laws and a phase-field 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 first 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.

1. 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]

$$\Psi \left(\mathit{u},\Gamma \right)={\Psi }^b+{\Psi }^{\Gamma }={{\displaystyle \int }}_{\Omega }\psi \left(\mathit{ϵ},\varphi \right)\mathrm{d}\Omega +G_c{{\displaystyle \int }}_{\Gamma }\mathrm{d}\Gamma$$

(1)

where the total potential energy functional $\Psi \left(\mathit{u},\Gamma \right)$ is assumed to include the surface energy ${\Psi }^{\Gamma }$ of the crack and the bulk energy ${\Psi }^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.

2. Phase-Field-Regularized Cohesive Zone Model

2.1. 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 $\varphi \left(x\right)\in \left[0, 1\right]$ can be introduced to describe the sharp crack topology; then, we can get Equation (2) based on the variational principle. The crack phase field $\varphi \left(x\right)$ located at $x$ = 0 can be derived by (3); see more details in [1,18].

$$A={{\displaystyle \int }}_{\Gamma }\mathrm{d}\Gamma ={{\displaystyle \int }}_{\Gamma }{{\displaystyle \int }}_{-\infty }^{\infty }\gamma \left(\varphi \right)\mathrm{d}x\mathrm{d}\Gamma ={{\displaystyle \int }}_{\Omega }\gamma \left(\varphi \right)\mathrm{d}\Omega =A\left(\varphi \right)$$

(2)

$$x\left(\varphi \right)=l_0{{\displaystyle \int }}_{\varphi }^1{\alpha }^{-\frac{1}{2}}(\widehat{\varphi })\mathrm{d}\widehat{\varphi }$$

(3)

where $A$ denotes the fracture surface area, and the crack surface density function is $\gamma \left(\varphi \right)≔\frac{1}{c_0}\left[{\ell }_0{\left({\varphi }^{\prime }\right)}^2+\frac{1}{{\ell }_0}\alpha \left(\varphi \right)\right]$; $\alpha \left(\varphi \right)=\xi \varphi +\left(1-\xi \right){\varphi }^2$ is denoted as the crack geometry function, $c_0=4{\int }_0^1\sqrt{\alpha \left(\varphi \right)}\mathrm{d}\varphi$ is the scaling parameter and ${\ell }_0$ is a length scale parameter.

As demonstrated by Wu [1,25], some conventional models can be obtained with different $\xi \in \left[0, 2\right]$ of the crack geometry function $\alpha \left(\varphi \right)$

• AT2 model [16,18]

$$\xi =0, \varphi \left(x\right)=\mathit{exp}\left(-\frac{\left|x\right|}{{\ell }_0}\right)$$

(4)

• AT1 model [45]

$$\xi =1, \varphi \left(x\right)={\left(1-\frac{\left|x\right|}{2{\ell }_0}\right)}^2$$

(5)

• WN model [1,46]

$$\xi =2, \varphi \left(x\right)=1-\mathit{sin}\left(\frac{\left|x\right|}{{\ell }_0}\right)$$

(6)

In this paper, we intend to use the WN model, i.e., $\xi =2$; thus, we have $\alpha \left(\varphi \right)=2\varphi -{\varphi }^2$, $c_0=\pi$, $\gamma \left(\varphi \right)=\frac{1}{\pi }\left[{\ell }_0{\left({\varphi }^{\prime }\right)}^2+\frac{1}{{\ell }_0}\left(2\varphi -{\varphi }^2\right)\right]$. Extending the 1D crack surface density function $\gamma \left(\varphi \right)$ into a multidimension gives Equation (7).

$$\gamma \left(\varphi ,\nabla \varphi \right)=\frac{1}{c_0}\left[{\ell }_0|\nabla \varphi {|}^2+\frac{1}{{\ell }_0}\alpha \left(\varphi \right)\right]=\frac{1}{\pi }\left({\ell }_0|\nabla \varphi {|}^2+\frac{1}{{\ell }_0}\left(2\varphi -{\varphi }^2\right)\right)$$

(7)

Two-dimensional sharp and diffusive crack topology is shown in Figure 2. The phase field model diffuses sharp crack $\Gamma$ into crack band $ℬ\subseteq \Omega$ with finite scale ${\ell }_0>0$. $\partial ℬ$ is the external boundary of the crack band $ℬ$, The outward unit normal vector of $\partial ℬ$ is denoted as ${\mathit{n}}_{ℬ}$. The body force is denoted as ${\mathit{b}}^{*}$. The regularized surface energy can be obtained by Equation (8).

$${\Psi }^{\Gamma }=G_c{{\displaystyle \int }}_{\Gamma }\mathrm{d}\Gamma =G_c{{\displaystyle \int }}_{\Omega }\gamma \left(\varphi ,\nabla \varphi \right)\mathrm{d}\Omega =G_c{{\displaystyle \int }}_{ℬ}\gamma \left(\varphi ,\nabla \varphi \right)\mathrm{d}\Omega$$

(8)

where $G_c$ is the critical energy release rate.

Correspondingly, the bulk energy is

$${\Psi }^b={{\displaystyle \int }}_{\Omega }\psi \left(\mathit{ϵ},\varphi \right)\mathrm{d}\Omega$$

(9)

in which $\psi \left(\mathit{ϵ},\varphi \right)=\omega \left(\varphi \right){\psi }_0\left(\mathit{ϵ}\right)$ is the bulk energy density, ${\psi }_0\left(\mathit{ϵ}\right)=\frac{1}{2}\mathit{ϵ}:{\mathbb{E}}_0:\mathit{ϵ}=\frac{1}{2}\mathit{\sigma }:{ℂ}_0:\mathit{\sigma }$ is the elastic strain energy, $\mathit{\sigma }$ and $\mathit{ϵ}$ are the stress tensor and the strain tensor, respectively. ${\mathbb{E}}_0$ and ${ℂ}_0$ are the elastic stiffness tensor and the flexibility tensor, respectively. On the energy degradation function $\omega \left(\varphi \right)\in \left[0, 1\right]$, 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.

$$\omega \left(\varphi \right)≔\frac{1}{1+g\left(\varphi \right)}=\frac{{(1-\varphi )}^p}{{(1-\varphi )}^p+Q\left(\varphi \right)},g\left(\varphi \right)=\frac{Q\left(\varphi \right)}{{(1-\varphi )}^p}$$

(10)

$$Q\left(\varphi \right)=a_1\varphi \cdot P\left(\varphi \right),P\left(\varphi \right)=1+a_2\varphi +a_2{a}_3{\varphi }^2$$

(11)

2.2. Principle of Virtual Work and Governing Equations

By using the results from the previous section, the total potential energy functional can be written as

$$\Psi \left(\mathit{u},\varphi \right)={{\displaystyle \int }}_{\Omega }\omega \left(\varphi \right){\psi }_0\left(\mathit{ϵ}\left(\mathit{u}\right)\right)\mathrm{d}\Omega +{{\displaystyle \int }}_{\Omega }{G}_c\gamma \left(\varphi ,\nabla \varphi \right)\mathrm{d}\Omega$$

(12)

Here, the compact tensor is not used. For a linear-elastic material, ${\psi }_0\left(\mathit{ϵ}\right)=\frac{1}{2}\lambda {\epsilon }_{kk}{\epsilon }_{ll}+\mu {\epsilon }_{ij}{\epsilon }_{ij}$, with $\lambda$ and $\mu$ being the Láme constants; $i,j,k,l=1,\dots ,d,$ with $d$ being dimensional numbers. For small strains, i.e., ${\epsilon }_{ij}=\frac{1}{2}\left(\partial u_i/\partial x_j+\partial u_j/\partial x_i\right)$. The governing equations are obtained as follows: we first define the variation of the external work increment as

$$\delta W_{\mathrm{ext} }={{\displaystyle \int }}_{\Omega }{b}_j^{*}\delta u_j\mathrm{d}\Omega +{{\displaystyle \int }}_{\partial {\Omega }_t}{t}_j^{*}\delta u_j\mathrm{d}A$$

(13)

On the other hand, the variation of the internal energy increment is given by

$$\delta W_{\mathrm{int} }=\delta \Psi =\frac{\partial \Psi }{\partial {\epsilon }_{ij}}\delta {\epsilon }_{ij}+\frac{\partial \Psi }{\partial \varphi }\delta \varphi$$

(14)

which for the case of (12) yields

$$\begin{gathered} \delta \Psi ={{\displaystyle \int }}_{\Omega }\omega \left(\varphi \right){\sigma }_{ij}\delta {\epsilon }_{ij}\mathrm{d}\Omega +{{\displaystyle \int }}_{\Omega }{\omega }^{\prime }\left(\varphi \right)\delta \varphi {\psi }_0\left(\mathit{ϵ}\right)\mathrm{d}\Omega \\ +{{\displaystyle \int }}_{\Omega }{G}_c\frac{1}{c_0}\left(2{\ell }_0\frac{\partial \varphi }{\partial x_i}\frac{\partial \delta \varphi }{\partial x_i}+\frac{1}{{\ell }_0}{\alpha }^{\prime }\left(\varphi \right)\delta \varphi \right)\mathrm{d}\Omega \end{gathered}$$

(15)

By appropriate constant transformation with divergence theorem, it can be found that the above expression is equivalent to

$$\begin{gathered} \delta \Psi ={{\displaystyle \int }}_{\partial \Omega }\omega \left(\varphi \right)n_i{\sigma }_{ij}\delta u_j\mathrm{d}A-{{\displaystyle \int }}_{\Omega }\omega \left(\varphi \right)\frac{\partial {\sigma }_{ij}}{\partial x_i}\delta u_j\mathrm{d}\Omega \\ +{{\displaystyle \int }}_{\Omega }{\omega }^{\prime }\left(\varphi \right)\delta \varphi {\psi }_0\left(\mathit{ϵ}\right)\mathrm{d}\Omega +{{\displaystyle \int }}_{\partial \Omega }{G}_c\frac{1}{c_0}2{\ell }_0\frac{\partial \varphi }{\partial x_i}{n}_i\delta \varphi \mathrm{d}A \\ +{{\displaystyle \int }}_{\Omega }{G}_c\frac{1}{c_0}\left(-2{\ell }_0\frac{{\partial }^2\varphi }{\partial x_i\partial x_i}\delta \varphi +\frac{1}{{\ell }_0}{\alpha }^{\prime }\left(\varphi \right)\delta \varphi \right)\mathrm{d}\Omega \end{gathered}$$

(16)

Combine the terms in (13) and (16), imposing that $\delta W_{int}=\delta W_{ext}$ should hold for the arbitrary values of $\delta u$ and $\delta \varphi$. This leads to the strong form of the governing equations:

$$\omega \left(\varphi \right)\frac{\partial {\sigma }_{ij}}{\partial x_i}+b_j^{*}=0 \mathrm{in} \Omega$$

(17)

$$\omega \left(\varphi \right){\sigma }_{ij}{n}_i=t_j^{*} \mathrm{on} \partial {\Omega }_t$$

(18)

$$u_j=u_j^{*} \mathrm{on} \partial {\Omega }_u$$

(19)

$$-\frac{2{\ell }_0}{c_0}{G}_c\frac{{\partial }^2\varphi }{\partial x_i\partial x_i}+\frac{G_c}{c_0{\ell }_0}{\alpha }^{\prime }\left(\varphi \right)+{\omega }^{\prime }\left(\varphi \right){\psi }_0\left(ϵ\right)=0 \mathrm{in} \Omega$$

(20)

$$\frac{\partial \varphi }{\partial x_i}{n}_i=0 \mathrm{on} \partial \Omega$$

(21)

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 ${\psi }_0\left(\mathit{ϵ}\right)$, in order to consider the different mechanical behaviors under tensile and compressive stress.

$$\overline{Y}=\frac{1}{2E_0}\langle {\sigma }_1\rangle {}^2$$

(22)

where $E_0$ is Young’s modulus; the Macaulay bracket 〈$\cdot$〉 is defined as $\langle x\rangle =max\left(x,0\right)$, and ${\sigma }_1$ is the maximum principal stress of the stress tensor $\mathit{\sigma }$.

2.3. Damage Irreversibility

To deal with the irreversibility of phase field evolution $(\dot{\varphi }\ge 0)$, we adopt the history field $H$ to prevent crack healing, which can ensure

$${\varphi }_{t+\mathsf{\Delta }t}\ge {\varphi }_t$$

(23)

where ${\varphi }_{t+\mathsf{\Delta }t}$ represents the phase field value in the current time increment while ${\varphi }_t$ denotes the value of the previous increment. The history field variable must satisfy the following Kuhn–Tucker conditions [48]

$$\overline{Y}-H\le 0,\dot{H}\ge 0,\dot{H}\left(\overline{Y}-H\right)=0$$

(24)

Therefore, the history field variable at the current moment can be taken as

$$H=\underset{\tau \in \left[0,t\right]}{max}\left(\overline{Y_0},\overline{Y}\left(\tau \right)\right),\overline{Y_0}\frac{1}{2E_0}{f}_t^2$$

(25)

Then, the evolution equation of the phase field is

$$-\frac{2{\ell }_0}{c_0}{G}_c\frac{{\partial }^2\varphi }{\partial x_i\partial x_i}+\frac{G_c}{c_0{\ell }_0}{\alpha }^{\prime }\left(\varphi \right)+{\omega }^{\prime }\left(\varphi \right)H=0 \mathrm{in} \Omega$$

(26)

2.4. Phase-Field-Regularized Cohesive Zone Model

The geometric crack function $\alpha \left(\varphi \right)$ and the energetic degradation function $\omega \left(\varphi \right)$ 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\left(\varphi \right)$ in the geometric crack function $\alpha \left(\varphi \right)$.

$$\{\begin{array}{l}\alpha \left(\varphi \right)=\xi \varphi +\left(1-\xi \right){\varphi }^2 \xi \in \left[0,2\right]\\ \omega \left(\varphi \right)=\frac{1}{1+g\left(\varphi \right)} g\left(\varphi \right)=\frac{a_1\varphi \left(1+a_2\varphi +a_2{a}_3{\varphi }^2\right)}{{(1-\varphi )}^p}\end{array}$$

(27)

in which the exponent $p\ge 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 $\sigma \left(w\right)$, respectively.

$$k_0=-\frac{c_0}{4\pi }\frac{f_t^2}{G_c}\frac{{\left[\xi \left(a_2+p+1\right)-1\right]}^{3/2}}{{\xi }^2}$$

(28)

$$w_{\mathrm{c}}=\frac{2\pi G_c}{c_0{f}_t}\sqrt{\xi P\left(1\right)}\underset{{\varphi }^{*}\to 1}{lim}{\left(1-{\varphi }^{*}\right)}^{1-p/2}\{\begin{array}{l}0 p<2\\ \frac{2\pi G_c}{c_0{f}_t}\sqrt{\xi P\left(1\right)} p=2\\ +\infty p>2\end{array}$$

(29)

where $P\left(1\right)=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 $\omega \left(\varphi \right)$ is taken in the range of $p\ge 2$. For the PF-CZM, $\xi =2$ is taken; then, we have $\alpha \left(\varphi \right)=2\varphi -{\varphi }^2$ and $c_0=\pi$.

$$a_1=\frac{4}{\pi l_0}{l}_{\mathrm{ch}}, a_2=2{\left(-2k_0\frac{G_c}{f_t^2}\right)}^{2/3}-\left(p+\frac{1}{2}\right)$$

(30)

$$a_3=\{\begin{array}{l}0 p>2\\ \frac{1}{a_2}\left[\frac{1}{8}{\left(\frac{w_c{f}_t}{G_c}\right)}^2-\left(1+a_2\right)\right] p=2\end{array}$$

(31)

for Griffith’s characteristic length $l_{\mathrm{ch}}:=E_0{G}_c/f_t^2$. 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:

• Linear softening curve: $p=2, a_2=-0.5, a_3=0$
• Exponential softening curve: $p=2.5, a_2=2^{5/3}-3, a_3=0$
• Petersson [49] softening curve: $p=2, a_2=0.3114, a_3=0.9910$

2.5. 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].

3. 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 fiber-reinforced 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. Material parameters for notched concrete beams.

Strength Level	Young’s Modulus ${\mathit{E}}_0$ (GPa)	Poisson’s Ratio ${\mathit{v}}_0$	Failure Strength ${\mathit{f}}_{\mathit{t}}$ (MPa)	Fracture Energy ${\mathit{G}}_{\mathit{c}}$ (N/m)
C60	35.7	0.17	4.43	128.7
C100	41.4	0.17	5.71	147.5

. 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 ${\ell }_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 $\mathsf{\Delta }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. Predicted peak loads using the PF-CZM with different softening laws.

Strength Level	Experimental Mean Peak Load (kN)	Max Predicted Peak Load (kN)	Min Predicted Peak Load (kN)	Optimal Predicted Peak Load (kN)	Optimal Softening Law
C60	3.501	4.055	3.618	3.618	Exponential
C100	4.761	4.967	4.418	4.553	Petersson

. 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.

4. Numerical Examples and Discussion

4.1. Analysis of the Fracture Mode and Mechanical Response of Concrete Beams under Three-Point Bending with Different Notch Offsets

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.5 GPa, Poisson’s ratio $v_0$ = 0.17, failure strength $f_t$ = 3.19 MPa, fracture energy $G_c$ = 99.67 N/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={\ell }_0/5$, while a larger mesh size is used for the elastic domain, as shown in Figure 6. In addition, the length scale parameter ${\ell }_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 $\mathsf{\Delta }u$ 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.

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 $\delta =\frac{\left|P_1-P_2\right|}{P_2}\times 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.

4.2. 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 $\mathsf{\Delta }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 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. Predicted peak loads of double-notched concrete beams with different $d$ (30 mm, 60 mm) and $a_1$ (30 mm, 45 mm, 60 mm) values.

Notch Depth ${\mathit{a}}_1$, $\mathit{d}$ = 30 (mm)	Peak Load (kN)	Notch Depth ${\mathit{a}}_1$, $\mathit{d}$ = 60 (mm)	Peak Load (kN)
30	4.062	30	4.061
45	4.069	45	4.065
60	4.205	60	4.133

, 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.

5. 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 three-point 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.