Crack nucleation in brittle and quasi-brittle materials: A peridynamic analysis

: Peridynamic (PD) models of bodies without pre-cracks, based on a single fracture parameter (associated with the critical fracture energy), produce different strengths when different horizon sizes are used to simulate crack nucleation under quasi-static conditions. To maintain the same strength and fracture energy under different horizon sizes, extra parameters have to be introduced in the failure model. Bilinear and trilinear bond force-strain relationships have been proposed in the literature for crack propagation in quasi-brittle materials. In this paper we study crack nucleation in a plate with a hole under quasi-static loading using bilinear and trilinear PD models. We provide analytical formulas to calibrate the models to measurable material properties. We show convergence for both strength and fracture toughness. The bilinear PD constitutive model works well for both brittle (e.g. ceramics) and quasi-brittle (e.g. concrete) systems, while the trilinear version is more suited for quasi-brittle fracture behavior. We also find that for quasi-brittle fracture, a model that accounts, stochastically, for the presence of small-scale pores/defects performs better than a homogenized model. A wedge-splitting test in concrete and crack nucleation in a quasi-isotropic composite plate with a circular hole are used to demonstrate the model’s performance. In contrast with other models, the current formulation does not depend on the sample geometry.


Introduction
Brittle and quasi-brittle crack propagation under static or dynamic loading have been successfully analyzed using peridynamics (PD) (e.g.see [1][2][3][4][5][6]). PD models converge in terms of crack path and strength (maximum load before failure, in the load-displacement curve) displayed as the nonlocal region (horizon) size goes to zero, when pre-cracks are present.In certain problems, where material or physical length-scale effects are significant in the material's mechanical behavior, an upper bound for the horizon size can be determined based on experimental results (see [7,8]).As we explain in Section 3, existing PD models based on a single fracture parameter (calibrated to match a material's critical fracture energy), lead to strength values that do not converge as the horizon decreases, for problems without pre-cracks or other defects.The reason for this behavior is the absence of a parameter linked to crack nucleation.The prototype microelastic brittle (PMB) model (see [9]) contains only the critical bond strain that is related to the energy release rate for a growing crack, not nucleation of a crack.To address this issue, one can introduce strength-dependent parameters into the peridynamic bond behavior, leading to, for example, bilinear or trilinear bond force-strain relationships.
Recent works based on such bilinear/trilinear behavior of PD bonds focused on aspects of crack nucleation from a material stability point of view and crack propagation [10][11][12][13][14]. Reference [10] presented a bilinear law for PD bond force-strain curve, alongside the original linear bond force-strain curve presented in [15] for the PMB model.In [10], a condition associated with material stability was suggested for spontaneous emergence of a discontinuity ("nucleation" of a crack) in a peridynamic body, and a numerical example of a plate with a centered hole was tested.A study on strength as the horizon size decreases was not provided.
A bilinear model was also employed in [11,12] in PD modeling of static crack growth from a pre-crack in brittle and quasi-brittle materials, such as concrete.The focus was only on crack propagation from existing pre-cracks, not on crack nucleation or the issue of convergence for strength as the horizon decreases.We also note that in [10][11][12], the extra parameter in the bilinear model is determined via curve fitting to an experimental force-displacement curve, with a goal of matching the location where the load-displacement curve departs from linearity as the moment for crack nucleation.In [13], the bilinear model is used to study quasi-brittle behavior, and the model parameters are calibrated to match the experimental force-displacement curves.An alternative approach to computing the model parameters by curve-fitting experimental data, as pursued in [10][11][12][13], is to directly connect these to some measurable quantities (e.g.ultimate strength).In [14], the extra parameter in the bilinear model is set to the ultimate tension strain, i.e.    ⁄ , in which   is the ultimate strength, and E is the elastic modulus.The other parameter in the bilinear model is then evaluated from the critical fracture energy (similar to the calculation of the critical bond strength in the PMB model in [9]).A PD model with a trilinear bond force-strain relationship for quasi-brittle behavior was proposed in [14] with the softening part of the bond force-strain curve being calibrated to the corresponding experimental softening curve.Note that PD models with bilinear force-bond strain behavior have also been used to simulate elasto-plastictype behavior ( [16,17]).In these references, the focus was not on crack nucleation.
Two other approaches to crack nucleation in PD were recently presented in [18,19].In [18], a refinement overlay technique is proposed to decouple the model strength from the grid resolution, so that very small horizon size (with a constant ratio of horizon size to grid resolution) is no longer required to match the actual strength.The special treatment is only needed for PD bonds that are close to reaching their breaking strain, but the extra cost may be significant, especially for problems with many cracks.In [19], an approach based on finite fracture mechanics (FFM) concepts was proposed.FFM has been used before to study crack nucleation in a plate with a circular hole in [20], and for several other simple geometrical settings, such as: elliptical hole [21], V-notch [22,23], straight interfaces in adhesive joints [23,24] and laminated composites [22,23].It appears that an FFM formulation for the general case (e.g.interacting nucleating cracks) is not available.Moreover, FFM approaches are limited to modeling only brittle failure since the basic assumption is that failure happens only when both stress and fracture energy conditions simultaneously reach their critical values.PD models based on bilinear bond-level force-strain relationship cover brittle and quasi-brittle behavior, and the bilinear and trilinear PD models are independent of the geometry of the sample.
The paper is organized as follows: in Section 2 we give a brief review of peridynamics; and then present the numerical discretization used; in Section 3 we study the convergence (in terms of the horizon size going zero) for the global load-displacement curve in a crack nucleation problem; in Section 4 we first derive, based on material properties, parameters for a bilinear bond force-strain model, and then apply the model to study its convergence behavior; here we use the fully homogenized peridynamic (FH-PD) model and the intermediately (or partially) homogenized peridynamic (IH-PD) model that takes into account microstructural defects; in Section 4, we describe a trilinear bond force-strain model and apply it to study crack nucleation; validation against experimental data and comparisons with other PD models from the literature are provided in Section 5; conclusions are given in Section 6.

