A Particle-Based Cohesive Crack Model for Brittle Fracture Problems

Numerical simulations of the fracture process are challenging, and the discrete element (DE) method is an effective means to model fracture problems. The DE model comprises the DE connective model and DE contact model, where the former is used for the representation of isotropic solids before cracks initiate, while the latter is employed to represent particulate materials after cracks propagate. In this paper, a DE particle-based cohesive crack model is developed to model the mixed-mode fracture process of brittle materials, aiming to simulate the material transition from a solid phase to a particulate phase. Because of the particle characteristics of the DE connective model, the cohesive crack model is constructed at inter-particle bonds in the connective stage of the model at a microscale. A potential formulation is adopted by the cohesive zone method, and a linear softening relation is employed by the traction–separation law upon fracture initiation. This particle-based cohesive crack model bridges the microscopic gap between the connective model and the contact model and, thus, is suitable to describe the material separation process from solids to particulates. The proposed model is validated by a number of standard fracture tests, and numerical results are found to be in good agreement with the analytical solutions. A notched concrete beam subjected to an impact loading is modeled, and the impact force obtained from the numerical modeling agrees better with the experimental result than that obtained from the finite element method.


Introduction
Fracture is a common failure mode for engineering structures and structural components. When structures are subjected to severe loading or large deformation, new fracture surfaces are created and cracks occur. Effective numerical modeling of the fracture process is important for assessing the safety, reliability, and structural performance of engineering structures. For modeling the fracture process, a number of numerical methods are often employed, such as the finite element method (FEM) [1], the extended finite element method (XFEM) [2], meshless methods [3], the particle finite element method (PFEM) [4], molecular dynamics [5], the particle method [6], and atomistic methods [7,8]. Among them, the atomistic method is able to provide great insight into the nanoscopic mechanism of fracture initiation and propagation since fracture problems essentially take place at the atomic level of materials by means of the breakage of bonds [9,10]. Modeling at the atomistic scale, however, is computationally intensive because extremely small length and time scales have to be adopted [11]. Therefore, atomistic modeling is generally not applicable to practical engineering problems. CZM, and a linear softening relation is employed for the traction-separation law upon fracture initiation. The criteria of both the fracture initiation and the fracture propagation are constructed in the displacement space. A damage parameter is used to record the damage state and, thus, the fracture process is irreversible. Note that some CZMs are also used in conjunction with the DEM for fracture modeling in homogeneous materials [18,[34][35][36] and laminated materials [30,37,38], but the cohesive formulation differs.
The remainder of the paper is organized as follows: the DE connective model is firstly outlined to describe the intact stage of elastic solids in Section 2. Then, a particle-based cohesive crack model is proposed in Section 3 to detail the fracture process of materials, where mixed-mode fracture initiation and propagation criteria are derived. Once the decohesion is completed, the conventional DE contact model is applied to describe the interaction of DE particles in Section 4. The model transition process and the implementation issue of the model are discussed in Sections 5 and 6, respectively. To validate the proposed model, a number of numerical examples are presented in Section 7. Lastly, some conclusions are drawn.

Connective Model: Representation of Isotropic Elastic Solid
The DE connective method is favorable for modeling elastic-brittle materials, such as glass [13], because this model can be conveniently switched to the contact model. As the DE connective model is capable of modeling the response of materials which are initially continuous but eventually cracked [19], this model is employed herein to represent the elastic continuum used in the domain of particular interest. It is worth noting that, in most DE connective models, all DE particles have the same geometrical shape and size for easy characterization of material properties.
In the DE connective model shown in Figure 1, spherical particles and virtual springs are used at microscale to formulate homogeneous and isotropic solids, generally described by the FEM using a phenomenological material constitutive law at macroscale. The formulation is based on the energy equivalence, i.e., the energy stored in springs is equal to that stored in elastic solids. The spring stiffness is derived based on the Cauchy-Born rule [39], and the details of this derivation are given in Reference [40].  A number of DE connective models were reported, for which the unit cell is either in a hexagonal structure [18,19,41] or in a cubic structure [20,[41][42][43]. In a hexagonal structure, the DE particles are stacked in a denser manner; however, the boundary of the model is not flat and, thus, not favorable for the coupling approach with the FEM. The cubic model proposed by Yu [20] was employed because of its computational accuracy and unique structure. This particular structure can result in flat boundaries with the DE model, which is practically very desirable for domain decomposition of coupling models because the burden of pre-processing can be effectively eased.
The unit cell structure of this model is shown in Figure 1b. It is apparent that this model is composed of 27 spherical particles of the same size which are located in a cubic structure. The central particle of this model is connected to 26 neighboring particles, which can be categorized into three groups according to their distance to the central one, as depicted by different numbers in Figure 1c.
The interaction forcef int between any two adjacent particles in the local coordinate system is calculated based on their relative displacement δ and spring stiffness K (see Figure 1d) as follows: where forcef int = f nfs1fs2 T and δ = [δ n , δ s1 , δ s2 ] T , while K is given as follows: Inside each pair, as shown in Figure 1d, there are one orthogonal (k n ) and two tangential (k s1 and k s2 ) linear springs invisibly connecting them. Note that subscripts n, s1, and s2 denote the unit directions of the local coordinate, as shown in Figure 1d. Their stiffness is determined based on the energy equivalence between that stored in the springs and that stored in solid elasticity, as given by Equation (3) [20].      k 1 n = k 2 n = 2Er 5(1−2ν) k 1 s1 = k 1 s2 = k 2 s1 = k 2 s2 = 2Er(1−4ν) 5(1−2ν)(1+ν) where E, ν, and r are the Young's modulus, Poisson's ratio, and the radius of DE particles, respectively. Note that superscripts in Equation (3) denote the type of connecting springs between each particle pair, as shown in Figure 1c. Let n, s 1 , and s 2 be the unit base vectors of the local coordinate system to be expressed in the global coordinate system, which is shown in Figure 1d. Then, the transformation matrix φ from the global frame to the local frame can be expressed by where n 1 , n 2 , and n 3 are the components of n, while s 1 1 , s 1 2 , and s 1 3 are the components of s 1 ; note that, here, the subscript replaced by superscript for clarity. The same interpretation and treatment of s 1 apply to s 2 .
Therefore, the internal force f int in the global coordinate system can be obtained by and the moment m int is calculated as where r denotes the effective radius vector, originating from the particle's center to the middle point between the two particles of a pair. Taking Figure 1d as an example, the effective radius vector for the particle i is r i = x c − x i , where x c and x i denote position vectors of the particle i and middle point c, respectively. Based on Equation (1), the relation between the traction T = [T n , T s1 , T s2 ] T and the relative displacement δ can be written as T n = k n δ n /A, T s1 = k s1 δ s1 /A, where T n , T s1 , and T s2 are the traction components along the normal and the two shear directions in the local frame, respectively, δ n , δ s1 , and δ s2 are the relative displacement components (or separation components) along the normal and the two shear directions, respectively, and A is an effective area between a particle pair; for the first and second nearest particle pairs, it is πr 2 /4 and 2πr 2 /9, respectively [40]. Note that, as the spring stiffness of the third type of connecting particle pair is zero as shown in Equation (3), the effective area of this type is not useful and, thus, ignored. Since spring stiffness in two shear directions are the same, as seen in Equation (3), i.e., k s = k s1 = k s2 , the two relative displacement components δ s1 and δ s2 can be considered together as Then, Equations (8) and (9) can be combined as T s = k s δ s /A. (11) The linear relation between the traction and relative displacement are illustrated in Figure 2.
Materials 2020, 13, x FOR PEER REVIEW 5 of 35 for the particle is ̅ , where and denote position vectors of the particle i and middle point , respectively.
Based on Equation (1), the relation between the traction T , T , T and the relative displacement can be written as / , where T , T , and T are the traction components along the normal and the two shear directions in the local frame, respectively, δ , δ , and δ are the relative displacement components (or separation components) along the normal and the two shear directions, respectively, and is an effective area between a particle pair; for the first and second nearest particle pairs, it is π /4 and 2π /9, respectively [40]. Note that, as the spring stiffness of the third type of connecting particle pair is zero as shown in Equation (3), the effective area of this type is not useful and, thus, ignored.
Since spring stiffness in two shear directions are the same, as seen in Equation (3), i.e., k k k , the two relative displacement components δ and δ can be considered together as .

(11)
The linear relation between the traction and relative displacement are illustrated in Figure 2. The critical values of the relative displacement at the elastic limit, δ and δ , can be calculated as / , (13) where T and T are the material strengths in the normal and shear directions, respectively. For isotropic solids, the material strengths are assumed to be the same, i.e., T T σ .

General Description
The general CZM used in the FEM is formulated between the FE surfaces in three dimensions. The FE formulation of the cohesive zone has two different forms in general, the widely used continuum CZM [28,31,44] as shown in Figure 3a and the discrete CZM [45][46][47] as shown in Figure   Figure 2. Linear relation between traction and relative displacement in (a) the normal direction and (b) the shear direction.
The critical values of the relative displacement at the elastic limit, δ n0 and δ s0 , can be calculated as where T n0 and T s0 are the material strengths in the normal and shear directions, respectively. For isotropic solids, the material strengths are assumed to be the same, i.e., T n0 = T s0 = σ c .

General Description
The general CZM used in the FEM is formulated between the FE surfaces in three dimensions. The FE formulation of the cohesive zone has two different forms in general, the widely used continuum CZM [28,31,44] as shown in Figure 3a and the discrete CZM [45][46][47] as shown in Figure 3b. The continuum CZM treats the fracture process zone between any two continuous FE surfaces as a cohesive element, and the cohesive traction is related to the displacement jump between the quadrature points on the separate surfaces by means of a defined cohesive law. Therefore, this cohesive formulation is surface-wise, and the cohesive traction between each quadrature point pair is dependent on the nodal displacements of the cohesive element. This surface-wise cohesive formulation may suffer from non-convergence [48], and great care must be taken with the numerical integration scheme for the cohesive elements [49,50]. By contrast, the discrete CZM constructs cohesive elements at adjacent nodes by means of a nonlinear spring [48]. The cohesive traction of this discrete spring-type formulation depends only on the displacement jump between the node pair and, thus, this cohesive formulation is point-wise. In contrast to the surface-wise cohesive elements, the point-wise cohesive elements have better convergence and are insensitive to mesh size and loading rate.
Materials 2020, 13, x FOR PEER REVIEW 6 of 35 a cohesive element, and the cohesive traction is related to the displacement jump between the quadrature points on the separate surfaces by means of a defined cohesive law. Therefore, this cohesive formulation is surface-wise, and the cohesive traction between each quadrature point pair is dependent on the nodal displacements of the cohesive element. This surface-wise cohesive formulation may suffer from non-convergence [48], and great care must be taken with the numerical integration scheme for the cohesive elements [49,50]. By contrast, the discrete CZM constructs cohesive elements at adjacent nodes by means of a nonlinear spring [48]. The cohesive traction of this discrete spring-type formulation depends only on the displacement jump between the node pair and, thus, this cohesive formulation is point-wise. In contrast to the surface-wise cohesive elements, the point-wise cohesive elements have better convergence and are insensitive to mesh size and loading rate.
. Figure 3. Schematic of different cohesive zone models in two dimensions: (a) continuum cohesive zone model with interaction through quadrature points; (b) discrete cohesive zone model with interaction through nodes; (c) particle-based cohesive zone model with interaction through DE particles.
In line with the point-wise cohesive formulation, the cohesive zone associated with the lattice DE model is inserted into each adjacent particle pair, shown in Figure 3c as a bridging zone between the DE connective model and the DE contact model. Each particle pair of the connective model constitutes a cohesive element and, thus, the CZM is formulated at microscale. It is evident from Figure 3c that the cohesive traction is only related to the displacement jump of the particle pair. Note that Gao and Klein [51] also used a cohesive formulation for material particles at microscale, but their formulation relied on a virtual internal bond model.
As shown in Figure 3c, the connection orientations of the particle-based CZM vary from pair to pair and, thus, the opening directions of crack fronts can never be uniform at microscale, which is different from the case of the continuum CZM ( Figure 3a) and the discrete CZM (Figure 3b).
A mixed-mode propagation criterion for the particle-based CZM is firstly proposed to describe the model transition from a cohesive model to a contact model, and then a mixed-mode initiation criterion is proposed for the model transition from a connective model to a cohesive model.

Mixed-Mode Fracture Propagation Criterion
Based on the fracture model originally developed by Tvergaard and Hutchinson [52] and modified by Espinosa and Zavattieri [53] and Song et al. [54], a potential-based mixed-mode formulation is proposed for the particle-based CZM at microscale. To formulate this mixed-mode cohesive relation, a non-dimensional scalar λ of the effective displacement jump is defined as In line with the point-wise cohesive formulation, the cohesive zone associated with the lattice DE model is inserted into each adjacent particle pair, shown in Figure 3c as a bridging zone between the DE connective model and the DE contact model. Each particle pair of the connective model constitutes a cohesive element and, thus, the CZM is formulated at microscale. It is evident from Figure 3c that the cohesive traction is only related to the displacement jump of the particle pair. Note that Gao and Klein [51] also used a cohesive formulation for material particles at microscale, but their formulation relied on a virtual internal bond model.
As shown in Figure 3c, the connection orientations of the particle-based CZM vary from pair to pair and, thus, the opening directions of crack fronts can never be uniform at microscale, which is different from the case of the continuum CZM ( Figure 3a) and the discrete CZM (Figure 3b).
A mixed-mode propagation criterion for the particle-based CZM is firstly proposed to describe the model transition from a cohesive model to a contact model, and then a mixed-mode initiation criterion is proposed for the model transition from a connective model to a cohesive model.

Mixed-Mode Fracture Propagation Criterion
Based on the fracture model originally developed by Tvergaard and Hutchinson [52] and modified by Espinosa and Zavattieri [53] and Song et al. [54], a potential-based mixed-mode formulation is proposed for the particle-based CZM at microscale. To formulate this mixed-mode cohesive relation, a non-dimensional scalar λ of the effective displacement jump is defined as where < > is the Macaulay bracket defined by < >= (( ) + | |)/2, δ n and δ s are current opening and shear separation components, and δ nc and δ sc are critical values for opening mode and shear mode at complete cracking points, respectively. Equation (14) indicates that both opening displacement and shear separation displacement are taken into account for the evaluation of fracture propagation initiation. Note that the use of the Macaulay bracket on the opening crack indicates that only the tensional fracture is accounted for, while the fracture due to compression is not considered in the present cohesive law. The resistance to compression is formulated as a penalty-like force with k n as the penalty stiffness.
An effective traction T(λ) is constructed in the form of linear softening as follows: where the cohesive strength σ c is the material tensile strength. As illustrated in Figure 4, this traction reduces linearly from the material cohesive strength at the fracture initiation point λ cr to zero at the decohesion point when Equation (16)  where 〈∎〉 is the Macaulay bracket defined by 〈∎〉 ∎ |∎| 2 ⁄ , δ and δ are current opening and shear separation components, and δ and δ are critical values for opening mode and shear mode at complete cracking points, respectively. Equation (14) indicates that both opening displacement and shear separation displacement are taken into account for the evaluation of fracture propagation initiation. Note that the use of the Macaulay bracket on the opening crack indicates that only the tensional fracture is accounted for, while the fracture due to compression is not considered in the present cohesive law. The resistance to compression is formulated as a penalty-like force with k as the penalty stiffness.
An effective traction T λ is constructed in the form of linear softening as follows: , where the cohesive strength σ is the material tensile strength. As illustrated in Figure 4, this traction reduces linearly from the material cohesive strength at the fracture initiation point λ to zero at the decohesion point when Equation (16) is satisfied.  The fracture initiation point λ denotes a critical value of the effective displacement jump at which the effective traction T λ reaches the maximum. λ is determined by the initiation criterion introduced in next section. Note that the critical value λ is dependent on the physical parameters of the particle-based CZM rather than a user-defined value as used to adjust initial loading stiffness [53,54].
Using the effective displacement jump and the effective traction, a potential function Φ λ is constructed as follows [52]: . (17) Furthermore, the total dissipated energy within the linear softening stage per unit area of newly created fracture surfaces is defined as the critical energy release rate , which is given by The fracture initiation point λ cr denotes a critical value of the effective displacement jump at which the effective traction T(λ) reaches the maximum. λ cr is determined by the initiation criterion introduced in next section. Note that the critical value λ cr is dependent on the physical parameters of the particle-based CZM rather than a user-defined value as used to adjust initial loading stiffness [53,54].
Using the effective displacement jump and the effective traction, a potential function Φ(λ) is constructed as follows [52]: Furthermore, the total dissipated energy within the linear softening stage per unit area of newly created fracture surfaces is defined as the critical energy release rate G c , which is given by The first derivative of the potential function yields the opening component T n and shear component T s of the cohesive traction as follows: Apparently, both traction components, T n and T s , are dependent on two separation components, δ n and δ s . When λ = 1, both components vanish without any extra enforcements, indicating that the decohesion at both modes is completed simultaneously. This is consistent with the physical phenomenon [55]. Note that this particular feature is sought in formulating the initiation criterion in next section. Once λ ≥ 1, fracture propagates to the next particle pair, and the current particle pair is transitioned to be in a contact condition. The specific treatment on contact is illustrated in Section 4.
Since the fracture process is irreversible, the irreversibility of the cohesive law is realized by recording the maximum λ as The damage can then be defined by a scalar as As λ * is monotonically increasing, the history of damage, as shown by Equations (21) and (22), is monotonic as well, and this is consistent with the damage mechanics. This damage definition provides a good description of the damage state as the damage process varies in a linear form, starting from 0 for a state without damage to 1 for a state with full damage [56].
The unloading is assumed to be back to the origin as shown in Figure 4, and the reloading is assumed to be back to the original point as determined by λ * without any energy dissipation. During the course of unloading/reloading, the loading ratio is defined as Then, the effective traction can be generally rewritten as It is apparent that the effective traction relies on the damage history as defined by Equation (23). Furthermore, the traction components as shown in Equation (19) can be rewritten in general form as follows: Materials 2020, 13, 3573 9 of 35 Therefore, two different shear traction components can be determined by T s1 = T s δ s1 δ s and T s2 = T s δ s2 δ s . (27) Note that these traction components are determined in the local frame.
The traction vector T is used to determine the local cohesive forcef coh aŝ The cohesive force f coh at the global frame can then be determined using the transformation matrix φ Equation (4) as Moreover, the moment m coh related to the cohesive force is obtained as where r denotes an effective radius vector. One limitation of the cohesive law is that the critical energy release rates of the opening mode G Ic and the shear mode G IIc are required to be the same in any mixed-mode conditions [52], i.e., G Ic = G IIc = G c . This is because a single G c is used to define the critical energy release rate as shown in Equation (18). Note that the current DE model is only applicable to homogeneous and isotropic elastic solids. For isotropic materials, the critical energy release rates at the two modes can be reasonably assumed to be the same [55,57,58]. For an extension to two different critical energy release rates, the reader can refer to Reference [53].

Mixed-Mode Fracture Initiation Criterion
Analogous to the formulation of the propagation criterion, a non-dimensional scalar ψ, which is similar to the scalar λ in Equation (14), is defined as It is used to develop a mixed-mode initiation criterion for the onset of fracture initiation in the displacement space. Once ψ = 1, fracture initiates, and Figure 5 describes this ellipse initiation criterion. One limitation of the cohesive law is that the critical energy release rates of the opening mode and the shear mode are required to be the same in any mixed-mode conditions [52], i.e., . This is because a single is used to define the critical energy release rate as shown in Equation (18). Note that the current DE model is only applicable to homogeneous and isotropic elastic solids. For isotropic materials, the critical energy release rates at the two modes can be reasonably assumed to be the same [55,57,58]. For an extension to two different critical energy release rates, the reader can refer to Reference [53].

Mixed-Mode Fracture Initiation Criterion
Analogous to the formulation of the propagation criterion, a non-dimensional scalar ψ, which is similar to the scalar λ in Equation (14), is defined as It is used to develop a mixed-mode initiation criterion for the onset of fracture initiation in the displacement space. Once ψ 1, fracture initiates, and Figure 5 describes this ellipse initiation criterion. This initiation criterion was also used by Alfano and Crisfield [55], Mi et al. [57], and Qiu et al. [59] to design a mixed-mode fracture model with the intention of avoiding the pre-solution of the mixed-mode ratio δ δ . A pre-solution of the mode ratio or a known mode ratio is generally needed by a number of mixed-mode formulations [56,[60][61][62]. These formulations are based on the energy release rate and use the power law criteria [63,64] or the B-K criterion [65] as the fracture This initiation criterion was also used by Alfano and Crisfield [55], Mi et al. [57], and Qiu et al. [59] to design a mixed-mode fracture model with the intention of avoiding the pre-solution of the mixed-mode ratio δ s /δ n . A pre-solution of the mode ratio or a known mode ratio is generally needed by a number of mixed-mode formulations [56,[60][61][62]. These formulations are based on the energy release rate and use the power law criteria [63,64] or the B-K criterion [65] as the fracture propagation criterion. This kind of formulation also requires the mode ratio to be constant during the course of the fracture process. The mode ratio, however, changes along the fracture process zone [62] or throughout the loading [66], which might be due to the displacement jump [67] or impact events [68].
The change of the mixed-mode ratio is illustrated in Figure 5, where the mode ratio determined at point A can possibly change to a different value when the loading equilibrium moves to point B at the next time step upon a displacement jump, for instance. According to Pinho et al. [66], once the mode ratio changes, the damage history and the strategy to restore the cohesive energy are not clear. A solution was also proposed, where they resorted to the maximum mixed-mode displacement and employed a circle propagation criterion for the mixed-mode propagation criterion, regardless of the mode ratio.
This idea was adopted as the present initiation criterion. When associated with the non-dimensional scalar ψ, Equation (31) is rewritten as with δ 0 = min{δ n0 , δ s0 }. The maximum of ψ is recorded as When ψ * = 1, fracture initiates. Furthermore, two separation components can then be used to calculate the critical value λ cr along with Equation (14). Note that this strategy in recording the maximum value to indicate fracture initiation was also employed in Reference [55].
Therefore, a circle initiation criterion is developed, as illustrated in Figure 5. It should be mentioned that this initiation criterion is similar to the propagation criterion [66] in the sense that both opening and shear modes can be activated simultaneously, which is more consistent with the fracture mechanism. Additionally, this treatment takes into account the concern of safe design, which was employed by Balzani and Wagner [61] for a propagation criterion. Furthermore, this initiation criterion essentially indicates the onset of the fracture process, and it plays a similar role to that used in extrinsic cohesive zone models [31,32,69].
To sum up, before ψ * reaches unity, the reversible connective model works; after satisfying this initiation criterion, the irreversible cohesive model comes into effect.

Contact Model: Representation of Particulate Materials
Once decohesion in a particle pair is completed, these two particles become in contact. The contact force is calculated based on the Hertz contact theory [70][71][72], which describes the relation between the contact force and the penetration between two particles of a pair.
When contact occurs between the two particles i and j, as shown in Figure 6, the Hertz contact model is used to calculate the contact force on particle i along the normal direction as where n denotes a unit normal vector, κ denotes a user-defined penalty factor, and E is an equivalent Young's modulus, as determined by where ν i and ν j are the Poisson's ratios of the two contacting particles, and r is the equivalent radius, given by where denotes a unit normal vector, κ denotes a user-defined penalty factor, and E is an equivalent Young's modulus, as determined by , (35) where ν and ν are the Poisson's ratios of the two contacting particles, and r is the equivalent radius, given by ̃ .

(36)
Lastly, δ is the penetration depth as determined by , (37) where the distance d between the centers of two particles can be calculated with the use of two position vectors as • .
(38)  Lastly, δ n is the penetration depth as determined by where the distance d between the centers of two particles can be calculated with the use of two position vectors as When contact is persistent, Mindlin's theory [71,73], which is based on the assumption that there is no partial slip between contacting particles, is used to calculate the tangential contact force on particle i as where G is the equivalent shear modulus as determined by G is the shear modulus, determined as G = E/2(1 + ν), while δ s denotes relative tangential displacement during the period τ con of persistent contact, expressed as where relative tangential velocity v c(s) at the contact point c is associated with the relative velocity v c , determined by v c(s) = v c − (v c ·n)n.
The relative velocity v c at the contact point c can be determined as where where v i and v j are translational velocities of the particles, ω i and ω j are angular velocities of the particles, and r i and r j are effective radius vectors of the particles, as shown in Figure 6.

Model Transition: Monotonically from Connection to Contact
As fracture initiates and propagates, the model can transition from a connective model to a cohesive model and further to a contact model. The DE connective model, representing elastic solids at microscale, indicates the phase of a continuum. The contact model represents particulate materials. The particle-based CZM bridges the gap between them and describes the transition process of materials. The three models and the transition in between are illustrated in Figure 7. As the fracture process is irreversible, the model transition is monotonic toward the contact model. This monotonic model transition is consistent with the material transition. con , (41) where relative tangential velocity at the contact point c is associated with the relative velocity , determined by • . (42) The relative velocity at the contact point c can be determined as , (43) where , , (45) where and are translational velocities of the particles, and are angular velocities of the particles, and ̅ and ̅ are effective radius vectors of the particles, as shown in Figure 6.

Model Transition: Monotonically from Connection to Contact
As fracture initiates and propagates, the model can transition from a connective model to a cohesive model and further to a contact model. The DE connective model, representing elastic solids at microscale, indicates the phase of a continuum. The contact model represents particulate materials. The particle-based CZM bridges the gap between them and describes the transition process of materials. The three models and the transition in between are illustrated in Figure 7. As the fracture process is irreversible, the model transition is monotonic toward the contact model. This monotonic model transition is consistent with the material transition. Specific transition procedures are shown in Table 1. Initially, all particle pairs are in a connective relation and are employed to represent solids. The scalar ψ in the fracture initiation Specific transition procedures are shown in Table 1. Initially, all particle pairs are in a connective relation and are employed to represent solids. The scalar ψ in the fracture initiation criterion of Equation (32) is used as the first transition criterion, and, once ψ * ≥ 1.0, the DE model changes its phase and becomes a cohesive model. Similarly, the scalar λ in the fracture propagation criterion of Equation (14) is used as the second transition criterion, and, once λ * ≥ 1.0, the DE particles become in contact condition. Table 1. Transition procedures for discrete element (DE) particle pairs.

Initial: All Pairs Being Connective
If (connective and ψ * ≥ 1.0) then connective Materials 2020, 13, x FOR PEER REVIEW 12 of 35 criterion of Equation (32) is used as the first transition criterion, and, once ψ * 1.0, the DE model changes its phase and becomes a cohesive model. Similarly, the scalar λ in the fracture propagation criterion of Equation (14) is used as the second transition criterion, and, once λ * 1.0, the DE particles become in contact condition.  Figure 7b, a particle might be in different interaction manners with different surrounding particles in the meantime, while a particle pair, e.g., particle and particle , can only exhibit one certain interaction behavior, whether connective, cohesive, or contact.

Implementation: Explicit Update of Kinematics
The equations of motion of DE particles are expressed as where m and I are the particle mass and the moment of inertia, respectively, and are cohesive else if (cohesive and λ * ≥ 1.0) then cohesive Materials 2020, 13, x FOR PEER REVIEW 12 of 35 criterion of Equation (32) is used as the first transition criterion, and, once ψ * 1.0, the DE model changes its phase and becomes a cohesive model. Similarly, the scalar λ in the fracture propagation criterion of Equation (14) is used as the second transition criterion, and, once λ * 1.0, the DE particles become in contact condition.  Figure 7b, a particle might be in different interaction manners with different surrounding particles in the meantime, while a particle pair, e.g., particle and particle , can only exhibit one certain interaction behavior, whether connective, cohesive, or contact.

Update of Kinematics
The equations of motion of DE particles are expressed as It should be mentioned that, as shown in Figure 7b, a particle might be in different interaction manners with different surrounding particles in the meantime, while a particle pair, e.g., particle i and particle j, can only exhibit one certain interaction behavior, whether connective, cohesive, or contact.

Implementation: Explicit Update of Kinematics
The equations of motion of DE particles are expressed as where m and I are the particle mass and the moment of inertia, respectively, v and ω are translational and rotational velocities, respectively, and ω are their first derivatives with respect to time; f ext denotes external force and may include the contribution of coupling force if this particle is in a coupling interaction with the FE domain, and m denotes the moment induced from f ext ; f int , f coh , and f con are the internal force, cohesive force, and contact force contributed by the above-mentioned DE connective model, cohesive model, and contact model, respectively; N i is the number of surrounding particles around particle i; r denotes an effective radius and, in Figure 7c, it is a vector from the center of particle i to point c, which is located at the common plane between particle i and particle j; lastly, f is f int , f coh , or f con because a specific particle pair can be in connection, in cohesion, or in contact.
The calculations of f int , f coh , or f con are the main concern here, and the specific procedures are briefly summarized in Table 2. Note that tractions can be transformed to forces. Table 2. Computational procedures of the interaction force between a particle pair.

If (Connective) then
Use Equations (1) and (5) to calculate internal force f int else if (cohesive) then Use Equations (25), (26) and (29) to calculate cohesive force f coh else if (contact) then Use Equations (34) and (39) to calculate contact force f con end if The central difference method was adopted here for explicit time integration. Since the central difference method is conditionally stable, the selected time step should be less than the critical one. For the selection strategy of the time step, the reader can refer to [74] for detail. Velocities and relative displacement are assumed to be known at time t (n) . At any time, the translational acceleration and rotational acceleration can be calculated from Equations (46) and (47), respectively. Supposing the incremental time step from time t (n) to t (n+1) is ∆t, the translational and rotational velocities at the mid-point of that period, i.e., t (n+1/2) , can be respectively calculated by The relative velocity v c at the contact point c, which is on the common plane (see Figure 7c), between a particle pair (i and j) can, therefore, be calculated as where Then, the incremental relative displacement vector ∆δ (n+1) from time t (n) to t (n+1) can be updated by and the total relative displacement vector δ (n+1) at time t (n+1) is updated as which is essential for the calculation of internal, cohesive, and contact forces.

Numerical Simulations
The particle-based cohesive crack model was validated using three numerical simulations of standard fracture tests and applied to the impact fracture of a notched concrete beam. These standard fracture tests included a double cantilever beam (DCB) model, an end notched flexural (ENF) model, and a mixed-mode bending (MMB) model. These three fracture tests correspond to mode-I (opening mode), mode-II (shear mode), and a mixed mode between I and II, respectively.
It is assumed that the material's macroscopic behavior is determined by the material's micro-structure and the interplay between particles at its associated scale. Therefore, the proposed microscopic cohesive crack model can be validated using the macro-structural loading-deflection relation of these fracture tests, which can be analytically determined based on the linear elastic fracture mechanics (LEFM). For these three models, cracks occur only on their mid-plane and, thus, cohesive elements are pre-defined as a weak layer.
The primary material properties for the isotropic solids were adopted from References [62,75], where E = 120 GPa, G Ic = 0.26 N/mm for the opening mode, and G IIc = 1.002 N/mm for the shear mode. In addition, density ρ = 2500 kg/m 3 , Poisson's ratio ν = 0.2, and cohesive strength σ c = 70 MPa. The DE particle radius was set to r = 0.125 mm.
In the numerical simulations, loading was in displacement control with a velocity prescribed at selected DE particles. The prescribed velocity v went up to a certain value v 0 within a duration t 0 , as expressed in the following form: Note that v 0 = 0.075 m/s and t 0 = 0.5 ms apply to the numerical simulations of DCB, ENF, and MMB models; otherwise, an explicit statement would be made. It also should be mentioned that these simulations are for the quasi-static loading condition. Hence, the dynamic relaxation strategy via global damping is employed to make sure that the kinetic energy is less than 0.001% of the total energy [68].
The definition of the critical energy release rate on the DE model needs to be modified because of its special spatial structure. According to Irwin's energy approach for fracture [76], the energy release rate is a measure of the energy available for an increment of crack extension. On the basis of the linear elastic fracture mechanics (LEFM), this energy release rate is also equal to its critical value G c . Therefore, for the present beam analysis with a crack width B and length a, the critical energy release rate can be determined in association with the potential Π of the beam as follows [9]: Note that the potential is related to the strain energy and the external work on the beam, and the kinematic energy is ignored for the quasi-static loading condition [9].
Considering the geometry of the present DE model, as shown in Figure 1b, it is found that the boundary DE particles in the model only have a connection with bulk particles and lose half their interaction on the other side. Thus, the width B needs to be corrected to B − 2r when used in Equation (49). This effect, however, diminishes with the increase of width.

Mode-I Validation: Numerical Analysis of a DCB Model
The DCB model shown in Figure 8 is widely used to test fracture behavior in opening mode. The dimensions of the model were as follows: length L = 30 mm, breadth B = 2 mm, height 2h = 2 mm, and initial crack length a 0 = 9 mm. According to the beam theory, the general form of the beam load-deflection relation is given as where E is the Young's modulus, ℐ = Bh 3 /12 is the second moment of area of one arm of the cantilever, and χ I is a mode-I correction factor for calculating an effective crack length to account for the root rotation and the shear deformation of the beam [77][78][79]. For isotropic solids used in the present study, this correction factor χ I associated with mode-I fracture is estimated as follows [77]: where shear modulus G = E/2(1 + ν), and Γ is the elastic modulus correction parameter, which is given by . Therefore, for the present beam analysis with a crack width B and length a, the critical energy release rate can be determined in association with the potential Π of the beam as follows [9]: . (56) Note that the potential is related to the strain energy and the external work on the beam, and the kinematic energy is ignored for the quasi-static loading condition [9].
Considering the geometry of the present DE model, as shown in Figure 1b, it is found that the boundary DE particles in the model only have a connection with bulk particles and lose half their interaction on the other side. Thus, the width B needs to be corrected to B 2r when used in Equation (49). This effect, however, diminishes with the increase of width.

Mode-I Validation: Numerical Analysis of a DCB Model
The DCB model shown in Figure 8 is widely used to test fracture behavior in opening mode. The dimensions of the model were as follows: length L 30 mm, breadth B 2 mm, height 2h 2 mm, and initial crack length a 9 mm. According to the beam theory, the general form of the beam load-deflection relation is given as where E is the Young's modulus, ℐ Bh 12 ⁄ is the second moment of area of one arm of the cantilever, and χ is a mode-I correction factor for calculating an effective crack length to account for the root rotation and the shear deformation of the beam [77][78][79]. For isotropic solids used in the present study, this correction factor χ associated with mode-I fracture is estimated as follows [77]: where shear modulus G E/2 1 ν , and is the elastic modulus correction parameter, which is given by When the crack propagates beyond the initial crack length, i.e., a a , the crack length is associated with the total energy release rate . For the fracture problem of mode-I, . Therefore, the crack length can be obtained by a B 2r Eℐ P ⁄ χ h based on the corrected beam theory [80]. Then, the line BCD (stage-II) in Figure 9 can be described by The initial loading line OA (stage-I) in Figure 9 corresponds to a stationary crack with an initial crack length a = a 0 and, thus, Equation (57) can be rewritten as When the crack propagates beyond the initial crack length, i.e., a > a 0 , the crack length is associated with the total energy release rate G T . For the fracture problem of mode-I, G T = G Ic . Therefore, the crack length can be obtained by a = (B − 2r)EG Ic /P − χ I h based on the corrected beam theory [80]. Then, the line BCD (stage-II) in Figure 9 can be described by Materials 2020, 13, x FOR PEER REVIEW 15 of 35 After the crack propagates to the end of the initial weak layer, i.e., a L, each cantilever beam behaves independently. Therefore, the line OE (stage-III) in Figure 9 is described by the following equation according to the beam theory: In numerical simulations of the DCB model where , loads are in displacement control along both upward and downward directions, as seen in Figure 8. Load is applied to both top and bottom DE particles on the left side. After crack propagation is completed, the final configuration can be seen from Figure 10. The numerical results agree with the analytical solution at all three stages, as can be seen from Figure 9. Slight noise is observed since the fracture model was implemented in an explicit time-integration code. Different cohesive strengths were employed to investigate the load-deflection behavior as shown in Figure 9, where the difference was generally indistinguishable except for the peak loads. It is understood from Equation (18) that, for a given , the size of the cohesive zone decreases with increasing cohesive strength, which results in materials tending to be more brittle. As a result, numerical results can better match the LEFM-based analytical solution.
The influence of the critical energy release rate on the load-deflection relation was investigated, and the results are shown in Figure 11. Numerical results computed from the proposed model are in good agreement with the analytical solution in all three stages for all three cases, i.e., 0.13 N/mm 2 , 0.26 N/mm 2 , and 0.39 N/mm 2 . With the increase of , only the line in the second stage changes, as determined by Equation (54). Since the initial loading stage and the complete split stage are with a stationary crack, the load-deflection relation at these two stages is irrelevant to . This is predicted by Equations (53) and (55) and confirmed by Figure 11. After the crack propagates to the end of the initial weak layer, i.e., a = L, each cantilever beam behaves independently. Therefore, the line OE (stage-III) in Figure 9 is described by the following equation according to the beam theory: In numerical simulations of the DCB model where G c = G Ic , loads are in displacement control along both upward and downward directions, as seen in Figure 8. Load is applied to both top and bottom DE particles on the left side. After crack propagation is completed, the final configuration can be seen from Figure 10.  The numerical results agree with the analytical solution at all three stages, as can be seen from Figure 9. Slight noise is observed since the fracture model was implemented in an explicit time-integration code. Different cohesive strengths were employed to investigate the load-deflection behavior as shown in Figure 9, where the difference was generally indistinguishable except for the peak loads. It is understood from Equation (18) that, for a given G c , the size of the cohesive zone decreases with increasing cohesive strength, which results in materials tending to be more brittle. As a result, numerical results can better match the LEFM-based analytical solution.
The influence of the critical energy release rate G Ic on the load-deflection relation was investigated, and the results are shown in Figure 11. Numerical results computed from the proposed model are in good agreement with the analytical solution in all three stages for all three cases, i.e., G Ic = 0.13 N/mm 2 , 0.26 N/mm 2 , and 0.39 N/mm 2 . With the increase of G Ic , only the line in the second stage changes, as determined by Equation (54). Since the initial loading stage and the complete split stage are with a stationary crack, the load-deflection relation at these two stages is irrelevant to G Ic . This is predicted by Equations (53) and (55) and confirmed by Figure 11. The size effect from different DE particle radii was examined, and results are shown in Figure  12. The results with r 0.10 mm match the results with r 0.125 mm, as well as the analytical results. This indicates that the size effect is insignificant. The loading rate's effect on the loaddeflection (~∆) relation is shown in Figure 13. The numerical results are in good agreement with the analytical solution, but they slightly deviate at the later second stage. It is also found that a larger loading rate leads to a larger discrepancy. The discrepancy is probably due to the fact that, when a larger loading rate is adopted in an explicit time-integration code, a larger inertia effect should be accounted for, which is against the quasi-static assumption of the LEFM-based analytical solution. The size effect from different DE particle radii was examined, and results are shown in Figure 12. The results with r = 0.10 mm match the results with r = 0.125 mm, as well as the analytical results. This indicates that the size effect is insignificant. The loading rate's effect on the load-deflection ( P ∼ ∆) relation is shown in Figure 13. The numerical results are in good agreement with the analytical solution, but they slightly deviate at the later second stage. It is also found that a larger loading rate leads to a larger discrepancy. The discrepancy is probably due to the fact that, when a larger loading rate is adopted in an explicit time-integration code, a larger inertia effect should be accounted for, which is against the quasi-static assumption of the LEFM-based analytical solution. The size effect from different DE particle radii was examined, and results are shown in Figure  12. The results with r 0.10 mm match the results with r 0.125 mm, as well as the analytical results. This indicates that the size effect is insignificant. The loading rate's effect on the loaddeflection (~∆) relation is shown in Figure 13. The numerical results are in good agreement with the analytical solution, but they slightly deviate at the later second stage. It is also found that a larger loading rate leads to a larger discrepancy. The discrepancy is probably due to the fact that, when a larger loading rate is adopted in an explicit time-integration code, a larger inertia effect should be accounted for, which is against the quasi-static assumption of the LEFM-based analytical solution. The CZM was constructed at the particle bonding springs. It is found from Equation (3) that the micromechanical parameters of spring stiffness k n and k s are correlated to the macroscopic material parameters of Young's modulus E and Poisson's ratio ν, apart from particle radius r. Hence, the effect of the variation of E and ν was also examined. The effect of Young's modulus E on the load-deflection relation is shown in Figure 14, where it can be seen that numerical results computed from the proposed model agree well with the analytical solutions for all three stages. It can be seen from Equation (3) that the increase of macroscopic E results in an increase of the micro-structural spring stiffness. The change of micromechanical parameters leads to varying macroscopic performance, albeit consistent with the respective analytical predictions. Moreover, Figure 14 shows a greater E leading to a higher peak value and a shorter fracture process, which is consistent with the prediction by Equations (60)- (62). The CZM was constructed at the particle bonding springs. It is found from Equation (3) that the micromechanical parameters of spring stiffness k and k are correlated to the macroscopic material parameters of Young's modulus E and Poisson's ratio ν, apart from particle radius . Hence, the effect of the variation of E and ν was also examined. The effect of Young's modulus E on the load-deflection relation is shown in Figure 14, where it can be seen that numerical results computed from the proposed model agree well with the analytical solutions for all three stages. It can be seen from Equation (3) that the increase of macroscopic E results in an increase of the micro-structural spring stiffness. The change of micromechanical parameters leads to varying macroscopic performance, albeit consistent with the respective analytical predictions. Moreover, Figure 14 shows a greater E leading to a higher peak value and a shorter fracture process, which is consistent with the prediction by Equations (60)-(62). The effect of Poisson's ratio ν on the load-deflection relation is shown in Figure 15. The computed results with various values of ν are generally in good agreement with the analytical solution, particularly the results on fracture propagation in stage-II. For stage-I and stage-III, it is observed that a smaller magnitude of ν leads to a slightly later start of the fracture process, as well as a later closure. This behavior indicates that the beam with a smaller ν possesses a slightly lower  The CZM was constructed at the particle bonding springs. It is found from Equation (3) that the micromechanical parameters of spring stiffness k and k are correlated to the macroscopic material parameters of Young's modulus E and Poisson's ratio ν, apart from particle radius . Hence, the effect of the variation of E and ν was also examined. The effect of Young's modulus E on the load-deflection relation is shown in Figure 14, where it can be seen that numerical results computed from the proposed model agree well with the analytical solutions for all three stages. It can be seen from Equation (3) that the increase of macroscopic E results in an increase of the micro-structural spring stiffness. The change of micromechanical parameters leads to varying macroscopic performance, albeit consistent with the respective analytical predictions. Moreover, Figure 14 shows a greater E leading to a higher peak value and a shorter fracture process, which is consistent with the prediction by Equations (60)-(62). The effect of Poisson's ratio ν on the load-deflection relation is shown in Figure 15. The computed results with various values of ν are generally in good agreement with the analytical solution, particularly the results on fracture propagation in stage-II. For stage-I and stage-III, it is observed that a smaller magnitude of ν leads to a slightly later start of the fracture process, as well as a later closure. This behavior indicates that the beam with a smaller ν possesses a slightly lower The effect of Poisson's ratio ν on the load-deflection relation is shown in Figure 15. The computed results with various values of ν are generally in good agreement with the analytical solution, particularly the results on fracture propagation in stage-II. For stage-I and stage-III, it is observed that a smaller magnitude of ν leads to a slightly later start of the fracture process, as well as a later closure. This behavior indicates that the beam with a smaller ν possesses a slightly lower macro-structural bending stiffness. The discrepancies in stage-I and stage-III are due to the spring stiffness. The stiffness of springs changes with the variation of ν and, thus, it leads to a slightly different macro-structural response on the current beam, which uses a limited number of DE particles (4096 particles on each beam). Nonetheless, this discrepancy may be reduced with the increasing number of DE particles used in the beam.
stiffness. The stiffness of springs changes with the variation of ν and, thus, it leads to a slightly different macro-structural response on the current beam, which uses a limited number of DE particles (4096 particles on each beam). Nonetheless, this discrepancy may be reduced with the increasing number of DE particles used in the beam.

Mode-II Validation: Numerical Analysis of an ENF Model
The ENF model, as shown in Figure 16, is widely used to simulate fracture processes under shear-mode loading. The dimensions of the model were as follows: length 2L 30.25 mm, breadth B 2 mm, height 2h 2 mm, and initial crack length a 9 mm. In the context of simple beam theory [57,81], the analytical solutions for the line OA and the line BCD in Figure 17 take the same form. When considering the corrected beam theory [78], the analytical form can be written as where χ is a correction factor similar to χ , associated with mode-II fracture. This factor only takes into account the shear factor and does not include root rotation [82]; hence, its value is less than that used for the DCB model. It is estimated to be χ 0.42χ [82].

Mode-II Validation: Numerical Analysis of an ENF Model
The ENF model, as shown in Figure 16, is widely used to simulate fracture processes under shear-mode loading. The dimensions of the model were as follows: length 2L = 30.25 mm, breadth B = 2 mm, height 2h = 2 mm, and initial crack length a 0 = 9 mm. In the context of simple beam theory [57,81], the analytical solutions for the line OA and the line BCD in Figure 17 take the same form. When considering the corrected beam theory [78], the analytical form can be written as where χ II is a correction factor similar to χ I , associated with mode-II fracture. This factor only takes into account the shear factor and does not include root rotation [82]; hence, its value is less than that used for the DCB model. It is estimated to be χ II = 0.42χ I [82].
macro-structural bending stiffness. The discrepancies in stage-I and stage-III are due to the spring stiffness. The stiffness of springs changes with the variation of ν and, thus, it leads to a slightly different macro-structural response on the current beam, which uses a limited number of DE particles (4096 particles on each beam). Nonetheless, this discrepancy may be reduced with the increasing number of DE particles used in the beam.

Mode-II Validation: Numerical Analysis of an ENF Model
The ENF model, as shown in Figure 16, is widely used to simulate fracture processes under shear-mode loading. The dimensions of the model were as follows: length 2L 30.25 mm, breadth B 2 mm, height 2h 2 mm, and initial crack length a 9 mm. In the context of simple beam theory [57,81], the analytical solutions for the line OA and the line BCD in Figure 17 take the same form. When considering the corrected beam theory [78], the analytical form can be written as where χ is a correction factor similar to χ , associated with mode-II fracture. This factor only takes into account the shear factor and does not include root rotation [82]; hence, its value is less than that used for the DCB model. It is estimated to be χ 0.42χ [82]. For the initial loading line OA (stage-I) in Figure 17, the crack is stationary and a a . Therefore, the analytical solution of the load-deflection relation can be obtained straightforwardly after substituting a into Equation (63). By contrast, the crack length a for the line BCD in Figure 17 is related to the critical energy release rate , i.e., a 64 B 2r Eℐ /3P χ h [57]. When For the initial loading line OA (stage-I) in Figure 17, the crack is stationary and a = a 0 . Therefore, the analytical solution of the load-deflection relation can be obtained straightforwardly after substituting a 0 into Equation (63). By contrast, the crack length a for the line BCD in Figure 17 is related to the critical energy release rate G IIc , i.e., a = 64(B − 2r)EG IIc /3P 2 − χ II h [57]. When substituting it into Equation (63), the analytical solution of the line BCD (stage-II) is obtained. Note that, at this stage, the crack length is limited to the mid-span, i.e., a < L.
When the crack propagates beyond the mid-span, i.e., a > L, the analytical solution of the line EF (stage-III) in Figure 17 is given by Equation (64) [81].
Once the beams are completely split, the load-deflection relation of the line OG (stage-IV) is simply given by Equation (65) [57].
In numerical simulations of the ENF model where G c = G IIc , load P is applied at the mid-span, and deflection ∆ is also measured there. The deformed configuration is shown in Figure 18. The overall performance agrees very well with the analytical prediction. The crack propagation is completed after the deflection of the mid-span ∆ exceeds 2.0 mm. Then, the two beams deform as an entity, and the load-deflection relation matches the beam theory, which is described by the line OG. Different cohesive strengths are used, and the difference can only be observed from the transition between stage-I and stage-II. It is not surprising to see that a higher cohesive strength results in a better agreement since materials behave in a more brittle way. Additionally, it is worth noting that, once cracks start initiate, the fracture behavior associated with the present particle-based CZM is generally insensitive to the cohesive strength variation, which was also observed by Harper and Hallet [75].   The influence of initial crack length on the load-deflection relation can be seen from Figure 19. Various initial crack lengths (a 0 = 6 mm, 9 mm, and 12 mm) were adopted and, overall, the computational results from the proposed model are in good agreement with the analytical ones. Obviously, with a shorter initial crack length, the beam is of a higher bending stiffness and, thus, the initial loading tends to be more resistant. This is consistent with the beam theory. It should be mentioned that, if the initial crack is too short (a 0 = 6 mm in this case), the instability phenomenon of snap-back occurs (see Figure 19). Because the present loading is in a displacement control, the mid-span deflection would monotonically increase, which violates the analytical prediction that deflection can decrease, as shown in Figure 19. This dynamic snap-back was also postulated to exist in reality by Mi et al. [57], consistent with the experimental observation [83] that instability takes place when ENF tests are conducted with an insufficient initial crack length. It should also be noted that, even though snap-back occurs, the remaining load-deflection relation still follows the correct path until both beams completely split. in reality by Mi et al. [57], consistent with the experimental observation [83] that instability takes place when ENF tests are conducted with an insufficient initial crack length. It should also be noted that, even though snap-back occurs, the remaining load-deflection relation still follows the correct path until both beams completely split. The influence of the critical energy release rate on the load-deflection relation can be seen from Figure 20, where results in stage-II and stage-III are affected but remain almost unaffected in stage-I and stage-IV. This observation is in agreement with the analytical prediction, as seen in Figure 20 and the theory of Equations (63) and (64), since only the crack length a is affected. Furthermore, it is confirmed that the critical energy release rate plays a significant role in fracture behavior. The influence of the critical energy release rate on the load-deflection relation can be seen from Figure 20, where results in stage-II and stage-III are affected but remain almost unaffected in stage-I and stage-IV. This observation is in agreement with the analytical prediction, as seen in Figure 20 and the theory of Equations (63) and (64), since only the crack length a is affected. Furthermore, it is confirmed that the critical energy release rate plays a significant role in fracture behavior.

Mixed-Mode I & II Validation: Numerical Analysis of an MMB Model
The MMB test (shown in Figure 21) was designed by Reeder and Crew [84] to provide a wide range of mixed-mode ratios by adjusting the length of the loading arm. The dimensions of the model were as follows: length 2L 30.25 mm, breadth B 2 mm, height 2h 2 mm, and initial crack length a 9 mm. An MMB model is equivalent to the superposition of a DCB model and an ENF model. Therefore, the analytical solution of the MMB model can be attained via superposition. Given the loading P and the loading arm C, as shown in Figure 21, the pure mode-I loading P and mode-II loading P can be determined by Equation (66) [84]. ; .

Mixed-Mode I & II Validation: Numerical Analysis of an MMB Model
The MMB test (shown in Figure 21) was designed by Reeder and Crew [84] to provide a wide range of mixed-mode ratios by adjusting the length of the loading arm. The dimensions of the model were as follows: length 2L = 30.25 mm, breadth B = 2 mm, height 2h = 2 mm, and initial crack length a 0 = 9 mm. An MMB model is equivalent to the superposition of a DCB model and an ENF model. Therefore, the analytical solution of the MMB model can be attained via superposition. Given the loading P and the loading arm C, as shown in Figure 21, the pure mode-I loading P I and mode-II loading P II can be determined by Equation (66) [84].
model. Therefore, the analytical solution of the MMB model can be attained via superposition. Given the loading P and the loading arm C, as shown in Figure 21, the pure mode-I loading P and mode-II loading P can be determined by Equation (66)  The mixed-mode ratio β / , where and and are energy release rates corresponding to the opening mode and shear mode, respectively. This ratio changes with the adjustment of the loading arm, and the specific relation is expressed by Equation (67) [84]. The mixed-mode ratio β = G II /G T , where G T = G I + G II and G I and G II are energy release rates corresponding to the opening mode and shear mode, respectively. This ratio changes with the adjustment of the loading arm, and the specific relation is expressed by Equation (67) [84].
The P ∼ ∆ relation, as seen in Figure 21, is only decided by the DCB model component [81], and it can be specifically given by In the initial loading with a stationary crack a 0 , the analytical solution of the line OA (stage-I)23, can be obtained by replacing a by a 0 in Equation (68).
When the crack occurs and proceeds, but with a crack length less than L, the length of crack, a, can be determined by the following equation: Then, the P ∼ ∆ relation of the line BCD (stage-II) can be obtained by substituting a into Equation (68).
When the crack propagates beyond the mid-span, i.e., a > L, the crack length can be determined by solving the following equation [57]: Thus, the line EF (stage-III) can be obtained by substituting a into Equation (68). The experimental set-up of an MMB model is schematically illustrated in Figure 21. The model is loaded by a lever with two loading points, where the end loading is in a tied condition and the middle one is a contact condition. Since the lever is considered as a rigid body, the load P and its displacement can be related to the loads and displacements at the left end and at the middle points. Their specific relations are given in Reference [60]. This strategy was also adopted in the present numerical model, where the lever was not simulated for simplicity (see a deformed configuration of the numerical model in Figure 22). The consequence of this numerical set-up is that loads at the left end and at the middle have to be tied with the model. Therefore, the loading condition at the middle is different from the contact condition used in the experimental set-up. This treatment is reasonable when the loading arm is not too short.   Three different cases with different mixed-mode ratios, i.e., β = 20%, 50%, and 80%, were modeled, and different lengths of loading arm, i.e., C = 32.590 mm, 13.226 mm, and 8.443 mm, were set correspondingly. The computed results (with assuming G c = G IIc ) from the proposed model are shown in Figures 23-25, and they are compared with the analytical solution mentioned above. When β = 20% and 50%, the computed results are in good agreement with the analytical solution in all three stages (a = a 0 , a < L, and a > L). The computed results with β = 80%, as shown in Figure 25, overall agree with the analytical solution, particularly for the first two stages. However, they gradually deviate from the analytical results in stage-III with the increase of deflection. This discrepancy is caused by the approximate treatment of the loading condition at the middle, being tied instead of in contact. Specifically, when the loading arm C becomes smaller, the deflection at the mid-span tends to be larger for a certain deflection ∆ at the left end. Therefore, the middle loading point has a bigger trend moving toward the left end because of the rotational effect between the left-end loading and the middle-point loading. The tied condition adopted in this research apparently resists this motion, but the rolling contact condition usually adopted in the experiment does not. It should also be mentioned that the trend of this left motion is pretty small in the models with β = 20% and 50%; thus, the tied loading condition at the middle is effective, and their performance in stage-III is close.
In each case with different mixed-mode ratio, the influence of cohesive strength on the load-deflection relation was investigated. As observed from Figures 23-25, the influence on the fracture behavior (at lines BCD and EF) is overall insignificant, but computed results at the initial loading stage approach the analytical solution with a greater cohesive strength. This is because material tends to be more brittle with a higher cohesive strength, which better matches the assumption of the LEFM-based analytical solution.
Materials 2020, 13, x FOR PEER REVIEW 23 of 35 the mid-span tends to be larger for a certain deflection ∆ at the left end. Therefore, the middle loading point has a bigger trend moving toward the left end because of the rotational effect between the left-end loading and the middle-point loading. The tied condition adopted in this research apparently resists this motion, but the rolling contact condition usually adopted in the experiment does not. It should also be mentioned that the trend of this left motion is pretty small in the models with β 20% and 50%; thus, the tied loading condition at the middle is effective, and their performance in stage-III is close.
In each case with different mixed-mode ratio, the influence of cohesive strength on the loaddeflection relation was investigated. As observed from Figures 23-25, the influence on the fracture behavior (at lines BCD and EF) is overall insignificant, but computed results at the initial loading stage approach the analytical solution with a greater cohesive strength. This is because material tends to be more brittle with a higher cohesive strength, which better matches the assumption of the LEFM-based analytical solution.  . . The effect of different fracture energies was also examined with β 50%. It can be seen from Figure 26 that all three cases with different match their analytical results. Furthermore, varying fracture energy has an impact on crack propagation stages but not on the initial loading stage. This is because crack length a is different (see Equations (69) and (70)) as a result of varying .

Application to the Impact Fracture of a Notched Concrete Beam
Concrete is a quasi-brittle material, and concrete beams are often used to investigate the mechanical performance and the fracture behaviors of the material when subjected to dynamic loading conditions. A significant number of studies on this were carried out experimentally [85,86] and numerically using the FEM associated with the CZM [87] and with continuum damage mechanics [88]. Here, the dynamic fracture behavior of a plain concrete beam, which has an initial notch at the mid-span as shown in Figure 27, was studied using the proposed cohesive crack model. The effect of different fracture energies was also examined with β = 50%. It can be seen from Figure 26 that all three cases with different G c match their analytical results. Furthermore, varying fracture energy has an impact on crack propagation stages but not on the initial loading stage. This is because crack length a is different (see Equations (69) and (70)) as a result of varying G c . The effect of different fracture energies was also examined with β 50%. It can be seen from Figure 26 that all three cases with different match their analytical results. Furthermore, varying fracture energy has an impact on crack propagation stages but not on the initial loading stage. This is because crack length a is different (see Equations (69) and (70)) as a result of varying .

Application to the Impact Fracture of a Notched Concrete Beam
Concrete is a quasi-brittle material, and concrete beams are often used to investigate the mechanical performance and the fracture behaviors of the material when subjected to dynamic loading conditions. A significant number of studies on this were carried out experimentally [85,86] and numerically using the FEM associated with the CZM [87] and with continuum damage mechanics [88]. Here, the dynamic fracture behavior of a plain concrete beam, which has an initial notch at the mid-span as shown in Figure 27, was studied using the proposed cohesive crack model.

Application to the Impact Fracture of a Notched Concrete Beam
Concrete is a quasi-brittle material, and concrete beams are often used to investigate the mechanical performance and the fracture behaviors of the material when subjected to dynamic loading conditions. A significant number of studies on this were carried out experimentally [85,86] and numerically using the FEM associated with the CZM [87] and with continuum damage mechanics [88]. Here, the dynamic fracture behavior of a plain concrete beam, which has an initial notch at the mid-span as shown in Figure 27, was studied using the proposed cohesive crack model. The same model was also studied experimentally by Zhang et al. [86] and studied numerically by Bede et al. [88] using the FEM associated with continuum damage mechanics.
Materials 2020, 13, x FOR PEER REVIEW 25 of 35 The same model was also studied experimentally by Zhang et al. [86] and studied numerically by Bede et al.
[88] using the FEM associated with continuum damage mechanics. The geometry and loading condition are depicted in Figure 27. The notched concrete beam was supported at its two ends and subjected to an impact at the mid-span. The length, width, and depth of the concrete beam were 400 mm, 10 mm, and 100 mm, respectively. Note that the width of the original geometric model in References [86,88] was 100 mm. Since this model is essentially a two-dimensional problem because there is no need to distinguish the difference between the front side and the back side, a 10% proportion was used instead to save computation time. An initial notch of height of 50 mm was located at the mid-span in the bottom of the beam. Two supporters were located 50 mm away from their respective edges of the beam. The impactor was a rigid cylinder featuring a diameter of 8 mm and a mass of 12.06 kg with an initial velocity of 1.76 m/s. The concrete beam had a density of 2368 kg/m 3 , Young's modulus of 31 GPa, Poisson's ratio of 0.18, tensile strength of 5.4 MPa, and critical energy release rate of 141 N/m. The supports were regarded as elastic material with a density of 7900 kg/m 3 , Young's modulus of 210 GPa, and Poisson's ratio of 0.20. The Young's modulus and Poisson's ratio of the impactor were 210 GPa and 0, respectively.
The combined finite-discrete element method was used to model this impact problem. The concrete domain with a range of 64 mm just below the impactor was discretized by DEs to better characterize the fracture behavior, and FEs were used to discretize the two sides of the beam to save computation time. Different element patterns are shown in Figure 28a. The radius of DE particles was 0.5 mm, and there were a total of 64,000 DE particles and 6048 FE hexahedrons. The time step was automatically calculated to keep computation stable, as required by the adopted central difference method. The geometry and loading condition are depicted in Figure 27. The notched concrete beam was supported at its two ends and subjected to an impact at the mid-span. The length, width, and depth of the concrete beam were 400 mm, 10 mm, and 100 mm, respectively. Note that the width of the original geometric model in References [86,88] was 100 mm. Since this model is essentially a two-dimensional problem because there is no need to distinguish the difference between the front side and the back side, a 10% proportion was used instead to save computation time. An initial notch of height of 50 mm was located at the mid-span in the bottom of the beam. Two supporters were located 50 mm away from their respective edges of the beam. The impactor was a rigid cylinder featuring a diameter of 8 mm and a mass of 12.06 kg with an initial velocity of 1.76 m/s. The concrete beam had a density of 2368 kg/m 3 , Young's modulus of 31 GPa, Poisson's ratio of 0.18, tensile strength of 5.4 MPa, and critical energy release rate of 141 N/m. The supports were regarded as elastic material with a density of 7900 kg/m 3 , Young's modulus of 210 GPa, and Poisson's ratio of 0.20. The Young's modulus and Poisson's ratio of the impactor were 210 GPa and 0, respectively.
The combined finite-discrete element method was used to model this impact problem. The concrete domain with a range of 64 mm just below the impactor was discretized by DEs to better characterize the fracture behavior, and FEs were used to discretize the two sides of the beam to save computation time. Different element patterns are shown in Figure 28a. The radius of DE particles was 0.5 mm, and there were a total of 64,000 DE particles and 6048 FE hexahedrons. The time step was automatically calculated to keep computation stable, as required by the adopted central difference method.  As the domain was composed of separate regions, the interaction between different regions was enforced by the coupling approach, as detailed in Reference [40]. The impactor and the two supporters were discretized by hexahedron FEs. The contact between the FE regions and the supports was of FE/FE contact type, while the contact between the impactor and the DE region was of FE/DE contact type. The possible coupling and contact interfaces are shown in Figure 28a. The DE region of the concrete beam initially employed the connective model to represent a continuum. Once this region was under certain deformation and fracture occurred, some DE particles changed their connecting state from connection to cohesion, and they could even be in contact. This model transition was simulated using the proposed cohesive crack model. All above-mentioned contact problems were addressed using the proposed unified contact algorithm, as detailed in Reference [40].
Upon the impact from the cylinder, damage appeared on the top surface of the DE region, as shown in Figure 28b. Note that the damage parameter was defined in Equation (22). This indicates that some DE particles changed their connecting state from the connective model to the cohesive models, but no full decohesion was observed before the damage initiated from the pre-set notch. Note that the occurrence of this damage was just because of the local failure due to an instantaneous impact.
Typical fracture patterns and damage evolutions in the beam at different times are displayed by a series of snapshots in Figure 29. When the contact between the impactor and the DE region persisted, damage initiated from the tip of the pre-set notch, as shown in Figure 29a, and it propagated upward along a straight direction. This crack propagation behavior is consistent with the experimental observation [86] and the numerical simulation using the FEM [88]. This crack pattern was a mode-I fracture at the macroscale. Furthermore, it can be observed that there was a conical damage pattern emerging at the impact point. This damage pattern initiated upon the impact and gradually propagated in a cone form when the contact between the impactor and DE region persisted, and it stopped until the release of contact. The computed impact force versus time is compared to the experimental results [86] and the numerical results from the FEM (based on continuum damage mechanics) [88] in Figure 30, and the peak force and curve shape were the two most important aspects for the evaluation of the proposed model. It should be mentioned that the penalty factor 0.5 was used, calibrated as seen in Reference [40]. As can be seen in Figure 30, the numerical results agree with the experimental results, although the peak value was slightly larger. Compared with the numerical results using the FEM [88], based on continuum damage mechanics, the results from the proposed model are in better agreement with the experimental results [86]. Note that, due to the difficulty in exact experimental measurement, there were some oscillations appearing in the experimental results in the post-impact stage (after about 0.35 ms). However, these oscillations gradually diminished and tended to converge to the numerical results obtained by the proposed model. A numerical case with an initial impact velocity of 0.881 m/s was also modeled, and the computed impact force-time relationship from the proposed model was compared with that from the experiment [86] and the numerical analysis using the FEM associated with continuum damage mechanics [88]. It can be found from Figure 31 that the current numerical results agree much better with the experimental results [86] than those obtained from the FEM [88].  To show the mixed-mode fracture behavior at macroscale, the initial notch was pre-set to 28 mm left of the mid-span. The initial velocity for this numerical case was 1.76 m/s. A series of snapshots of computed fracture pattern at different times are shown in Figure 32. Damage initiated from the tip of the pre-set notch and initially propagated with an inclined angle. This indicates that the mixed-mode behavior dominated at this stage. Then, cracks moved vertically, indicating that mode-I dominated. After a few transitions of moving direction, cracks propagated to the impact point.
The crack propagation behavior which initiated from the notch tip to the area around the loading point agrees with the physical observations [85,89].

Conclusions
A particle-based cohesive crack model was developed for the connective DE model to model the fracture process of brittle and quasi-brittle materials so as to formulate the material transition from a solid phase to a particulate phase. Because of the particle characteristics of the connective DE model, the cohesive crack model was constructed at inter-particle bonds in the connective stage of the model at microscale. A potential formulation was adopted for the CZM, and a linear softening relation was employed for the traction-separation law upon fracture initiation. The criteria of both the fracture initiation and the fracture propagation were constructed in the displacement space. A damage parameter was used to record damage state and, thus, the fracture process was irreversible. An important feature of this cohesive crack model is that the elastic connective stage of the lattice discrete element model was used as the initial loading stage of the CZM, rather than specifying an artificial loading stiffness as often used by conventional intrinsic CZMs when incorporated into the FEM. This particle-based CZM bridges the microscopic gap between the DE connective model and the DE contact model, and it is suitable to describe the material separation process from solids to particulates.
The proposed model was validated by a number of numerical examples, including standard fracture tests of mode-I, mode-II, and mixed-mode I and II. The close agreement between numerical results and the analytical solution confirmed the effectiveness of the proposed particle-based cohesive crack model in modeling the mixed-mode fracture process. The proposed model was also applied to a notched concrete beam subjected to an impact loading. The fracture processes of both mode-I and mixed-mode I and II were reproduced, and crack propagation was consistent with reality. Furthermore, the impact force obtained from the proposed model was in good agreement with the experimental result and more accurate than the result obtained from the FEM, which employed a continuum damage model.
A limitation of the proposed model is that it can only describe the fracture process of homogeneous and isotropic solids, because the DE connective model can only represent isotropic materials. However, this limitation may not compromise the application of the model in modeling the mixed-mode fracture process of a large range of brittle and quasi-brittle materials, such as glass and ceramics. The DE connective model and the particle-based cohesive crack model could be further extended to anisotropic materials in future work.
Author Contributions: H.C.: Literature search, conceptualization, methods, investigation, algorithm, validation, writing, revising; Y.X.Z.: Study design, data interpretation, investigation, supervision, writing and revising; L.Z.: Data analysis, data interpretation, review and editing; F.X.: Data analysis, data interpretation, review and editing; J.L.: Data analysis, data interpretation, review and editing; W.G.: Data interpretation, data analysis, investigation, review and revising. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Symbol Description
f int : (f n ,f s1 ,f s2 ) Interaction force vector and its components in local coordinates. δ: (δ n , δ s1 , δ s2 ) Relative displacement vector and its components in local coordinates. K: (k n , k s1 , k s2 ) Spring stiffness and its components of the normal and two shear spring stiffness. E, E Young's modulus and equivalent Young's modulus. G Equivalent shear modulus. ν Poisson's ratio. r, r Particle radius and equivalent radius. n, s 1 , s 2 The unit base vectors of the local coordinate system to be expressed in global coordinates. n 1 , n 2 , n 3 Components of n. s 1 1 , s 1 2 , s 1 3 Components of s 1 . φ Transformation matrix from the global frame to the local frame. f int , m int Internal force and associated moment. r The effective radius vector. x i Position vector of particle i. T: (T n ,T s1 ,T s2 ) Traction vector and its components in local coordinates. A Effective area. δ n0 , δ s0 Critical values of the relative displacement in the normal and shear directions at the elastic limit. T n0 , T s0 Material strengths in the normal and shear directions. σ c Material tensile strength. λ Non-dimensional scalar of the effective displacement jump.