Original peridynamic model for the PMB material
In this part we briefly review the original PD theory and discuss its numerical discretization using the one-point Gaussian integration, also called the "meshfree" discretization approach.

Brief review of peridynamics
Peridynamics was introduced as a nonlocal form of continuum mechanics by Silling in 2000 [15] for modeling fracture.Since then, it has been extended to a variety of other problems in which domain changes/discontinuities are part of the problem [5,6,8,[25][26][27][28][29][30][31][32][33][34][35].In this theory, each material point is connected through peridynamic bonds to other points within a certain neighborhood region called "the horizon".The peridynamic bonds transfer forces between points (or mass or heat, as in [27,30]) and their failure defines damage at each point.In peridynamics, one replaces the equation of motion by an integro-differential equation in which spatial derivatives are eliminated.This allows peridynamics to avoid the mathematical difficulties and inconsistencies present in the classical theory when cracks, for example, develop in the domain.The PD equations for quasistatic problems are: ∫ ( ̂− , ( ̂, ) − (, )) ⅆ  Ĥ + (, ) = 0 (1) where f is the pairwise force function in the peridynamic bond that connects point  ̂ to x, u is the displacement vector field, and b(x, t) is the body force.The integral is defined over a region Hx called the "horizon region", or simply the "horizon".The horizon is the compact-supported domain of the pairwise force function around a point x (see Fig. 1).The horizon region is taken here to be a circular disk (in 2D or sphere in 3D) of radius δ.We refer to δ also as the "horizon", and, from the context, there should be no confusion whether we refer to the region or its radius.A micro-elastic material is defined as one for which the pairwise force derives from a potential: where  =  ̂−  is the relative position in the reference configuration and  = ( ̂, ) − (, ) is the relative displacement between  and  ̂.A micropotential that leads to a linear microelastic material is given by: where  = ‖‖ , and  = ‖+‖−‖‖ ‖‖ is the relative elongation of a bond, or bond strain.The function () is called the micro-modulus and has the meaning of bond's elastic stiffness.The pairwise force corresponding to the micropotential given above has the following form: where the function  is a history-dependent binary function and  0 is the critical bond strain, calibrated to the material's fracture energy (see below). has the same units as the micro-modulus function () (see Eqs. ( 5) and( 6)).This bond force-bond strain relationship is the constitutive model for a prototype microelastic brittle (PMB) material in a PD formulation (see [9]).
Several choices for the particular form of the micromodulus function () can be used (see [1,37]), but at the macro-scale the differences between these options should be minimal ( [1]).Here, to ease the analytical calculations in Section 4, we employ the "constant" micromodulus function, in which the stiffness of bond (,  ̂) does not depend on its length.The constant micromodulus function for 3D, and 2D plane stress/plane strain are obtained from enforcing a match between the strain energy in PD and the classical strain energy, under a homogeneous deformation (see [1,9]): where E is Young's modulus of the material.We note that with the bond-based version of PD, Poisson ratio ν is fixed to 1/3 in 2D plane stress and to 1/4 in 2D plane strain and 3D conditions.The material model in Eq. ( 4) is equivalent to the kernel function using  = 1 for the peridynamic kernel ( ̂, )/| ̂− |  (see [36]).No significant differences on crack patterns were observed between models with  = 1 and  = 2 [8].In this work, to simplify the analytical calculations in Section 4, we use the model with  = 1.
Failure is introduced in peridynamics by considering that the peridynamic bonds break irreversibly and no longer sustain a force [9,15] when they are deformed beyond a critical value, called the critical relative elongation or critical bond strain, s0, computed based on the material's fracture energy G0.In this study, once a peridynamic bond breaks, it remains broken [27,28].The connection between s0 and G0, for 2D ( [1]) and 3D ( [9]) cases, is: Substituting the constant micro-modulus functions from Eq. ( 6) into Eq.( 7), one obtains s0 as: A damage index quantity is then defined as: In the discrete version, the damage index is computed as the ratio between the number of broken (or failed) bonds (  ) and the total number of bonds (N) originally associated with a node at time (or load step, for quasi-static problems) t (d(x, t) =

𝑁
).When all bonds connected to point x are broken, d(x, t) = 1 and point x becomes a free point.

Numerical discretization
In principle, Eq. ( 1) can be discretized using the finite element method (FEM) [16,38], meshfree direct discretization [9], a combination of both [16,39,40], pseudo-spectral methods [41], or any other method suitable for numerically computing the solution to an integro-differential equation (or an integral equation for the static case).Here we use the meshfree discretization, which makes it easiest to handle damage and fracture [42,43].Using a mid-point integration scheme (see [9]), the discretized version of Eq. (1) at a node xi is: +   = 0 (10) where Fam(i) is the family of nodes j with their area (or volume in 3D) covered, fully or partially, by the horizon region of node i,   is the bond length between nodes i and j,   is the relative elongation for the bond connecting nodes i and j, and   is the area (or volume) of node j estimated to be covered by the horizon of node i.Note that node j may not be fully contained within the horizon of node i, so a "partial volume" integration scheme is used to improve the accuracy [42,44,45].Note also that the system of equations obtained from Eq. ( 10) is linear in bond strains (i.e. = ‖+‖−‖‖
For a fixed horizon size, the ratio  = /Δ describes how accurate the numerical quadrature for the integral in Eq. ( 1) will be.We call this ratio "the horizon factor" [37].We recall that in m-convergence we consider the horizon δ fixed and take  → ∞ (see [37]).The numerical PD approximation will converge in this case to the exact nonlocal PD solution for the given δ.In the case of δ-convergence, the horizon δ → 0 while m is fixed or increases with decreasing δ.For δ-convergence, the numerical PD solutions are expected to converge to the classical local solution (as m increases also) (see [37,47,48]).
For domain discretization, both uniform [9] and non-uniform grids [49][50][51][52] are possible.Nonuniform grids (in which node density does not vary significantly over the domain so that the quadrature error is minimized) are better at modeling rounded shapes compared to uniform grids, and they can be easily created from finite element meshes [53].For the non-uniform grids used in this work, the PD nodes are located at the centroid of their nodal areas (see [29] for a discussion on other options).We use ANSYS to create non-uniform meshes to better match the shape of the hole in plate used to study crack nucleation problems: each PD node is the centroid of a finite element and the element area is assigned to be the nodal area of the corresponding PD node.An ANSYS APDL code exports element centroids and areas (see Appendix B in [33]), which is the only information from the finite element discretization that is used in a PD model.Schematic pictures of element centroid for uniform and non-uniform mesh are shown in Fig. 2. In each case, for both uniform and non-uniform grids, we use constant horizon size over the domain.A good practice is to always select the horizon size for a problem, and then decide on the grid to be used.For the non-uniform meshes mentioned above, however, we first create the grids, and the compute the constant horizon size  by selecting a generic Δ value, related to the average finite element size, multiplied by the horizon factor m.This works in our case because the density of nodes in these grids is fairly uniform, and there are regions where the grid is very close to being uniform (see Fig. 4).A comparison of results obtained by uniform and non-uniform grids with those obtained by experiments for crack nucleation problems are provided in Section 5.All simulations performed in this work are quasi-static.The nonlinear (in displacements) system in Eq. ( 10) can be solved via the energy minimization method using the nonlinear conjugate gradient (NCG) [25,33], or by the adaptive dynamic relaxation (ADR) method [14,39,54].A complete discussion of the NCG method is available in, e.g., [55,56], and a brief account is given in Appendix A of reference [33].The micropotential function  for the PMB model in Eq. ( 3) is quadratic (and therefore convex), while for the bilinear and trilinear bond force-strain models it is not.It is well known that minimizing non-convex functions (see Fig. 7) makes convergence to the global minimum more difficult (see [55,57]).Because of this, here we use the ADR method for the bilinear and trilinear PD constitutive models.A discussion of the ADR method is provided in Appendix A. The algorithm we use for simulating quasi-static fracture with PD is shown in detail in Appendix B.

Crack nucleation for the PMB model
The original PD model based on a single fracture parameter relies on obtaining the critical bondstrain from the measured critical fracture energy (see Eq. ( 8)).This means that one can match the energy release rate with a PD model, but not, at the same time, the material strength (related to crack initiation phase).This is not an issue for cases in which there is a pre-crack, because in these cases, the PD-computed strength is controlled by (or coupled with) the fracture toughness, whereas in problems in which there is no pre-crack, the PD computed strength is independent from (not controlled by) the fracture toughness (see [10,[18][19][20]).Because of this, PMB models for problems with pre-cracks do show convergence (see [37], and Section 2.2 for types of convergence in PD models) for both strength and dissipated fracture energy.Appendix C discusses -convergence for the PMB model for a plate with a pre-crack.Note that while microstructural details in a material generate dependencies between a material's strength and its fracture toughness, there are no mathematical theories directly connecting these quantities.Because of this, we will attempt to introduce models that match some given values of strength and toughness.A description of how microstructural arrangements can lead to high toughness without sacrificing the strength and stiffness is given in [58].
In this section, we use the original PD model to study -convergence in problems that do not contain pre-cracks or other defects.Reference [18] showed that for a plate without a hole, the original PD model gives higher and higher strengths as the horizon size decreases.We consider a plate with a circular hole under quasi-static loading.This geometry is selected in order to control the location for crack nucleation, which should happen from the points with maximum stress concentration factor.
Geometric parameters used are shown in Fig. 3.The thickness of the plate is considered to be 0.5 mm and used in the post-processing step to compute the total load on the plate for plotting load-displacement curves.Material properties are chosen as: E = 192 GPa and G0 ≈ 83 kJ/m 2 (Section 9.1 of [54]).Note that we calculated G0 through Eq. ( 8) with the critical bond strain ( 0 =0.02), the horizon size and the elastic modulus given in Section 9.1 of [54].The material properties are hypothetical values that under certain conditions may lead to brittle or quasi-brittle fracture behavior (in section 4.2 we will show how we control the fracture behavior to be brittle and quasi-brittle by involving an extra material property) A vertical displacement-controlled type loading with a constant rate of 0.275µm per step is applied to each end of the plate (see Fig. 3).2D plane stress conditions are considered.A dynamic relaxation method (ADR) solver with a time step of  = 0.01 s is used (see Appendix A).
We use non-uniform grids generated by ANSYS APDL (see Section 2.2) for an eighth of the structure, as shown in Fig. 4 and Fig. 5.We use symmetries to then construct the grid for the entire plate and solve the problem over the entire plate, without using symmetry, in order to observe fracture without forcing symmetry conditions into the model, other than the natural ones present in the geometry itself.This option also leads to creating more uniform element sizes over the sample, which helps avoiding zones with few nodes (and therefore with low quadrature accuracy) that could happen if we were to generate the FE mesh directly for the entire sample.PD nodes are the centroids of corresponding finite elements in the grid and nodal areas are the element areas (see Section 2.2).We consider regions of fictitious nodes with the width of δ at the top and bottom of the plate to apply the displacements that control the loading ( [59,60]).For example, for the horizon factor of m = 5, five rows of fictitious nodes in the top and bottom of the plate are considered.These rows are on uniform grids, since the non-uniform mesh becomes uniform far enough from the hole.Over these fictitious nodes we enforce the constant values of the applied displacements, with the goal of mimicking the imposition of the corresponding local boundary conditions.Note that imposing local-type boundary conditions may be desired/needed because, in reality, such conditions can only be measured at the surfaces of a body, not through a finite layer near the surface, which what nonlocal boundary conditions would entail.Different implementations of the fictitious nodes method (FNM) are possible to impose local-type boundary conditions in PD models, see [48].
A comparison between displacements obtained by the PMB PD model (solved with ADR) with those from FEA, for the elasticity problem (when no failure is allowed), is shown in Section 9.1 of [54].Fig. 4(a) shows the force-displacement curves from the current PD simulations of failure with the PMB model (see Eq. ( 4)) and different horizon factors m.The non-uniform grids for several horizon factors are shown in Fig. 4(b-d).These m-convergence results allow us to use m = 5 for the remaining simulations in this section and in Section 4. Note that we compute the total load at the top end of the plate (i.e. at the cross-section separating the top of the plate from the fictitious nodes above it), by summing up the vertical components of nodal forces.
In Fig. 5 we show the results with the PMB model for different horizon sizes.As the horizon decreases, we reach higher and higher strengths (load at which failure initiates).The reason for this behavior is the absence of an independent parameter linked to crack nucleation.The PMB model contains only the critical bond strain, which is calibrated to match the fracture energy release rate.To address this issue, we will introduce strength-dependent parameters into the peridynamic bond behavior, leading to bilinear or trilinear bond force-strain relationships (see Section 4).Some alternatives have appeared recently [18,19], but they come with some limitations/deficiencies (see discussion in Section 1).
To analyze the influence spatial discretization has on strength values for our problem, we also use uniform grids.In Fig. 6 we show the increase of strength values with a decreasing horizon size and we compare the results obtained with uniform and non-uniform grids.As expected, strength values calculated by using uniform grids are lower than corresponding ones from non-uniform grids because a uniform grid around a circular hole shape will be staggered and thus lead to artificial stress concentration points (see also [49,53]).

PD models with softening
One way to reach convergence for strength values under decreasing horizon values is to introduce an independent parameter related to ultimate strength.Instead of the sudden bond-failure in the PMB model, one can use a bond force-bond strain constitutive model with softening.This allows for additional parameters besides the critical bond failure one (see Fig. 7 for different options for PD bond force-strain curve).One of the simplest versions for such an approach is a bilinear PD model (see Fig. 7b).In this section we show how to calibrate this to strength values and study crack nucleation for the plate with a circular hole.A trilinear model is also discussed in this section.
The bilinear and trilinear constitutive models for PD bonds allow for the simulation of both nucleation of a crack (primarily associated with the peak of the PD bond force-strain curve) and crack propagation through full failure (primarily associated with the bond breakage point in the PD bond force-strain curve).

4.1.
Bilinear PD softening model Methods using a bilinear bond force-bond strain model have appeared in, for example, [10][11][12][13][14], where some approximate methods were used to obtain the parameters involved.Here, we pursue an analytical approach to connect the model parameters to the ultimate strength (associated with crack nucleation) and fracture energy (associated with crack propagation).
Similar to the approach used for the calibration of the PMB model (see Eq. ( 7)), we enforce a match between the material's fracture energy and the energy required in the PD model to break a sample in two (see Fig. 8), as follow: This leads to: For the additional equation needed here, we use the ultimate strength value (a material property available from experimental data).The ultimate strength value is controlled by small pores, defects, microcracks, and it is reasonable to assume that most of these grow and lead to either softening or sudden failure because of isotropic expansion (zero deviatoric components) local conditions.We connect the maximum bond force cs1 (see Fig. 7) to the ultimate strength by assuming that PD stress normal to the plane along which a crack initiates is equal to the ultimate tensile strength (see [16]): For each point A along 0 ≤ z ≤δ, the force for bonds connecting A to each point B (assumed all to be at their maximum value, for crack initiation) is integrated by Eq. ( 13) (redrawn, with modifications, from [1] and [9]).
From Eq. ( 13), one obtains: By substituting the constant micromodulus functions from Eq. ( 6) into Eq.( 14), the bilinear parameter is obtained as: and by substituting s1 from Eq. ( 15) into Eq.( 12), the other bilinear parameter is obtained as: For 2D plane stress and plane strain conditions, the same sc is obtained analytically.Note that similar formulae to Eq. ( 15) are found in 3D ( [16]) and 2D plane strain ( [17]) but in the different context of elasto-plastic-type behavior.
Note that parameter s1 is independent of horizon size, but sc is inversely proportional with the horizon size.For the same PD material model, using two different horizon sizes δ1 > δ2, the bilinear model will look as shown in Fig. 9.We note that the condition sc > s1 is required to construct the bilinear model with a softening part.This condition imposes an upper limit on δ, assuming some given material properties, as follows (see Eqs. ( 15) and ( 16)): One can also use Eq. ( 17) to determine the range of material properties which, for a particular horizon size, can be analyzed with the bilinear constitutive model.

Crack nucleation for the bilinear PD model
In this section we utilize the bilinear model (see Fig. 7(b)) to study nucleation of a brittle or quasibrittle crack in the plate with hole under tensile loading (see Fig. 3).The PD models discussed so far are homogeneous.Real materials are heterogeneous, and a closer representation of them is given by the IH-PD model ( [5,6,33,60]), which is a stochastic PD model able to take into account, for example, the presence of pores/defects/phases into the elastic and failure behavior of materials ( [5,6]) without having to include the detailed representation of such pores/defects/phases.Our previous studies showed such models are able to capture initiation of damage and failure in materials with critical pre-notches, matching experimental observations (see [5,6]).In this section we will also test the IH-PD model for crack initiation from a hole and observe if there are advantages compared with homogenous PD formulations.
In our example, with the IH model, to account for the presence of defects, we assume a porosity of 1% for the brittle and quasi-brittle cases.Many ceramic materials, which are nominally brittle materials, have a porosity around 1% (see Table 5.5 in [61]).Crystalline rock, such as granite, a quasi-brittle material, can have a very low porosity (<1%) [62].For a large variety of low-porosity brittle and quasi-brittle materials, the critical porosity   (the porosity beyond which the material can exist only as a suspension [63]), can be higher than 0.9 (or 90%) [63][64][65].Here we consider   ≈ 1 as input in the simulations with the IH model.
To determine the material properties needed in the IH-PD model of poroelastic materials, the following relationship can be used to find  sp (see [5]): where  sp is a particular property of the solid phase (sp) constituent in the porous medium (e.g., elastic modulus, fracture energy, and ultimate strength),  eff is the effective material property for the porous medium,  is material's porosity, and   is critical porosity.The last three quantities are taken from experimental measurements.Note that, material properties to be used in the bond force-strain models (see Eqs. ( 6), ( 8), (15), and ( 16)) for the IH-PD model are the properties of the solid phase constituent material of the porous medium,  sp ,  sp ,  ,sp , not the effective properties that were used in the FH-PD model above.For our example, the solid phase Young's modulus and fracture energy for    = 1% and using the values given in Section 3, are obtained as (see Eq. ( 18))  sp = 195.9GPa, and  sp =84.6 kJ/m 2 .
In this section we consider two materials with different strengths but the same Young's moduli and fracture energies (with values given in Section 3).In Table 1 we give the corresponding  1 values.If we consider   = 0.576 GPa, we find that the bilinear bond model leads to a global quasi-brittle fracture behavior (see Fig. 10(a)), while with   = 1.44 GPa, we get a macroscale brittle fracture behavior (see Fig. 10(b)).These values for   are arbitrarily chosen so that they lead to a global fracture behavior that appears to be brittle (sudden failure) or quasibrittle (failure with a softening curve).The horizon sizes and the material properties given in Table 1 satisfy Eq. ( 17), which is required for the bilinear model.Fig. 10 and Fig. 11 show the PD simulation load-displacement results for the crack nucleation problem from the FH-PD and IH-PD models, respectively, when the bilinear bond constitutive model is used.Fig. 10 (a) and (b) show the bond-level constitutive models responsible for the global responses shown in Fig. 10 (c) and (d), respectively.Softening behavior in the constitutive law at the bond level can lead to an abrupt/sudden global failure (brittle behavior) under certain selections of bond parameters (see Table 1 and Fig. 10 (b)).To transition to a quasibrittle global response, the softening part in the bond-level model needs to be extended significantly compared to one that leads to brittle failure (see Fig. 10 (b)).These observations correlate well with the physical understanding that in brittle fracture, the extent of the fracture process zone (FPZ) is, in general, smaller than the extent of the FPZ for quasi-brittle fracture.Note that the quantity plotted in Fig. 10 (a) and (b), fdA, is the integral kernel in Eq. ( 1) (dA is the nodal area).
Fig. 12 shows the variation of the ultimate load with respect to the horizon size.These results show that the PD model with a bilinear constitutive model for bonds leads to δ-convergence for crack nucleation, for both brittle and quasi-brittle behavior.Fig. 13 compares the FH-PD and IH-PD results in terms of the relative difference between a base value for the ultimate load (obtained with the smallest horizon size used in this work  * =0.375 mm; see below) and the corresponding value computed for a certain horizon size, .This relative difference, for each model in part, is computed by: |P u () − P u ( * )| P u ( * ) where P u ( * ) is the computed ultimate load obtained with the smallest horizon size (here,  * = 0.375 mm).As seen in Fig. 13, for both brittle and quasi-brittle behavior, the relative difference in Eq. ( 19) decreases with the horizon size shrinking.For brittle-type fracture behavior, there is little difference between the two models.However, for a quasi-brittle behavior, the IH-PD model, even with a large horizon size, computes an ultimate load value that is much closer to the "converged" value computed with a significantly smaller horizon size.This means that for quasi-brittle fracture, the IH-PD will be more efficient to use than a corresponding FH-PD model because we can reach a reasonable ultimate load using a much coarser horizon size (and coarser corresponding discretization grid).This issue is discussed in detail in [5,6] in the context of materials with high and very high porosities.
In Appendix D, we also compare the computed ultimate load with an empirical estimate for the ultimate load often used in engineering design.The formula for that estimate was based on the photoelastic method, which is supposed to serve as a conservative estimate for a variety of materials.

Trilinear PD model
In this section we review the trilinear PD model, introduced in [14], and then use it to model crack nucleation in the plate with a hole described in Section 3. We then compare the results with those obtained in the previous section by the bilinear model.In Section 5.1, we will compare the performance of these models in an example of quasi-static fracture in concrete, also studied in [14].
Unlike brittle materials, quasi-brittle materials have a complex strain-softening behavior due to the existence of a larger fracture process zone (FPZ) ( [66,67]).The global stress-strain curve for quasi-brittle materials is generally not easily obtained from direct tension tests because it is strongly influenced by the specimen's dimensions, which is a manifestation of size effects in quasi-brittle fracture.Instead, a softening constitutive curve has been used to investigate failure of quasi-brittle materials, especially in concrete ( [14,66,67]).The softening constitutive curve gives the relationship between the "cohesive" stress (stress normal to the propagation direction at the crack tip) and the critical crack opening displacement (COD), which can be regarded as a material property.Reference [14] used the softening constitutive curve to introduce a trilinear constitutive PD model, in which the force-strain curve is controlled by three parameters.
For the trilinear model in the PD constitutive model, parameter s 1 is the same as that in the bilinear model (see Fig. 7(b) and (c), and Eq. ( 15)).By equating  0 from bilinear and trilinear models (corresponding to equating the areas below the bond force-strain curves; see Fig. 7 (b) and (c), and Eq. ( 11)), one obtains the other two parameters in the trilinear model as follow: In order to compute  2 and ĉ, β and  need to be determined.These can be obtained from a bilinear approximation of the experimental softening constitutive curve, as suggested in [14], if such a curve is available from fracture tests.In the absence of an experimental softening curve, reference [14] suggests the following prescribed values for a generic concrete material:  = .In Section 5.1, we will show that using these values for β and  leads to poor performance in reproducing some experimental fracture data.
In the trilinear model, the condition s2 > s1 needs to be satisfied, and according to Eq. ( 20), this leads to the condition sc > s1.This condition is the same with that for the bilinear model which imposes upper limits on δ for given material properties (see Eq. ( 17)).Note that other conditions for the trilinear model, i.e. the conditions  < 1 and  > 1) are satisfied by choosing the abovementioned values for β and  (see Fig. 7(c)).
Note that the constitutive models for PD bonds we introduced here should also work for dynamic problems, since no assumptions were made in the derivations of the mathematical formulas obtained that would restrict their application to static problems.

Crack nucleation for trilinear PD models
Here, we apply the trilinear model for the crack nucleation problem discussed in Section 3. We choose β and  as suggested in [14] (see Section 4.3).Similar to Section 4.2, we here utilize the both FH-and IH-PD models.Fig. 14 shows the FH-and IH-PD results with the trilinear model for a quasi-brittle material with   = 0.576 GPa (see Section 4.2).We note that the trilinear PD model is defined based on the softening curve, and has been used for modeling failure of quasi-brittle materials (e.g.concrete) (see [14]).
From Fig. 14, we observe that using the trilinear PD model, the force displacement curves for smaller and smaller horizons converge throughout the loading sequence.We also see that utilizing IH-PD model alongside with the trilinear model leads to a faster convergence.

Validation against experimental results and comparisons with other PD models
In this part we compare the results obtained by our PD models with softening (bilinear and trilinear) against two experimental tests and against those obtained by other PD models.We first use a wedge-splitting test in concrete (Section 5.1), and then look at crack nucleation in a woven glass/epoxy composite plate with a center-hole (Section 5.2).

Validation for a wedge-splitting test in concrete
We consider the wedge-splitting test of a concrete specimen ( [68]) (see Fig. 15), with material properties given as: E=28.3 GPa,   = 2.27 MPa, and G0=490 N/m.We perform PD simulations under 2D plane strain conditions.Quasi-static loading conditions are mentioned in [68], [69].In our simulations, we apply a horizontal displacement-controlled type loading with a constant rate of 60 nm per loading step, as shown in Fig. 15.We use δ=0.193 m and m=3, as used in [14].
Reference [14] used this experiment to validate results obtained by bilinear and trilinear PD models.A conclusion drawn in [14] was that results with the trilinear PD model matched better the experimental data compared with those from a bilinear PD model (see Fig. 15(a) in [14]).
Fig. 15.Wedge-splitting test specimen (unit: cm) (redrawn from [14,68]) In this section we will test whether the bilinear PD model we introduced here (with analytically-derived parameter, see Section 4.1) leads to results that are comparable, or not, to those from the trilinear model.For the trilinear model, we will employ two different ways of computing its parameters (see Fig. 16(a)).
Fig. 16 shows the comparison of load versus crack mouth opening displacement (CMOD) in the crack propagation problem for the wedge-specimen, utilizing different bilinear and trilinear PD models.We compute the solution with two different bilinear models (see Fig. 16(a)): one in which parameter s1 is determined by curve fitting to an experimental force-displacement curve [10][11][12][13] (or setting it equal to    ⁄ ( [14], see Section 1); a second one in which parameter s1 is obtained from Eq. (15).We also provide the results from two trilinear models: one in which the parameters in the trilinear model (see Eq. ( 20)) are calculated based on an available experimental softening curve; the other in which the parameters are the prescribed values suggested in [14].Reference [68] presented the experimental softening curve for the concrete specimen used in the wedge-splitting test, and then presented a bilinear approximation of the curve.Thus, parameters in the trilinear model are directly evaluated as β=0.29/2.27=0.1277,and γ=1.42/0.25=5.68[14,68].If the experimental curve was not available, we would choose β=1/4 and γ=28/3 suggested in [14] (see Section 4.3).Fig. 16 shows that the second trilinear model leads to poor performance in reproducing experimental data.The bilinear PD model with analytically-derived parameter leads to compromising results with those from the trilinear model suggested in [14].Since for the trilinear model, an experimental softening constitutive curve (as an extra material properties) is needed for an accurate fracture analysis, we suggest using the analytically-derived bilinear model introduced here to study quasibrittle fracture, e.g. in concrete materials.

Validation for a crack nucleation problem
We now compare the results obtained by our bilinear PD bond model (derived analytically; see Section 4.1) with the experimental data in [70] for crack nucleation in center-hole plate specimens of woven glass/epoxy composite.This experimental data was also used to compare with the prediction of the PD model based on finite fracture mechanics (FFM) in [19].
We note that the composite specimens have similar elastic properties in vertical and horizontal directions (E xx =E yy =23.6 GPa) corresponding to the quasi-isotropic materials, and they show brittle fracture behavior [70].
We consider the plate with thickness of 40 mm and with six radii for the centered hole (similar to Fig. 3).In order to compare the results obtained by our PD models with those by PD model proposed in [19], we conduct our PD simulations under the same parameters and conditions taken by [19].We use δ = 0.2 mm and m = 4 and we consider plane stress and displacementcontrolled loading conditions.In [19] a uniform grid is employed.We employ both uniform and non-uniform grids.As shown in Fig. 17, the bilinear PD bond model reproduces the strength in the experiment more accurately than the PD model presented in [19] does.
Since non-uniform grids conform better the region where a curvature exists (see Section 2.2), we use non-uniform grids for this problem.We generate non-uniform grids (see Section 3) for the plate with six radii of the hole.A comparison of strengths computed by bilinear PD bond model utilizing non-uniform and uniform grids with those measured from the experimental data is shown in Fig. 17.Note that the notched strength is the maximum strength that the plate with hole can reach, and this is different from the ultimate strength which is a material property.We calculate the relative deviation and maximum error of our PD simulation results with respect to the experimental data.We calculate the relative deviation by using ( [19]): =1 (21) and the maximum relative error by using: where  Exp (   ⁄ ) is the notched strength for hole radius   , and  PD (   ⁄ ) is the corresponding computed strength by PD. 2 is the plate width.The index  is one of the 6 configurations used in the experiment corresponding to a different hole radius.
Table 2 shows that deviations and maximum errors in notch strength of our PD results relative to the experimental data are significantly smaller than those obtained with the FFM-based PD model in [19].We also see that the results with non-uniform grids deliver the smallest errors relative to the experimental data ( PD and  rel,max ).FFM considers a finite characteristic length, in addition to the natural one given by the discretization size, related to geometrical parameters.In FFM one inserts a crack of that length at the nucleation site, and convert the crack nucleation problem, from computational viewpoint, into a crack propagation problem ( [19][20][21]).This extra length-scale in the FFM model might be the reason for the FFM results deviating from the experimental results in this problem (see Fig. 17).In contrast, our PD models with softening are able to track very well the experimentally-observed notched-strength dependency on sample geometry.

Conclusions
In this paper we showed that peridynamic (PD) models with two or three damage-related parameters lead to convergence, in terms of the nonlocal size, for crack nucleation and growth in brittle and quasi-brittle materials with no pre-cracks.These bilinear and trilinear "constitutive models" for bond force-strain relationships are calibrated to the fracture energy and the ultimate material strength.The original PD models based on a single fracture parameter (associated with the critical fracture energy) also converge if there is a pre-crack in the domain, but produce different strengths when different horizon sizes are used to model crack nucleation under quasistatic conditions in bodies without pre-cracks.
We provided analytical formulas for the calibration of the bilinear and trilinear models to measurable material properties.The bilinear form of the PD constitutive model showed good performance for both brittle (e.g.ceramics) and quasi-brittle (e.g.concrete) systems, while the trilinear version was more suited for quasi-brittle fracture behavior.We also found that for quasibrittle fracture, a PD model that accounts, stochastically, for the presence of small-scale pores/defects performs better than a homogenized model.Using these PD models, we conducted numerical simulations for a wedge-splitting test in concrete and crack nucleation in a woven-type composite plate with a circular hole.In the concrete fracture problem, we showed that even a bilinear PD model can capture well the failure and post-peak softening curve observed experimentally.The results for the plate with a hole match very well the experimentally measured dependence of strength on hole size.In contrast with other methods (e.g.those using Finite Fracture Mechanics approaches in PD), our model is more general, does not depend on the sample geometry, and the results match closer those obtained experimentally.the same load step) until no more bonds break and the physical system reaches equilibrium.If too many bonds break after one solution, numerical instabilities may happen [25].In the examples shown in this paper, the crack growth is stable, and as long as a reasonable number of load steps are used, numerical instabilities are avoided.

Appendix D
In this appendix, we estimate the ultimate load in Fig. 12 by a correlation that has been used for the purpose of engineering design.This correlation is derived through the use of a photoelastic method and is based on the definition of stress concentration factor.This factor is the ratio of maximum stress  max (  at θ = 0 o and at the edge of the hole) to nominal stress  nom (the average stress at the reduced cross-section due to the hole).Remote stress  ∞ is the stress applied on the plate ends.With these, one obtains ( [71]): We enforce  max =   .We use the two values for   chosen for brittle and quasi-brittle materials (see Table 1).Using Eq. (D1) we find the remote stress  ∞ by which crack nucleation is about to happen.The corresponding remote load related to crack nucleation (or the ultimate load shown in Fig. 12) is  ∞ =  ∞ .. , where W is the plate width and b is the plate thickness.The ultimate (remote) loads for the brittle and quasi-brittle materials is obtained as 9.73 kN and 3.89 kN, respectively.These estimated values are a bit smaller than the corresponding values shown in Fig. 12, since Eq.(D2) is mostly used for the engineering designs that seek allowable stress (or load) before failure/yielding being about to happen.For the engineering designs of elasto-plastic materials (such as steels), one may use Eq.(D1) with enforcing  max =   , where   is the yield strength, for the materials not being yielded.

Fig. 1 .
Fig. 1.A point x interacts directly with any point  ̂ in its horizon region Hx , here a disk of radius δ.

Fig. 2 .
Fig. 2. Finite element centroids for (a) uniform and (b) non-uniform meshes are chosen as PD nodes.The element area is used as the corresponding PD nodal area.

Fig. 3 .
Fig. 3. Geometry and loading configuration for the plate with a hole."Δ" denotes the displacement-controlledtype loading.

Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 4. In (a): force-displacement curves from PD simulation results with the PMB model with horizon size δ = 3 mm for different m-factors.Grids shown (used, with symmetry, to create discretizations for the plate in Fig. 3), are for horizon factors (b) m = 3, (c) m = 4, and (d) m = 5.Red circles show the horizon size.

Fig. 7 .
Fig. 7. PD bond micropotentials (top row) and force versus strain (bottom row) for (a) the PMB model, and softening models with (b) a bilinear behavior, and (c) a trilinear approximation.

Fig. 8 .
Fig. 8. Evaluation of ultimate tensile strength,   , for (a) 3D, and (b) 2D cases.For each point A along 0 ≤ z ≤δ, theforce for bonds connecting A to each point B (assumed all to be at their maximum value, for crack initiation) is integrated by Eq. (13) (redrawn, with modifications, from[1] and[9]).

Fig. 12 .
Fig. 12.Comparison of computed ultimate load for different horizon sizes using the PMB and the bilinear model (with FH-PD and IH-PD) for brittle and quasi-brittle behavior.

Fig. 13 .
Fig. 13.Relative difference in computed ultimate load (see Eq. (19)) for FH-PD and IH-PD with the bilinear model, for brittle and quasi-brittle behavior.

Fig. 14 .
Fig. 14.PD bond force (times nodal area) vs. bond strain (a) and global load-displacement curves obtained with the FH-PD (b) and IH-PD (c) using the trilinear constitutive model for quasi-brittle behavior.

Fig. 16 .
Fig. 16.Comparison of bilinear and trilinear PD bond-force (times nodal area) vs. bond strain (a) and the corresponding results obtained with these different models plotted against experimental data from Ref. [68] (b).

Fig. 16
Fig.16also shows that the analytically-derived bilinear model can match the peak in the experimental load-displacement curve, while the approximate bilinear model overshoots the peak.The bilinear PD model with analytically-derived parameter leads to compromising results with those from the trilinear model suggested in[14].Since for the trilinear model, an experimental softening constitutive curve (as an extra material properties) is needed for an accurate fracture analysis, we suggest using the analytically-derived bilinear model introduced here to study quasibrittle fracture, e.g. in concrete materials.

Fig. 17 .
Fig. 17.The comparisons of strengths computed by the bilinear PD model with those by the PD model based on FFM in [19], and with those measured from experiment ([70]).

Fig. B1 .
Fig. B1.Flowchart for the simulation process of PD quasi-static fracture modeling.Here t is the load step.
concentration factor   is obtained experimentally through the use of a photoelastic method for a variety of materials, as follows (see chart 4.1 of[71]):

Table 2 .
Relative deviations and maximum relative errors with respect to experimental data