## Abstract

The recently conceived gap test and its simulation revealed that the fracture energy G_{f} (or K_{c}, J_{cr}) of concrete, plastic-hardening metals, composites, and probably most materials can change by $\xb1100%$, depending on the crack-parallel stresses σ_{xx}, σ_{zz}, and their history. Therefore, one must consider not only a finite length but also a finite width of the fracture process zone, along with its tensorial damage behavior. The data from this test, along with ten other classical tests important for fracture problems (nine on concrete, one on sandstone), are optimally fitted to evaluate the performance of the state-of-art phase-field, peridynamic, and crack band models. Thanks to its realistic boundary and crack-face conditions as well as its tensorial nature, the crack band model, combined with the microplane damage constitutive law in its latest version M7, is found to fit all data well. On the contrary, the phase-field models perform poorly. Peridynamic models (both bond based and state based) perform even worse. The recent correction in the bond-associated deformation gradient helps to improve the predictions in some experiments, but not all. This confirms the previous strictly theoretical critique (JAM 2016), which showed that peridynamics of all kinds suffers from several conceptual faults: (1) It implies a lattice microstructure; (2) its particle–skipping interactions are a fiction; (4) it ignores shear-resisted particle rotations (which are what lends the lattice discrete particle model (LDPM) its superior performance); (3) its representation of the boundaries, especially the crack and fracture process zone faces, is physically unrealistic; and (5) it cannot reproduce the transitional size effect—a quintessential characteristic of quasibrittleness. The misleading practice of “verifying” a model with only one or two simple tests matchable by many different models, or showcasing an ad hoc improvement for one type of test while ignoring misfits of others, is pointed out. In closing, the ubiquity of crack-parallel stresses in practical problems of concrete, shale, fiber composites, plastic-hardening metals, and materials on submicrometer scale is emphasized.

## 1 Introduction

In computational mechanics, enormous attention has recently been paid to two numerical approaches—the phase-field (PF) models [1–9] and the peridynamic (PD) models [10–14]. Yet proper experimental justifications of both approaches by a sufficiently broad range of the experimental results have been missing despite the abundance and diversity of relevant experiments in the literature. A recent review [15] attempted to compare the performances of these models on various experimental data, but the scope of that review and the experiments examined was still far too limited and scrutiny to shallow. Serious discrepancies with experiments and their causes in the basic concepts continued to receive no attention. The present objective is to remedy this situation.^{2}

In the case of peridynamics, grave conceptual deficiencies have been identified in 2016 [16], though only by theoretical analysis. Even though several partial remedies of these deficiencies [17–19] have been attempted, no comparisons of these remedies with the full range of relevant experimental fracture data available in the literature have been made. Therefore, the predictions of PF and PD models are here compared to a broad range of classical laboratory experiments important for practical applications and also to the results of the recently developed gap test [20,21].

The experimental comparisons of PF and PD models that exist in the literature (e.g., Refs. [15,22–24]) involve only a few selected simple types of tests that can be fitted equally well by many different models. Such selective comparisons are insufficient to justify applicability in engineering practice. While peridynamics symposia have, for instance, featured impressive, realistically looking, videos of impact of projectiles, no comparisons have been made with the ample existing data on the penetration depths, wall exit velocities and crater dimensions and shapes for various entry velocities, wall thicknesses, and projectile dimensions, which would have revealed severe disagreements. Moreover, the material model used has not been verified by the available comprehensive laboratory test data on small specimens, which could provide the material constitutive properties. Such data cover the full range of uni-, bi-, and tri-axial experiments, proportional and nonproportional, with stabilized postpeak and rotating principal stress directions, under various confinements, tension-compression transitions, unloading and load cycles, etc. (e.g., Ref. [25]).

Here, we aim to fill this gap, for both PF and PD models. In each case, comparisons are also made with the predictions of the finite element crack band (CB) model with microplane damage vectorial constitutive model M7^{3} [26,27]. Our study focuses mainly on concrete, not only because of the dominance of this material in national economies and its huge CO_{2} footprint but also because the diversity of the available fracture and damage test data on this material is much greater than for any other. Besides, other quasibrittle materials [28,29], as well as polycrystalline metals and polymers at micrometer scale, are known to behave similarly [30] (the similarity includes widely ignored phenomena such as the vertex effect [31,32]).

Model M7 for concrete, which is the damage constitutive model used here [25,31], has been calibrated to fit virtually the full range of material test data that exist for concrete [25], which include about 20 different types of material tests. In the CB model, the size of the finite elements controls the crack band width and is proportional to the material characteristic length *l*_{0} that serves as a localization limiter. This length is best determined by fitting the test data on the size effect, while the determination of *l*_{0} from stabilized postpeak softening may be ambiguous [33]. A very rough estimate of *l*_{0} is one to three times the inhomogeneity size of the material.

## 2 Overview of the Gap Test

The recent discovery [20,21] of the simple gap test revealed that the fracture energy *G*_{f} and the effective width *c*_{f} of the fracture process zone (FPZ) depend strongly on the crack-parallel normal stress *σ*_{xx} (= *T*) in the direction *x* of crack propagation and, according to a calibrated CB model, also on the crack-parallel stresses *σ*_{zz} and *σ*_{xz}. In this test (Fig. 1(a)), elastic-perfectly plastic loading pads are used to apply onto a standard notched three-point-bend beam a uniform compression along the notch the end supports of the beam are installed with a certain gap. These supports would engage in contact only after a constant plastic plateau has been attained in the pads (Fig. 1(b)). Scaled notched concrete beams of size ratios 1: 2: 4: 8 have been tested to extract the fracture energy. From their peak loads (Fig. 1(c)), one can evaluate the fracture energy *G*_{f} and the effective size *c*_{f} of the FPZ by means of the size effect method [34] (which is a standard RILEM Recommendation T89-FMT [35] and is endorsed by ACI-446). This method can be reduced to linear regression, which yields both *G*_{f} and *c*_{f} (with their standard deviations).

The tests show that the crack-parallel compression *σ*_{xx} (for metals previously called the *T*-stress) has a strong effect on *G*_{f} (Figs. 1(d) and 1(f)). This effect has for a long time gone unnoticed, obviously because all the standard fracture specimens feature a zero or negligible *σ*_{xx}, *σ*_{zz}, and *σ*_{xz}. This effect was also unsuspected because in the linear elastic fracture mechanics (LEFM) and the cohesive crack model (CCM), the crack is assumed to be a line (of zero width at crack front), while a line cut along the direction of a uniform uniaxial stress field produces no stress change. The transverse normal stress *σ*_{zz} has also a large effect, as indicated by the simulations in Ref. [21]. A significant effect might also be expected for *σ*_{xz}.

In computations, one might be tempted to approximate the curves in Figs. 1(e) and 1(f) by explicit formulae. But this would not work because the *σ*_{xx} effect is enormously path-dependent (Fig. 1(d)). If there were no path dependence, the failure points highlighted in Fig. 1(d) for different loading paths shown in Fig. 1(d) would have to coincide. But they do not, by far, and so the modeling must be incremental. This rules out the use of line crack models with variable *G*_{f}, particularly the LEFM and CCM (and the extended finite element (XFEM), as well), except for the rare situations in which the *σ*_{xx} is a priori known to vanish.

The effective widths, *c*_{f}, of the FPZ obtained with the size effect method are shown in Fig. 1(f). The increase of this width is intuitively explained in Fig. 2(a). The fracture process zone in quasibrittle materials (as well as polycrystalline metals) consists of microcracks of varied orientations. Static friction under moderate compression *σ*_{xx} initially helps to resist slip, but once the inclined shear microcracks slip, the friction drops, and all the slips combined widen the FPZ.

The gap test results shown by data points in Figs. 1(e) and 1(f) have been simulated in Refs. [20,21] by means of the crack band model [26,36] in which the microplane model M7 for concrete [31] has been used as the damage constitutive equation (the other curves are explained later).

It should be pointed out, too, that the *σ*_{xx} effect occurs also in plastic-hardening polycrystalline metals [37]. Their fracture has so far been treated as a line, but the size of grain boundaries in metals such as aluminum and their alloys indicates that the FPZ width should be at least several micrometers (≈2 − 50 *μ*m, [38]). Although the FPZ width is much smaller compared to the width of the hardening yielding zone (roughly 10 mm when *σ*_{xx} = 0), it may produce a significant effect of *σ*_{xx} on *G*_{f} (or *J*_{cr}), in addition to the well-known effect of the yielding zone size.

Figure 3 shows the results of the on-going gap tests of aluminum, which is an extension of the tests reported in Ref. [37]. By using the size effect law (SEL) for fracture transition to small-scale yielding [37], one can identify the *G*_{f} values (equal to *J*_{cr}) for various levels of *σ*_{xx}, as shown in Fig. 3. Considering that the effective yielding zone *r*_{p} remains the same for specimens scaled geometrically to various sizes with a fixed *σ*_{xx}, the size effect curves in this figure indicate [37] that *G*_{f} of aluminum depends significantly on *σ*_{xx}, as shown in Fig. 3(a). The size effect also delivers the *σ*_{xx} dependence of the yielding zone effective radius *r*_{p} (Fig. 3(b)), which is probably combined with the effect of *σ*_{xx} on the *c*_{f} of the FPZ (to distinguish these two effects, different kinds of tests are necessary). A further discussion of metals is beyond the scope of this article.

We note that the fracture energy obtained from the size effect method agrees with the critical value of J-integral at the peak load for specimens with different sizes. However, to incoporate the J-integral calculations into the crack band model (with a proper constitutive law), a contour to domain integral conversion is necessary [37]. As the contour must not cross the FPZ, it needs to be drawn carefully, considering the changing size of the FPZ due to *σ _{xx}*.

## 3 Overview of Phase-Field Concept

The phase-field concept was introduced more than half a century ago in physical chemistry, as a way of “spreading out”, for numerical purposes (e.g., Ref. [39], Fig. 2.2), the Heaviside step function that represents a sharp interface between, for example, the solid and liquid phases of the material (see [39, Fig. 4.1]). Similarly, a sharp crack may be regarded as a Dirac delta function, and the idea is to “spread it out” via a smooth “phase field”, φ, as shown in Fig. 4(a) [1]. In the simplest form, the function φ is chosen so that the energy minimization would cause the φ to decay from 1 at the crack line to 0 as $\phi \u221de\xb1x/w0$. However, φ physically represents no material damage characterized by some constitutive law and there is no associated Irwin’s material characteristic length *l*_{c}. The length scale *w*_{0} merely serves the numerical purpose of anchoring a sharp line crack with a point-wise tip to the mesh of elastic finite elements, so as to achieve mesh independence of crack direction. The phase parameter φ must be applied as a stiffness reduction factor in the phase-field band of a width *w*_{0} spanning over a sufficient number of finite elements. Thus, *w*_{0} represents a “fictitious” damage (or stiffness loss) of the finite element system, as shown in Fig. 4(a) by the exponential decay across the imagined (shaded) damage field toward the line crack. Typically, the elasticity matrix of the material is multiplied by *c* = 1 − φ at all finite element integration points in that band.

*t*

_{n−1},

*t*),

*n*= 1, 2, …, one calculates the state vector $Xn=[un,\phi n]T$ from the preceding state vector $Xn\u22121=[un\u22121,\phi n\u22121]T$ by a variational algorithm (Fig. 4(a)) described as follows:

_{e}(Helmholtz’s, isothermal) of the finite element system with its applied loads and imposed displacements is minimized while keeping the phase-field parameter

*φ*=

*φ*

_{n−1}(or its free energy Ψ

_{c}

_{n−1}) constant. This yields a system of linear equations for $u$. Then, the free energy Ψ

_{c}of the phase-field φ is minimized while keeping Ψ

_{e}constant. This yields another system of equations for the discrete values of φ. These minimizations may but need not be iterated. The result is the new vector $Xn=[un,\phi n]T$ at the end of the loading step. This alternative minimization (which is usually referred to as the “staggered scheme”) makes the algorithm efficient.

*G*

_{f}is the fracture energy of the material which, together with the elasticity matrix, forms the complete input of material properties,

*V*is the volume of the body domain Ω, and $\u2207$ is the gradient operator. Upon minimization of Ψ

_{c}, in which the expression in the parenthesis was postulated in 1998 by Francfort and Marigo [1], leads to the vanishing of the following expression as a function of the transverse coordinate,

*x*:

To validate the phase-field models, good fits of a few selected experiments, such as propagating curved cracks, have often been presented as supportive evidence. For example, Ref. [23, Figs. 19,20] presented an experimental validation for a standard compact tension specimen with a hole drilled on the side of the notch extension line including (1) the up-and-down curve of load versus deflection along with the curve of crack length versus the crack opening displacement and (2) the crack path in the same specimen. The hole causes the crack to run toward it, and the phase-field predictions matched these observations qualitatively and quantitatively. However, the question of scaling was left unanswered.

Another validation attempt used a standard double-edge-notched tensile specimen under axial tensile loads at the ends. Without a check of instability breaking symmetry, two opposite curved cracks are predicted to propagate from both sides symmetrically [41]. Even though the load versus displacement curve agreed with the test, the overall failure prediction was unsatisfactory since the energy analysis of the stability of the of postbifurcation path [42,43] showed that only a crack from one side can propagate.

An obvious weakness of the phase-field model is that the existing validated constitutive laws cannot be used directly. Instead, they need to be approximated in terms of free or potential energy functions, which cannot be fully equivalent and lead to complicated material characterizations in terms of multiple inelastic functions and complex loading/unloading responses.

Another weakness of the phase-field models is the arbitrariness in the choice of boundary condition imposed upon the microscopic force (i.e., the state variable that is work-conjugated with the phase variable). Mostly, the Neumann boundary condition is chosen, but with no meaningful physical justification.

## 4 Overview of the Concept of Peridynamics

*t*= time,

**x**= coordinate vector of material point,

**u**= displacement vector at the center point,

*ρ*= mass density,

**b**= body force vector,

*H*= volume or area within the assumed horizon,

**f**= assumed function characterizing the central interacting force between material points. The bond-based version assumes that all central–force bond interactions are independent of each other. This assumption leads to a fixed Poisson’s ratio equal to 0.25 [44] and limits further generalizations.

To overcome this problem, a version of peridynamics called “state based’” [13] was created upon introducing the assumption that the interacting forces within a horizon are inter-dependent and could be characterized by means of the states of forces, **T**, and states of deformation vectors in all the bonds emanating from the same center, **Y**:

It has been generally overlooked that (as pointed out in Ref. [16]) peridynamics actually represents a mathematically rigorous refinement and generalization of the “network model” proposed in 1977 by Burt and Dougill [45]. In that model, a field of random nodes is created, and each node is connected by bars to all adjacent nodes up to a certain maximum distance (which is what is in peridynamics called the horizon radius). After criticisms at conferences, these authors promptly switched to another approach.

In 2016 [16], several fundamental problems with the concept of peridynamics were identified, as follows:

*Lattice microstructure:*The basic physical problem with peridynamics is that it implies a lattice microstructure (Fig. 5(a)), which leads to a particle-skipping potential and a neglect of shear interactions that resist particle rotations. The lattice microstructure is unrealistic even for the state-based version. Although that version allows the equivalent continuum strain and stress tensors for each center point to be calculated and thus any constitutive damage model including the microplane model M7 to be applied, the underlying lattice characterized by central forces exists nonetheless.*Particle skipping interactions*are another physically unrealistic feature of peridynamics. These fictitious interactions are particularly detrimental for the modeling of compression-shear damage and fracture. In reality, the material contains finite particles, which not only displace but also rotate (Fig. 5(b)). Above the atomic scale, no potential communicating particle-skipping interactions exists. Indeed, interactions with the second and farther neighbors are communicated through the intermediate particles, by interparticle normal and shear forces. The shear forces resist relative rotations of particles (which is the basis of the lattice discrete particle model (LDPM)) [46,47]. Another oddity of the particle skipping interactions is that, in the bond-based version of peridynamics, the micromodulus governing the overall material stiffness has the queer dimension of MN/m^{6}.*Boundary conditions:*The peridynamic horizon that skips several particles inevitably leads to unphysical boundary conditions, as well as crack face conditions and FPZ border conditions (Fig. 5). The horizon of the nodes located near the boundary protrudes outside the boundary or across the open crack. To prevent it, such interactions must be deleted. This alone would, of course, significantly decrease the stiffness and strength of the boundary layer or crack surface layer (of thickness*δ*_{B}= horizon radius).To clarify the last point graphically, consider the diagram in Fig. 5(a), which shows the central force interactions within a circular horizon that emanate from the central point. When the horizon reaches beyond the boundary, the protruding central force interactions must be deleted. But this causes a major decrease of the stiffness and strength of the boundary layer. Also, it degrades the accuracy of the calculation of the deformation gradient in the correspondence state-based method. To counter this problem, surface-correction factors have been applied, as reported in Ref. [44]. But again this is only a partial remedy.

Another problem arises in applying the free surface boundary conditions. Since the displacement derivatives are unavailable for PD, a surface traction needs to be smeared over a certain volume. The displacement and velocity boundaries, however, are usually applied indirectly onto an imagined no-fail zone [44] considered as an undamaged zone of the material.

*Boundary conditions at crack faces and the FPZ dilemma*: The boundary problem is more severe and more complicated at the crack faces, and even more so at the boundary of the FPZ or the cohesive crack zone. Ahead of the crack and of its FPZ (Fig. 5(c)), the central force interaction along segment AA″ crossing the crack extension line is, of course, undisturbed. The interaction along segments or CC″ or DD″ crossing the open crack must be deleted, and the same surface correction factor must be introduced. A tougher predicament arises for the interactions along segment BB″ crossing the FPZ, for which the central forces crossing the FPZ get weakened only partly—strongly near the crack tip, where the damage is high, but weakly near the front of the damage zone (dashed ellipse in the figure), where the damage is still developing. A satisfactory, physically realistic resolution of the FPZ predicament looks impossible (this will be shown in the following comparisons with experiments).*Homogenization, localization, wave dispersion, fatigue:*Another physical problem, highlighted in Ref. [16], is that the same material length, the horizon, governs both the elastic homogenization and damage localization. Dispersion of waves by damage or plastic zones and by innate heterogeneity cannot be distinguished and independently controlled. However, experimental comparisons with wave propagation experiments are beyond the scope of this study. So is the fatigue loading of materials such as concrete [48], which is another unresolved problem.

The fact that peridynamics preserves the mass automatically has often been cited as an advantage over other computational models for fracture and damage. However, various meshfree (or hybrid) methods such as extended element-free Galerkin method, cracking particle, material point method and, especially, lattice discrete particle method may be easily programmed to preserve the mass. It is also no problem to program the crack band model to lump the mass of failed elements into surrounding nodes. All crack band simulations of projectile penetration run on Abaqus at Northwestern University [49–51], including those with dynamic particle comminution, as well as those in the wavecodes in national labs, were implemented this way.

The previous criticism of peridynamics [16] was strictly theoretical. No support by comparisons with experiments was given. This is remedied here. Comparisons are given not only for PD but also for PF models. In many cases, different versions of PD models, i.e., bond-based and state-based, give very different predictions for the same experiment. Even worse, subcategories of state-based models, which are ordinary and nonordinary (the latter are usually misnamed as the “correspondence-based” models [44]), suffer from the same deficiency. In this study, both “state-based” model types have been investigated. The performance of a bond-based model [52] will be briefly discussed in Appendix B.

## 5 Comparisons With Critical Classical Experiments Important for Practice

Seven models, which represent three main model types—PF, PD, and CB, are here compared to the most important test data for quasibrittle materials. They are as follows:

- –
*CB-M7*is the crack band model [26] with the microplane model M7 [25,31] as the inelastic constitutive law. Both the explicit and implicit versions (the latter contained slight updates [53]) are freely downloadable.^{4}For scaled specimens, the size of elements in the damage zone, which represents a material characteristic length,*l*_{0}, is kept the same for all specimen sizes. - –
*bPF*refers to the widely used (*basic*) phase-field model developed originally by Francfort and Marigo [1]. Its implementation in abaqus was carried out by Martínez-Panẽda and coworkers [54,55].^{5} - –
*bPD*refers to the (*basic*) ordinary state-based model with the critical-stretch damage law developed by Silling and Askari [11]. The use of this model and peridynamic models in general is facilitated by peridigm, a C++ library originally developed at Sandia National Laboratories [56].^{6}

After a lecture at the ASME-IMECE annual meeting on Nov. 1, 2021, which revealed large deviations of peridynamic and phase-field predictions from experimental evidence (see^{7}), some discussers claimed that the poor predictions were not inherent to the basic concepts of the models but were caused by oversimplification of the constitutive law. It was also claimed that the peridigm code [56] at Sandia lab was not the state-of-the art and that it would fit the data well if it were combined with a certain new material law for concrete. But it turned out that this material law and its implementations for peridynamics were unpublished and unobtainable. Apparently, this and other models were designed to capture one important aspect of the targeted material, i.e., concrete [57,58], while ignoring others. For PF models, some improvements to fit few particular tests have also been published, but again, with a focus on demonstrations of some mixed-mode loading cases instead of checking the scaling and the full range of published experimental evidence [40,59]. Therefore, to examine such claims and purported improvements, four more models and their comparisons to experiments had to be added:

- –
*CB-Gr*is a tensorial damage constitutive model implemented within the same finite element framework as*CB-M7*, except that M7 has been replaced with the second-generation concrete damage plasticity model (CDPM2) developed by Grassl et al. [60]. This model is an update of Ref. [61] and represents arguably the best plastic-damage constitutive models of concrete formulated in the classical way—in terms of tensors, their invariants, and loading surfaces.^{8} - –
*PD-Gr*is the implementation of the CDPM2 constitutive law into Peridigm based on the nonordinary (correspondence) state-based formulation. - –
*PDba-Gr*is similar to PD-Gr. However, in this new development, Bazilevs, Behzadinasab and Foster [62–64] corrected the definition of the deformation gradient (1) by eliminating the zero-energy mode of instability plaguing the conventional PD theory; (2) by limiting the nonlocality of the model via the reduction of the horizon to the nearest neighbors, which was intended to avoid the overpredictions of the damage and plasticity in the comparisons that follow; and (3) by treating peridynamics as a discretization method rather than a model per se (see a brief summary of the bond-associated PD in Appendix E).^{9}

A further study dealt with two recently developed PF models [40,59,65] in which the constitutive laws were enriched to capture the cohesive stress in tension, shear and the frictional contact with dilatancy (see Appendices C and D). Their improvements, however, were aimed to correct poor fits of only one or a few types of experimental data. It was not reported whether the fits of other types of tests got improved or worsened. These models are as follows:

- –
*PF-Wu*is an enhanced phase-field model incorporating the cohesive zone theory developed by Wu et al. [40,59] (see Appendix C). This model can be downloaded from.^{10} - –
*PF-Choo*is an improved phase-field model proposed by Fei and Choo [62], which is intended to capture the mode-II fracture and the frictional slip between cracked surfaces (see Appendix E).^{11}This model, though, could be applied only to the confined compression of a slab experiment since it is two dimensional.

The PF and PD models investigated here are the current state-of-the-art representations of each model type (which are either publicly uploaded on open-source domains or have been generously shared by the authors). The purpose of the direct comparisons of *CB-Gr*, *PD-Gr*, and *PDba-Gr* is to clarify whether the deficiencies of PD are due to inadequacy of the material constitutive model for broader applications, as some developers now claim, or to its basic concept. For the comparisons to be fair, the *PD-Gr*, *PDba-Gr,* and *CB-Gr* simulations of all experiments share the same set of material parameters. The horizon radii in the peridynamic counterparts were taken as 1/2 of the characteristic element size in the crack band model. We chose the horizons for PD model in the same manner and set the critical stretch to match the fracture energy, *G*_{f}, according to Silling et al. [13]. The influence function in all PD models was considered constant. These requirements ensure that the results for various models be comparable. However, the horizon dependence of PD models and the *l*_{0} dependence of PF models will be addressed in detail later. The calibration process for the present models is summarized in Appendix A.

In addition, we posit that a model that is updated to agree with only one or few test types but is not compared to other relevant tests of different sizes is questionable. What is needed is a model that can adequately characterize various aspects of a given material for general applications in practice. Selective ad hoc model adjustments that improve one type of response but may compromise another are hardly useful.

Most of the comparisons presented here (except one) deal with *concrete*, for four reasons:

The diversity of available fracture test types is for concrete much greater than for any other material.

The size effect data are plentiful and cover a broad range of structural sizes.

Concrete is an archetypical quasibrittle material, and its fracture behavior is similar to all others (rocks, fiber composites, tough ceramics, sea ice, rigid foams, bone, stiff soils, etc., and most materials on the sub-micrometer scale).

In terms of nation-wide usage, financial outlays, and the associated CO

_{2}emissions, concrete outstrips every other material by far (in terms of CO_{2}emissions, the worldwide production of cement is about to exceed all fossil fueled transportation, while suppression of fracturing would diminish ingress of corrosive agents into concrete structures, which would extend their woefully inadequate lifetimes and thus mitigate the future demand for cement).

It is also worth noting that the hardening plastic metals exhibit at macro-scale a few phenomena that are similar to quasibrittle materials—the vertex effect [32,66–68] or, for small-scale yielding, the scaling of structural strength [37], to name a few. Therefore, a model that can capture accurately the plastic and damage responses is indispensable for any material.

### 5.1 Type 2 and Type 1 Size Effects.

Scaling is the quintessential characteristic of every physical theory. In fracture mechanics, most important is the size effect, which represents the dependence of the nominal structure strength *σ*_{N} on size *D* of geometrically scaled structures (*σ*_{N} = *P*/*A*, where *P* = maximum value of the applied load or load parameter, *A* = *bD* in the case of two-dimensional scaling, *b* = specimen thickness, *D*, *A* = characteristic cross section dimension and area, measured homologously). In elasticity with strength limit and in plasticity, there is no size effect (except statistical), but in fracture and damage mechanics, it is an essential feature.

*G*

_{f}, and to the FPZ size was found in 1987 [70] and in 1990 [35]. The SEL describes the gradual transition from strength-based failures to brittle failures at increasing structural size

*D*. In 1990, Eq. (6) was reformulated as Eq. (7), where

*g*(

*α*) is the dimensionless energy release function of LEFM [35]. This was achieved by means of second-order asymptotic matching [71] of the approaches to the large-size LEFM asymptote (dashed in Fig. 6(a)) and to the horizontal small-size asymptote corresponding to plasticity ($E=Young\u2032smodulus$,

*G*

_{f}= material fracture energy,

*f*

_{t}= tensile strength of material,

*B*= dimensionless geometry constant,

*c*

_{f}= material characteristic length, roughly equal to 40% of the FPZ length,

*α*=

*a*/

*D*,

*α*

_{0}=

*a*

_{0}/

*D*,

*a*= crack length including the cohesive zone,

*a*

_{0}= length of open crack or notch).

The type 2 size effect occurs when there is a long enough notch creating a positive geometry [71], or when the structure geometry allows large stable crack growth prior to the maximum load [27,71], as typical for RC, e.g., RC beam or slab shear [71]. The initial stable growth requires negative geometry (i.e., [∂*K*_{I}/∂*a*]_{P} < 0, where *K*_{I} = stress intensity factor and *a* = crack length). Failure under controlled load *P* occurs dynamically when the geometry changes to positive ([∂*K*_{I}/∂*a*]_{P} > 0). In type 2, the geometry is initially positive (as in 3PB specimen) and material randomness affects significantly only the coefficient of variation of *σ*_{N}, while the size effect on mean strength is deterministic (energetic) for all structure sizes.

*P*dynamically as soon as a macro-crack initiates from a fully formed FPZ. For small sizes

*D*, stress redistribution occurs due to the initiation of distributed damage within the representative volume element (RVE), which controls the small-size power-law asymptote. For very large sizes, the type 1 size effect transits to the power law of the Weibull statistical size effect, provided that fracture can initiate randomly from many nearly equally stressed locations (this does not happen for three-point bend beams and thus cannot affect the mean maximum load). In addition, such a transition from the deterministic to the statistical (Weibull or Gauss–Weibull) size effect occurs, if at all, at sizes larger than those of the present tests [71,72]. Consequently, the law of type 1 size effect [71,73], derived by the matching of these two asymptotes, reads:

*n*/

*m*= 1/12,

*r*= 1 in typical concrete structures. For small sizes considered here, the statistical size effect is negligible, which ensues by setting 1/

*m*= 0. The type 1 to 2 transition is more complicated [33] and is not needed here.

The transition of the size effect from one power law to another is very broad. It is governed by the material characteristic length *c*_{f} (typically about 0.4 of Irwin’s length). The *bPF* model involves no material characteristic length (*w*_{0} serves only as a numerical check for regularization). On the contrary, the horizon radius in the **bPD** model is the material characteristic length, which is what should govern its size effect transition.

Figure 6 shows, by circle points, the size effect test data of Hoover et al. [74,75] for type 2 and type 1. However, the individual data points are shown only for type 1. For type 2, the error bars (i.e., ± standard deviation) had to be used because low scatter makes the plot of data points too congested. No tests were made for *D* = 1000 mm, but the size effect was extrapolated up to *D* = 1000 mm, using the calibrated SEL. It should be mentioned that Grégoire et al. [76] also performed similar experiments on specimens of a different aspect ratio. They are not used here due to their higher scatter and narrower size range but will be considered in Appendix B.

To be comparable, all the predicted curves are calibrated to pass through at least the first point representing the results of the specimens with the smallest size. To illustrate more clearly the deviation from the size effect law, the simulations are run for *D* = 1000 mm beyond the size range of data (circle points), for both size effect types. The following points should be noted:

- –
For type 2, the solid curve shows the optimum fit by the SEL (Eq. (6)). The finite element predictions of the CB model with microplane model M7 are represented by the dotted lines passing through diamond-shaped points. They are satisfactory.

- –
Both phase-field predictions are shown by the crosses and triangle points connected by line segments with the same colors. For type 2, these predictions deviate significantly from the data (Figs. 6(c) and 6(d)). Note that the deviations in

*PF-Wu*are larger than those reported in Ref. [77] for the same model. This stems from the recalibration of the model to obtain more accurate results for small-size specimens and to strike a balance between the accuracy of type 1 and type 2 size effects, provided that, importantly, the same set of parameters is used. The fact that no transition exists and that these models both result in straight lines, corresponding to power laws, is incorrect. For the*bPF*, the slope in the log-log scale is −1/2 (Fig. 6(c)). This is correct only for LEFM, which is not the case for these test data. For*PF-Wu*, the straight line slope is less steep than −1/2, which is thermodynamically inadmissible because such a power law implies a zero energy flux into the crack tip [27,71,78]; see Appendix H. - –
The PD models for type 2 not only deviate from the experimental data but also produce wrong curvatures in the log-log scale plots (Fig. 6(e), note that in the linear scale plots such deviations can, and have been [52,77], misjudged as data scatter).

- –
For type 1 (Fig. 6(d)), the results of

*PF-Wu*model lie within the scatter range of the experimental results, and so no limitation to*PF-Wu*validity is indicated here. Similar to type 2,*bPF*tends to result in a more brittle material behavior (manifested by a smaller characteristic structure size). Moreover, the curve shapes of these models for type 1 size effect disagree with the mean data trend (Fig. 6(d)), and the*bPF*model even diverges away from the data scatter band (note that type 1 has a much smaller slope than type 2). - –
Surprisingly, the three peridynamic models produce results with greatly diverse trends and values (Fig. 6(f)). The

*bPD*model tends to capture accurately the trend, in proximity to its*CB-M7*counterpart, while*PDba-Gr*shows a significant deviation from experiments. However, Fig. 7(c) shows that the correct trend of the*bPD*model is associated with a wrong mechanism. The extra energy is dissipated via branched cracks instead of a gradual dissipation via a gradual softening law. The*PD-Gr*model, however, shows an increasing trend at larger*D*(Fig. 6(f)), which is physically impossible. - –
To explain what happens in the PD models, the principal strain contours that they deliver are compared with the CB counterparts. As obvious from Fig. 7, the breadth of the FPZ in

*PD-Gr*for type 1 is wider than it is in other models, even though the horizon is of the same size. This overestimates the nonlocality of the damage zone and leads to a spurious growth of the FPZ across scales. The root cause of this overestimation may be attributed to the zero-energy mode of instability discussed in detail in Refs. [16,42,63,79]. The*PDba-Gr*shows a more localized damage pattern, but it results in a longer crack at the peak load (Fig. 7(g)). - –
It should be here noted that similar phenomena for both type 2 and type 1 size effects were observed with a bond-based peridynamic model studied by Hobbs et al. [52]; see Appendix B, Fig. 24. For type 2, a power-law trend along with significant deviations from experimental data [76] was observed. For type 1, the strengths of scaled structures remained approximately the same across sizes, which was attributed to a lack of statistical size effect calculations. However, this cannot be true because the statistical size effect would require many possible potential failure locations with random material strength, while, in three-point bending, the failure can start only at

*one*location, the bottom midspan [71,80]. In addition, the specimen sizes in this comparison are close to the deterministic type 1 asymptote.

### 5.2 Concentrated Mode II Shear Fracture (In-Plane).

The double-notched four-point loaded Iosipescu beam shown in Fig. 8(a) was used in the 1980s to prove that Mode II shear fracture in concrete exists when large strain energy is getting released from a narrow strip of material. In addition, this test also reveals that the crack-parallel compression affects not only mode I but also mode II (shear) fracture although, contrary to mode I, it decreases the fracture energy, which is expected because a high-enough parallel compression can alone produce splitting fracture. Normally the opposite loads on Iosipescu beam are placed farther away from the notch, in which case the computations indicate two curved cracks called mixed-mode failure although they are of pure mode I at their tips. They diverge to opposite sides of the crack line and can be reproduced well by many models. In experiments, only one curved crack occurs because the path-stability criterion [81] excludes both from propagating simultaneously.

More revealing was the 1986 test [79] in which the loads were placed as close to the notch mouths as possible without shearing of the crack mouth corners. This loading was shown to produce a planar pure mode II shear crack [79]. Significant crack-parallel compression exists, and its role is documented by this test. The simulation results are presented in Figs. 8(c)–8(i), to be compared with the experimental observation in Fig. 8(b).

- –
The

*CB-M7*and*CB-Gr*’s stress–strain relations both perfectly reproduce the pure mode II crack creating sliding surfaces in contact. Both models predict the crack to start from the notch. Also, both predict the same strain level at the onset of damage localization and both reproduce the planar Mode II crack perfectly. Because of the finite width of the crack band, they, of course, cannot reproduce explicitly the sliding surfaces. These surfaces form not at the maximum load but only during the postpeak, after full softening. - –
The results of

*bPD*and*bPF*models (Figs. 8(e) and 8(g)) might look satisfactory except that the fracture damage band deviates to the side of notch under the loading pad, which is incorrect. The models do not distinguish shear from compression crushing, which stems from the inability to discern various damage mechanisms. Hence, they incorrectly predict splitting cracks under the loading pads. - –
The

*bPF*and*bPD*models result in similar incorrect cracking patterns (Figs. 8(e) and 8(g)). By using a single scalar fracture criterion, they both favor the opening mode cracks. This causes a diffused cracking pattern instead of a localized shear along the notch extension line. - –
The

*PF-Wu*model captures well the straight propagation of mode-II cracks subjected to concentrated load. This is so thanks to the consideration of multiple principal stress components in the damage phase variable, instead of one scalar stress value, or*G*_{f}(Fig. 8(f)). - –
In

*PD-Gr*, we encounter again a delocalized damage (Fig. 8(h)). In addition to the wide spread of the zone of principal damage strain, the damage zone grows, curiously, in a direction opposite that to that in*bPD*, away from the loads, and its breadth is far too big. - –
Among all the PD models, only the

*PDba-Gr*produces the correct crack path.

### 5.3 Compression-Torsion Fracture in Mode III (Antiplane).

This test (Fig. 9) was used in the 1980s to prove that mode III shear fracture in concrete does exist [82]. A circumferential notch was cut in a cylindrical specimen, and a torque was applied. It produced a perfectly planar mode III crack; see Fig. 9(b), left. Then axial compression strain *ɛ* = −0.001 (for which the axial force *P* ≈ was about 30% of the uniaxial compression strength) was applied. A simultaneous torque then produced a conical mode III shear failure shown in Fig. 9 (center). The axial strain level that marks the change in the fracture mode was not measured and has been extracted from calibrated *CB-M7* simulations. The simulation results of the two crack band models are shown in Figs. 9(c) and 9(d).

- –
The explanation of the change of mode is twofold: (1) Axial compression causes the energy dissipation by mode II fracture and frictional slip to be higher on the planar crack than conical surface, and (2) in the light of the gap test and the splitting fracture role in the mode II test, the compressive axial force component parallel to the conical surface must have reduced the mode III fracture energy (in fact, the analysis of this test might be used to sort out the role of this component quantitatively).

- –
None of the models, except the two CBs, reflects the plane-to-cone transition when

*ɛ*increases. - –
The flat mode III shear crack is reproduced perfectly in Figs. 9(c) and 9(d). For torsion with axial compression, there is an indication of mode III cone formation, clearer for

*CB-M7*than for*CB-Gr*, which indicates secondary fracturing to spread to the side of the cone (Figs. 9(c) and 9(d)) further shows the crack band prediction for a tripled axial compressive strain*ɛ*= −0.003, for which the torque, in a cylinder of finite length, should lead to a near-cylindrical failure surface, although this has yet to be checked by experiments). - –
The

*bPF*shows, for*ɛ*= 0, a curved cylindrical surface propagating from the notch to the upper half of the specimen and a detachment of the top part of the specimen. When*ɛ*increases to a higher level, materials are removed from both halves, but the mode III crack is missed in all the cases. All of this is incorrect.*PF-Wu*incorrectly shows that the damage occurs only near the circular edge of the notch and predicts no shear crack formation. At $\epsilon =\u22120.3%$, there is more diffused damage formed around the notch edge, but again no shear crack. - –
The

*bPD*incorrectly indicates the formation of an inclined crack farther away from the notch edge, both for*ɛ*= 0 and −0.001. Note that these cracks form abruptly due to the brittle nature of this model. For*ɛ*= −0.003, the results incorrectly suggest a plane crack propagating from the notch edge, accompanied by axial splitting. - –
Because the

*PD-Gr*model cannot form a localized damage surface, the notch in this model cannot propagate and the material points adjacent to the notch begin to suffer from spalling. At higher*ɛ*, this failure mode is suppressed. Instead, a curved cylindrical surface is formed, similar to the case of the*bPF*model. - –
The

*PDba-Gr*model, however, can capture the cracks emanating from the notch and forming a flat horizontal surface. Nevertheless, at high*ɛ*, the damage starts to spread out rather than concentrating to a localized damage pattern.

### 5.4 Unconfined Uniaxial Compression and Shear Band Growth.

The axial compression test of a cylinder (diameter 10.1 cm and length 20.2 cm) is used at construction sites as the standard quality check of concrete. The failure is triggered at *P*_{max} by the propagation of inclined shear bands consisting of dense axial splitting cracks (Figs. 10(a) and 10(b)), which look after *P*_{max} like shear cracks. The cylinder ends are assumed to be lubricated so that they can slide perfectly (if there is significant friction, the damage would spread widely and get combined with axial splitting). Figure 10(c) compares the curves of axial average stress versus average strain.

- –
Both

*CB-M7*and*CB-Gr*fit the response curves well (Fig. 10(c)) and both show the formation of inclined shear bands (Figs. 10(d) and 10(e), including the postpeak response (which can be observed only in a very stiff frame with fast servocontrol, as discovered in 1963 [83–86]). However, the softening slope indicated by*CB-Gr*is steeper than in experiments due to an underestimation of the frictional resistance. Note that the shear band is subjected to crack-parallel compression, and the fact that CB models can capture it helps getting correct results. - –
The

*bPD*gives the correct peak stress (Fig. 10(c)), but the strain at peak stress is about 70% of what it should be, and the steep postpeak stress drop is wrong. The problem is that*bPD*does not have a state variable that can control the gradual loss of material cohesion. Therefore, it predicts a rather abrupt release of energy in the cylinder when the peak load is attained. The specimen is wrongly predicted to fail by distributed fragmentation, with no sign of shear band formation or shear interaction between adjacent material points. Also, the lack of shear interaction in*bPD*prevents capturing the frictional slip, which leads to an incorrect sudden stress drop (such a drop used to be seen prior to 1963, but this was due to energy release from the soft testing machines used in those times, rather than to the material properties). - –
The

*bPF*model gives unlimited elastic response, with no sign of damage (or any growth of phase variable φ). This is not surprising since*bPF*can reproduce neither compression fracture nor shear slip. It cannot develop a localized shear band even when a small flaw is introduced in the FE system. The flaw produces a small bump (shown in Fig. 10(c)), but then the stress continues to rise without limits (and the phase variable φ increases to 1 throughout the entire specimen). Obviously, the fracture of quasibrittle materials under uniaxial compression can be captured only if shear boundaries are present in the material model (which will also be seen in the next experiment). - –
The

*PF-Wu*model, however, takes advantage of the deviatoric stress response to generate damage (which, however, accumulates mainly at the surface in contact with the loading plates). Nevertheless, this requires a higher*G*_{f}, i.e.,*G*_{f}= 1200 N/m, to attain the observed compressive strength. The lack of the uniaxial compressive boundary was shown in Ref. [59, Fig. 3]. To remedy this problem, the frictional and deviatoric stress would have to be somehow represented in the phase-field equation, which seems to be difficult. - –
The

*bPD*model incorrectly shows the uniaxially compressed specimen to fail by the removal of particles under the cones. This in turn is a consequence of the excessively brittle nature of the hypothesis of critical stretch of interparticle connections. - –
Both

*PD-Gr*and*PDba-Gr*show a dense network of diagonal crack surfaces, while the damage field in*PD-Gr*is more diffused, which is not quite realistic. - –
Even though the nominal stress–strain curves produced by

*CB-Gr*,*PD-Gr*, and*PDba-Gr*models are close to each other, the failure mechanisms are not the same (see Figs. 10(e), 10(i), and 10(j)). This might be explained by the material shear resistance that can arise either from a single inclined plane or multiple inclined mutually orthogonal planes (see Fig. 10(l)), as long as the total surface areas of frictional contact are the same.

### 5.5 Confined Compression of Slab.

Similar conclusions can be drawn for a plate of sandstone (Figs. 11 and 12) confined in two dimensions (2D). As described in Ref. [87], a plate specimen, confined in cross-thickness direction and sliding between two stiff steel plates, was compressed in vertical direction and subjected to controlled lateral confining pressures *p*_{c} = 0 and 8 MPa (Fig. 11(a)). In the case of PD models, due to the complexity of applying a uniform pressure transverse to the specimen, *p*_{c} is indirectly generated by a layer of elastic-perfectly plastic material that has a yield strength of 0.01 and 8 MPa.

- –
In the diagram of average axial stress versus strain (Fig. 11(b)), the down-and-up trend after the peak in the experiment at

*p*_{c}= 8 MPa shows a hardening due to the increase of confining pressure as the bolted plates resist the dilation of the specimen. Therefore, in this regime,*p*_{c}may have exceeded the expected pressure of 8 MPa. - –
When only the response for

*p*_{c}= 0 MPa is used for the calibration, the*CB-M7*prediction for 8 MPa is satisfactory (especially in view of the strange postpeak data). - –
While the

*CB-M7*model shows a ductile response, the*CB-Gr*exhibits at both*p*_{c}= 0 and*p*_{c}= 8 MPa an excessively brittle response, with the axial splitting cracks forming an inclined band and distributed damage forming at*p*_{c}= 8 MPa under the loading pads. Nonetheless, both show a change in the direction of the localization pattern when*p*_{c}is increased. - –
The

*bPF*model gives, again incorrectly, an unlimited stress rise for both*p*_{c}levels (Fig. 11(b)). Figure 11(e) shows the absence of localized bands despite the introduced void, which can be explained by the lack of a state variable controlling the deviatoric damage growth. - –
For

*PF-Wu*, there is a damage pattern similar to the uniaxial compression. Unlike the*bPF*model, the*PF-Wu*model considers the growth of deviatoric stress, but the absence of a variable mode II fracture again leads to a wrong failure mechanism, i.e., damage near the loading pads (Fig. 11(f)). A more serious problem is the average strength decrease when*p*_{c}increases. This can be explained by inappropriate constitutive representation of triaxial and biaxial stress boundaries, in which the presence of a lateral pressure increases the damage. - –
The foregoing poor predictions are improved for

*PF-Choo*(Fig. 11(g)). This is largely due to its capability to account for the growth of phase variable φ as a function of shear displacement and frictional contact (see Appendix D for more detail). However, the opening mode and antiplane shear responses now become questionable. - –
It appears that the common problem for all of the present phase-field simulations is their inability to account for all important damage mechanisms and for the transitions between them. By using only one phase-field variable, several damage mechanisms, if considered, must be tuned concurrently, which leads to spurious damage development. Recently, models with more phase-field variables [88,89] have been formulated. However, the questions of the appropriate number of phase-field variables to create a predictive model, and of the interaction of their respective damage mechanisms, remain unanswered.

- –
The damage patterns at peak load in all PD models are unreasonable for any pressure

*p*_{c}, even though the calibration ensures correct*P*_{max}at*p*_{c}= 0. For*bPD*, instability occurs once the compressive stress reaches the nominal strength value at*p*_{c}= 0. However, the typical diagonal damage pattern can still be observed. When*p*_{c}= 8 MPa, no clear pattern of cracking and damage emerges. Rather, the material points at the corners are stripped away by the biaxial loads. This causes an immediate drop of load at*p*_{c}= 8 MPa, which is unrealistic. - –
The

*PD-Gr*model produces an identical nominal strength at low and high*p*_{c}. However, the damage develops and remains concentrated around the initial flaw, and no localized damage band forms even after the peak load. This can again be explained by the emergence of zero-energy mode instability, leading to an excessively wide spread of inelastic strain. - –
The

*PDba-Gr*model overestimates the (mode-II) energy release due to its denser network of localized damages. In addition, the dense network results in no well-defined diagonal damage bands known to dominate the shear failure.

### 5.6 Confined Compression of Cylinder.

Although this is not a fracture test, it is important for some fracture simulations. When a projectile impacts a concrete bunker, first very high confining pressures briefly arise under the nose of the projectile. Under very high pressures, concrete becomes a plastic material allowing shear angles such as 70 deg without any fracturing [49,90]. The plastic deformation dissipates much impact energy prior to fracturing (and, during subsequent fracturing, the comminution of particles dissipates more energy [50,51]).

To test this phenomenon, a small cylindrical hole (19 mm in diameter) within a large diameter tungsten carbide vessel was filled with small aggregate concrete [91] (see Fig. 13). This setup provided a near-rigid confinement with virtually zero lateral strain. The cylindrical specimen confined in the hole was subjected to high axial load developing pressures up to 2000 MPa, 50 times higher than the uniaxial compression strength. Wall lubrication was provided by copper foil sliding under high pressure plastically. Under such perfect confinement, concrete exhibits no softening. The tangent modulus *E*_{t} remains always positive. It first decreases due to collapsing voids but later increases due to void closing, greatly exceeding the initial elastic value *E*_{l}. At 2000 MPa, the initial porosity is reduced roughly to one half. Therefore, the unloading modulus *E*_{u} becomes several times higher than its initial elastic value *E*_{l}. Missing this behavior, the impact penetration depths or wall exit velocities get grossly overestimated [51].

- –
Only the two CB models are able to capture the observed behavior. They succeed thanks to a constitutive model that involves a state variable simulating the change in void volume fraction. Yet, in

*CB-Gr*, the change of modulus*E*_{t}is modest, while*CB-M7*not. - –
The inability of

*bPD*and*bPF*to capture complex triaxial stress states suppresses the development of inelastic strain at material points. In these models, the material keeps compressing elastically. - –
The

*PF-Wu*model again shows a weakened resistance under triaxial stress states. The specimen is predicted to fail at excessively low*σ*even when an unrealistically high*G*_{f}, equal to 2500 N/m (instead of 100 N/m), is assumed. - –
Despite using the same constitutive law,

*PD-Gr*and*PDba-Gr*show two different damage patterns and differ greatly in the nominal stress–strain curves. In both models, damage accumulates on the surface of the specimen and comes to play right at the beginning of softening. During unloading, however,*PD-Gr*shows a significantly milder slope, or lower unloading modulus*E*_{u}, than the CB counterpart and the data.

### 5.7 Vertex Effect.

The term “vertex effect” stems from the theory of incremental plasticity, in which the response to a load increment tangential to the current yield (or loading) surface must be elastic unless the yield surface has a corner, or vertex, at the current state point in the stress space. In the case of a test cylinder gradually loaded by uniaxial compressive force *P* (Fig. 14(a)), all theories of plasticity based on tensor invariants (as well as most tensorial damage theories), the loading surface is normal to the *σ*_{xx} axis and crosses the axis smoothly, with no corner at the crossing. This means that if one suddenly applies a shear stress increment *σ*_{xy}, which is best done by means of torsional moment *M* (see the loading path in Fig. 14(a)), then the incremental stiffness, according to incremental plasticity, is supposed to be elastic. However, it was experimentally demonstrated, first for metals in the 1950s [32,92], and later for concrete [66], that the incremental response is much softer than elastic. This is called the vertex effect (for a detailed discussion, see Sec. 25.4 in Ref. [93]).

Generally, the microplane model is expected to work since it implies a very large number of stress–strain boundaries on the microplanes, analogous to the tensorial loading surfaces in incremental plasticity. Vectorial microplane boundaries of different orientations exist everywhere in the stress space (i.e., the space of nine tensorial stress components). They can be seen as the vectorial analogs of the tensorial loading surfaces. The fact that they are vectorial rather than tensorial provides tremendous simplification and helps physical insight. In the test of concrete, the vertex effect begins when the axial force reaches about one half of the compression strength. It becomes particularly strong in the postpeak softening (which has no parallel in plastic metals).

- –
The microplane model M4, when used in CB-FE, was shown in Ref. [25, Fig. 8] and earlier in Ref. [67] to reproduce closely the vertex effect data observed on concrete in [66]. Here, the same data (circle points in Fig. 14(b)) are compared with the

*CB-M7*[25,31,53]. The*CB-M7*vertex effect is seen to be satisfactory though somewhat too strong (probably due to microcracking zone being larger than in reality). - –
Thanks to the scalar damage variables, the

*CB-Gr*does exhibit the vertex effect, too, and does so slightly better at a higher compressive strain (black curve in Fig. 14(b)). Figures 14(c) and 14(d) shows that*CB-Gr*generates less microcracking damage than*CB-M7*at the same level of uniaxial compressive strain, and that the strain localization starts at multiple planes (Figs. 14(c) and 14(d)). - –
While the CB models predict correctly the inclined localized shear bands at progressively higher loads, both PF models predict damage bands to form on planes normal to the cylinder axis, i.e., cracks are generated at a low inclination at all levels of the compressive strain. This leads to the wrong damage patterns seen in Figs. 14(e) and 14(f). For both

*bPF*and*PF-Wu*models, the lack of frictional sliding and microcrack splitting mechanisms prevents them from reproducing the change in damage pattern at different compressive strains. This lack leads to an elastic loading increment when the torque is superposed. - –
Due to the abrupt failure of the

*bPD*model under compression, the torsional stiffness at strain $\epsilon >0.22%$ is not available. The cracking pattern of*bPD*disagrees with observations in experiments, showing a localization into a twisted surface instead of a tilted plane, even at low compressive strain. At higher compressive strain, most of the particles are abruptly removed before applying the torque, which causes stability loss (seen as a massive damage volume in Fig. 15(b)). - –
The localized damage pattern and the reduction of the initial torsional stiffness of concrete produced by

*PDba-Gr*agree with its CB counterpart and with experimental data. The success of this prediction is due to the fact that this experiment depends only on the correctness of the tensorial constitutive law and is independent of fracture properties as well as the characteristic length. - –
Even though the torsional stiffness declines as

*ɛ*increases, the wide spread of microcracks, due to the underlying lattice structure, causes the*PD-Gr*to underestimate the torsional stiffness reduction.

It is worth noting that the vertex effect is important for modeling many practical situations in which the shear loading follows high uniaxial or biaxial compression, and generally when the principal stress directions rotate during loading, which is commonplace (examples are the impact, or seismic shear loading of a vertically compressed column or wall). The computer codes for the existing popular plasticity models based on invariants miss the plastic vertex effect. As long as the computer analysts insist on tensorial models with invariants and one or two loading surfaces, this remains an insurmountable deficiency. After initial frustrations [93], the vertex effect has simply been ignored after 1980 in computational plasticity, although this is acceptable only for near-proportional loading. Only recently has the question been revived and discussed again [68].

### 5.8 Diagonal Shear Failure of Reinforced Concrete Beams and Slab Punching.

This has been a formidable problem, arguably the most difficult one in fracture mechanics. It has been studied for over a hundred years and by fracture mechanics for four decades. Many RC structures collapsed unexpectedly in this manner. The observed load capacity, with its size and shape effects, has never been successfully predicted by LEFM, nor by the cohesive crack model. Why?

The answer has recently been given by the gap test. The line crack models, such as LEFM or CCM, can predict correctly the mode I initiation of the main diagonal crack, which happens at about 0.5 *P*_{max}. But at the maximum load, *P*_{max}, at which the crack already extends through about 80% of the cross section (Fig. 16), the compressive stress *σ*_{xx} along the crack, transmitted along an imagined compression “strut” parallel to the crack (the darker strip above the main crack in Fig. 16(b)), nearly attains the uniaxial compression strength limit $fc\u2032$. At that moment, the crack-parallel compression reduces the mode I fracture energy to almost 0, as recently revealed by the gap test [20,21]. The failure finally occurs by the propagation of a compression-shear band of splitting cracks across the top of the compression “strut” and the push-up of the triangular wedge above the crack tip (of length *a*_{c} marked in Fig. 16(b)) [94,95]. This is a failure mode that follows closely the type 2 size effect, as confirmed by many tests and physically explained by the release of strain energy from the compression “strut,” which increases quadratically with the beam size.

Enormous funds have been spent worldwide on the tests of this kind of failure. The database, which was assembled in 1984 at Northwestern and was later more than doubled by the ACI Committee 445, contains now over 800 test series collected from the literature [96]. The scaling of the shear failure is particularly well documented by the recent tests of Syroka-Korol and Tejchman [97]. Geometrically scaled RC beams of three sizes (or depths), with *D* = 200, 400, 800 mm (Fig. 16(a)), were tested. The observed diagonal shear cracks are plotted in dimensionless coordinates in Fig. 16(a), in which all the beams coincide. The crack paths for different sizes are seen to almost coincide as well, and the *P*_{max} occurs when the crack tip reaches almost the same point in dimensionless coordinates.

As confirmed generally by the database, the RC beam shear strength *σ*_{N} follows the size effect law as closely as the data scatter permits one to discern; see Fig. 16(b) (based on this fact and on a similar analysis of slab punching fracture, the SEL has been incorporated into the ACI-318/2019 concrete design code). The simulations by the *CB-M7*, shown in Fig. 16(b), agree with the test data within their inevitable scatter.

The performance of the models is compared with the experiments in three aspects: the development of cracks when *P* increases, the consistency of the main crack paths at the peak load for various beam sizes, and the prediction of the structural strength of scaled beams.

- –
The fracturing starts by many distributed cracks at

*P*_{max}/3. At*P*_{max}, which eventually localize into one main diagonal crack (Fig. 17(a)). - –
Both CB models show the formation of distributed cracks at

*P*≈*P*_{max}/3 and of a large diagonal crack at*P*_{max}; see Figs. 17(b) and 17(c). However, the distance between the cracks in*CB-Gr*is larger than typical, and some of these cracks still contribute to the dissipation at*P*_{max}, instead of a single diagonal crack as in*CB-M7*. - –
A common problem of

*bPF*and*PF-Wu*is an unrealistic pattern of crack formation (Figs. 17(d) and 17(e)). Both incorrectly show the formation of a single bending crack at the middle cross section of the beam, followed by extensive debonding of steel bar from concrete, and only after that a wide inclined zone of diffused damage leads to failure. Moreover, for all sizes, the*bPF*model shows no peak load, wrongly predicting a continuing load increase without failure. This is a major deficiency. It gets only slightly mitigated in*PF-Wu*, but the lack of a major crack at the peak load still leads to gross overestimation of the beam shear strength (Fig. 17(d)). - –
The

*bPD*model incorrectly indicates the formation of several scattered cracks at*P*_{max}/3. At the maximum load,*P*_{max}, these cracks lead to widespread damage at the wrong place. The brittle nature of*bPD*further triggers multiple splitting damage under the loading pads. After that, debonding happens between concrete and steel. - –
The

*PD-Gr*model predicts a completely wrong crack pattern, in which neither the distributed bending cracks nor the major diagonal shear crack appears. Instead, an inelastic zone initiates under the loading pad and spreads out when the load increases. - –
The

*PDba-Gr*model predicts the formation of a narrow array of bending cracks and the failure then happens due to a major vertical crack instead of a diagonal one. Both are incorrect.

Figure 18(a), drawn in relative coordinates, shows again for convenience of comparison the diagonal crack locations and shapes, marked by discrete points, as obtained at *P*_{max} in the tests of geometrically scaled beams [97]. The points coincide remarkably well, as they should. This confirms the geometric similarity of the failure patterns.

- –
Both CB models are consistent in the location and inclination of the diagonal crack at the peak load, and match the experimental size effect well. However, for

*CB-Gr*, the inclination angle for beam depth*D*= 360 mm is incorrect. It is lower than for the other two sizes*D*and the location of the inclined crack is incorrectly shifted toward midspan (Figs. 18(b)–18(d)). - –
The load–displacement curves resulting from the

*bPF*model do not show any peak load; rather the load continues to increase without failure. Therefore, to compare and interpret the damage patterns and strengths, the peak had to be taken as the values at certain displacements marked as “fictitious”, particularly at*δ*= 10, 20, 40 mm for*D*= 200, 400, 800 mm (Fig. 18(e)), - –
Somewhat lesser errors occur in

*PF-Wu*. But the lack of a major crack at the peak load still leads to gross overestimation of the beam shear strength (Fig. 18(b)). Also, the*PF-Wu*model shows, incorrectly, a strength increase with increasing beam size. The culprit is obvious in Fig. 18(f), where the peak load is accompanied by a diffused zone of debonding and an inclined damage zone. The volume of this zone grows as*D*increases, which dissipates more energy. - –
Even greater problems are found in the

*bPD*model (Fig. 19(b)). The cracking patterns in beams of different sizes are inconsistent and incorrect. They evolve from inclined cracks to complete concrete-steel debonding at the peak load, and a major erosion of the material points. Even though the predicted strength values of scaled beams fall within the experimental scatter range, the mechanisms corresponding to these results are wrong (Figs. 19(a) and 19(b)). Also wrong is that the trend and curvature of the*bPD*size effect curve in Fig. 19(a) is opposite to that of the SEL. - –
Similarly, the

*PD-Gr*model does not produce a consistent trend of*σ*_{N}versus structure size*D*(Figs. 19(c)). This is in turn due to the inconsistency in the crack patterns predicted for various sizes. For example, while the bending damage dominates the failure mechanism for the size of*D*= 200 mm, the failure of beams with*D*= 400 and*D*= 800 mm shows a local accumulation of damage under the loading pads and nearby. - –
The data in

*PDba-Gr*model follow a trend opposite to other PD models, i.e., the nominal strengths are significantly lower than those observed experimentally. However, the cracking patterns are also inconsistent among beam sizes. While the failure mechanism at the small size is caused by localized bending cracks, the failure of beams with*D*= 400, 800 mm transits to damage accumulation under the loading pads, and bond breakage ultimately leads to a removal of material from this area.

The punching of RC slabs by a supporting column is a similar shear problem, except that the expansion of the punch zone is resisted by circumferential tension that provides radial confinement [98]. There exists an ACI-445 database of 440 test series conducted worldwide. The database confirmed that the nominal punching strength follows the SEL, Eq. (6), and this equation has been introduced also into the punching provisions of the ACI-318/2019 design code. As shown in Ref. [98], the CB model fits the data closely, while it is expected that PD and PF models would not, for the same reasons mentioned earlier. Because of the similarity to beam shear, the simulations of slab punching are here omitted.

### 5.9 Axial Double Punch of Cylinders.

This simple 1989 test [99] is an excellent check on the compression shear failure. The most important is the size effect when this test is scaled [99]. The double-punch loading of concrete cylinder produces a conical shear crack under each punch. Then, the penetration of the cone into the cylinder pushes the peripheral parts apart, and circumferential tension causes axial splitting cracks subjected to crack-parallel compression. The size effect is simpler and easier to evaluate than it is for the Brazilian split-cylinder test, for several reasons: (1) the scaling of the double-punch test is not complicated by changes of curvature of the Hertzian contact; (2) it does not use soft loading strips causing dissimilarity of Brazilian tests at increasing diameters; and (3) the plastic wedge zone forming under the loading strip in the Brazilian tests disturbs the geometrical similarity as the diameter increases.

- –
Both CB models predict a pattern agreeing well with the experiments of all sizes. The

*CB-Gr*model shows a more diffused damage than the*CB-M7*, but the main crack is correctly predicted by both. For larger*D*(from 76.2 to 1219 mm), the*CB-Gr*model overestimates the nominal strengths, which is due to excessive fracture energy. - –
Both the

*bPF*and*PF-Wu*models fail to reproduce the conical damage surface and show an incorrect failure mechanism. The*bPF*model leads to a cylindrical, rather than conical, damage surface, having the same diameter as the punch. The splitting starts, incorrectly, from this cylindrical surface. This also leads to excessive brittleness (i.e., sudden load drop). The*PF-Wu*shows only concentrated damage under the plate. For both PF models, the variation of the nominal strength across sizes follows, incorrectly, the power-law trend seen as a straight line in the log-log scale, similar to their Type 2 size effect behaviors for the crack of opening mode; see Fig. 6 (note that if data were plotted in the linear scale, such an obvious deviation would likely be ignored as the data fall within the scatter band). - –
All PD models either overestimate or underestimate the nominal strengths for every specimen size

*D*. For the case of*bPD*and*PDba-Gr*, not only they underestimate the fracture energy but also they show a wrong curvature of the type 2 size effect plot. These errors can be explained by Figs. 20(h)–20(j). The excessively brittle*bPD*model releases the strain energy abruptly by the removal of a near-cylindrical volume of material points, which ultimately leads to the splitting of the cylinder. The*PDba-Gr*, on the other hand, captures the conical surface under the loading pads. However, these surfaces localize to form, incorrectly, two elongated cones connected at apices while splitting cracks connecting the cones and separating cylinder parts are missing. - –
The increasing trend of the size effect curve produced by

*PD-Gr*(Fig. 20) is physically impossible. This stems from the increasing volume of damage predicted by this model as*D*increases. Compared with*CB-Gr*and*PDba-Gr*, this model shows a damage pattern resembling conical cracks, but the damage is spread out too widely.

### 5.10 Comparisons with Gap Test.

Finally the new test, the gap test, which has already been discussed in Fig. 1, and whose successful simulations by *CB-M7* was presented in Figs. 1(e) and 1(f). Due to differences in accuracy of various numerical models, the size effect method, when applied to different models, yields for the same material different values of fracture energy *G*_{f}_{0} at *σ*_{xx} = 0; see Fig. 21(a). Therefore, to facilitate comparisons, ratios *G*_{f}/*G*_{f0} are used as the coordinate in the diagram of Fig. 21(b).

- –
Both

*CB-M7*and*CB-Gr*follow the trend of gap test, which is an increasing*G*_{f}through low to medium*σ*_{xx}and weakening at high*σ*_{xx}. Quantitatively, though, the*CB-M7*captures the variation of fracture energy*G*_{f}with*σ*_{xx}better than*CB-Gr*(Fig. 21). - –
The PD models cannot predict the effect of crack parallel stress on

*G*_{f}at all. Furthermore, due to their excessive brittleness, the*bPD*incorrectly shows a monotonically decreasing*G*_{f}, explained by sudden bond breakages upon reaching their compression strength limit. This model predicts a premature failure (marked by ×) much before reaching the material compression strength*σ*_{c}. Among the PD models,*PDba-Gr*shows the least error and captures at least the increasing trend of*G*_{f}as*σ*_{xx}increases. - –
As

*σ*_{xx}increases, the*bPF*model shows incorrectly no change in*G*_{f}, which confirms it is a LEFM model. The formulation of*PF-Wu*, however, generates an FPZ with nonzero width that can interact with*σ*_{xx}. Yet this model gives an insufficient monotonic growth of*G*_{f}, which terminates prematurely (indicated by another × point) because, in using the size effect method, the gap test fails by compression before the gaps close. This is similar to what happens in the compression test of 2D confined specimens in which the biaxial stress state actually weakens the calculated response.

### 5.11 The Lack of Objectivity of Some Phase-Field and Peridynamic Models.

Structural analysis must be physically objective. This means that its results must be independent of human choice, particularly of the chosen quantities such as the mesh size (except for numerical errors that converge to zero with model refinement). When, however, the strain softening was introduced into the FE analysis, the mesh sensitivity of strain softening was found to be unobjective. The spurious mesh sensitivity was the initial reason for postulating the necessity of material characteristic length [100], in strain-softening material. It was justified by continuum smoothing of a material with microstructural heterogeneity and, in CB and cohesive crack models, by the interaction of fracture energy *G*_{f} with material strength *f*_{t}, two physical characteristics of different dimensions, N/m and N/m^{2}. When the stress attains the strength limit, a crack can form and propagate, but it will propagate only if a sufficient energy required by the fracture process zone is also supplied. Also, the results of the gap tests indicate the existence of at least two independent material lengths that determine the fracture energy *G*_{f} in quasibrittle materials [101].

In the *bPF* model, however, the crack is a line and its growth is fully characterized by *G*_{f}. The phase-field width *w*_{0} is an extra quantity with no physical justification. Its purpose is purely computational—to anchor the crack in the mesh without directional bias (note that *w*_{0} has nothing to do with Irwin’s material characteristic length $EGf/ft2$, characterizing the cohesive crack and CB models). However, the predictions of PF models are found to depend on the choice of ratio *w*_{0}/*h* (arbitrarily taken between 3 and 6). This is unobjective. For instance, in the *bPF* prediction of an unnotched 3PB specimen strength with *w*_{0} = 10 mm (or *w*_{0}/*h* ≈ 3), the peak load in Fig. 6 is different from that in Fig. 22(a), in which it is predicted with *w*_{0} = 20 mm (or *w*_{0}/*h* ≈ 6).

The sensitivity to the choice of *w*_{0} is also a serious conceptual problem of many other PF models, except that the *G*_{f0} value in the *PF-Wu* model for *σ*_{xx} = 0 is insensitive to *w*_{0} [102]. However, if *σ*_{xx} ≠ 0, then *PF-Wu*, too, is sensitive to *w*_{0} choice (note that, in Fig. 21(b), *PF-Wu* fails to capture the regime of weakening *G*_{f}). We posit that the modification in the phase-field functions (see Appendix C) inappropriately generates a finite-width FPZ. However, the FPZ expands in size as *w*_{0} increases at *σ*_{xx} ≠ 0, which leads to an excessive dissipated energy (Fig. 22(b)).

The crack band model overcomes the objectivity problem (except for a small numerical error) by adjusting the postpeak softening based on the element size [26,27,71]. For example, Grassl et al. made an adjustment to the softening slope corresponding to different element sizes in the CDPM2 (*CB-Gr*) model [60]. In addition, there is no objectivity problem if the element size in the damage zone is kept the same when the structure size is increased. This way is more accurate and is the preferred approach despite its computational burden.

Peridynamics, too, has an objectivity problem with respect to the choice of horizon radius. Similar to *PF-Wu*, PD models produce an FPZ scaled with the choice of the horizon size. Despite the inherent adjustment within the Grassl et al.’s model, which should suppress the spurious sensitivity with respect to different horizon radii in a way similar to the *CB-Gr* model, the *PDba-Gr* model still generates a larger FPZ width and predicts a higher peak load of the notched specimen, as shown in Fig. 22(c). This also applies to *PD-Gr*.

## 6 Discussion and Summary

### 6.1 Peridynamic Dilemma at Boundaries, Distorting the Size Effect.

The present comparisons confirm that the PD theory cannot be cured by introducing a better damage constitutive model. The problem is deeper.

As already pointed out, the need to cut off the potential or central force interactions in proximity to the boundary modifies the properties of the boundary unrealistically, due to sparser horizon connections in the boundary layer of a thickness equal to the horizon radius, *δ*_{h} (Fig. 5(a)). The cutoffs make this layer weaker and softer. To compensate for it at least partly, the properties within the horizons crossing the structure boundary or the crack face have usually been adjusted in some way, as already mentioned while discussing Fig. 5. The toughest problem is a horizon crossing the FPZ (or the cohesive zone) in front of a crack. In such a horizon, the cross-interactions must be weakened but not eliminated (Fig. 5(c)). The weakening should be modest near the front of the cohesive zone, and strong in its tail, close to the open crack tip. The boundary and crack line corrections of the horizon have different effects on different kinds of fracture. They are particularly disruptive for the size effect.

### 6.2 Comparison of Performance of 7 Models in Tests of 11 Types.

The results of the foregoing comparisons are summarized in Table 1. The performance of the *CB-M7* model is good (labeled OK) in all the 11 tests considered. The performances of *bPF*, *PF-Wu*, and *bPD* models range from poor to wrong in all the 11 tests except one, the concentrated shear test, in which their performances are satisfactory (or fair), which means less than good. The *PF-Wu* model can somewhat capture the type 1 size effect, but poorly. Among the PD models, *PDba-Gr* manages to capture more experimental observations than the others, but in 7 out of 11 tests, it is still incorrect. Notably, the *PD-Gr* model fails in all the tests, which implies a deeper conceptual problem for all PD models.

Test type | CB-M7 | CB-Gr | PF | PF-Wu | PD | PD-Gr | PDba-Gr |
---|---|---|---|---|---|---|---|

1. Size effect, type 1 | OK | OK | Wrong | Poor | Wrong | Wrong | Wrong |

2. Size effect, type 2 | OK | OK | Wrong | Wrong | Wrong | Wrong | Fair |

3. Concentrated shear fracture | OK | OK | Fair | OK | Fair | Wrong | OK |

4. Compression-torsion fracture (mode III) | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

5. Uniaxial compression | OK | Fair | Wrong | Wrong | Wrong | Wrong | Fair |

6. Confined compression of slab | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

7. Confined compression of cylinder | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

8. Vertex effect | OK | OK | Wrong | Wrong | Wrong | Wrong | OK |

9. Diagonal shear failure of RC beams | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

10. Double punch | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

11. Gap test | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

Test type | CB-M7 | CB-Gr | PF | PF-Wu | PD | PD-Gr | PDba-Gr |
---|---|---|---|---|---|---|---|

1. Size effect, type 1 | OK | OK | Wrong | Poor | Wrong | Wrong | Wrong |

2. Size effect, type 2 | OK | OK | Wrong | Wrong | Wrong | Wrong | Fair |

3. Concentrated shear fracture | OK | OK | Fair | OK | Fair | Wrong | OK |

4. Compression-torsion fracture (mode III) | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

5. Uniaxial compression | OK | Fair | Wrong | Wrong | Wrong | Wrong | Fair |

6. Confined compression of slab | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

7. Confined compression of cylinder | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

8. Vertex effect | OK | OK | Wrong | Wrong | Wrong | Wrong | OK |

9. Diagonal shear failure of RC beams | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

10. Double punch | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

11. Gap test | OK | Fair | Wrong | Wrong | Wrong | Wrong | Wrong |

### 6.3 Poor Performance: Where Does It Matter in Practice?.

Shear failure of RC beams, slabs, footings, prestressed or not

Horizontal shear of columns, shear walls, tall pylons in earthquake

Prestressed containment or pressure vessel failures

Safety of dams

Penetration of projectiles into concrete walls, exit speed

Safety of anchors in rock or concrete

Failure of composite beams

Fracking, especially with poromechanical stress transfer to solid

Hydraulic fracturing for geothermal energy or CO

_{2}sequestrationSlow subcritical growth of crack systems in geology

Rupture and growth of earthquake faults

Sea ice sheet pushing on a fixed structure, icebreaking

Load-bearing capacity of floating sea ice

Fiber composite airframes—fuselage crack, wing box, wing shear, rudder

Tunnel collapse, rock burst in mine stopes and deep boreholes

Many cases of fracture of fatigued metals, bone, and other biomaterials

Stiff soils, rigid foams, printed materials, etc.

### 6.4 Crack-Parallel Stresses: Where Do They Matter in Practice?.

Shear failure of RC beams and slabs

Cracks in prestressed concrete

Longitudinal crack in pressurized aircraft fuselage

Crack in the casing of solid-propellant rockets

Shear crack in aircraft wing, wing box, rudder, and stabilizer

Fracking, especially with poromechanical stress transfer to solid

Sea ice sheet pushing on a fixed structure

Pressure vessel fractures

Fracture in arch dam or arch bridge abutments, in footings

Crush cans for automobile crashworthiness

Cracks caused by projectile impact

Cracks in inflatable shells

Most thermal cracks

Cracks in geology, in seismic events

Pullout fracture of anchors from rock or concrete

Fuselage cracks, shear cracks in fiber composite wings, wing box, rudder

Shear crack in fiber composite wind turbine leafs

Particle comminution in projectile impact

Fracture of fatigued plastic-hardening polycrystalline metals

Burst of mine stopes, borehole breakout, failure of tunnels, excavations

Fracture of bone, dental materials, other biomaterials, etc.

## 7 Closing Comments

The PF model is currently the best approach for simulating curved or branching sharp LEFM cracks, though *only* in situations with negligible crack parallel stresses, and with no inelastic behavior outside the crack line. Unfortunately, such situations are *rare* in practice.

To become viable in general situations, the PF models would need to be fundamentally modified to consider:

A fracture process zone (FPZ) of finite width and finite length;

An FPZ characterized by a sufficiently realistic (multiparameter) tensorial constitutive model with multiparameter damage (such as microplane model M7);

Inelastic behavior outside the crack line;

The lack of objectivity with respect to the choice of the phase-field band width

*w*_{0}. The modified model would have to represent sufficiently well the experimental stress-strain and failure data as well as the crack-parallel stress effect.

For PD to become viable as a general predictive model, one would need to:

enhance the microstructure with particle rotations;

introduce dilatant frictional inter-particle shear;

avoid a particle-skipping force potential;

find a way to implement realistic boundary conditions, including those at the crack faces and at the fracture process zone boundary;

overcome the lack of objectivity with respect to the choice of peridynamic horizon (as well as the wave dispersion problems identified in 2016) [16].

These improvements seem impossible without a fundamental change of concept, which would inevitably lead to something like the LDPM [47,103].

The crack band model with the M7 microplane constitutive model agrees satisfactorily with all the existing experimental evidence and is predictive. Nevertheless, a further improvement is desirable to make the band run through the finite element mesh in any direction with zero mesh bias.

## Footnotes

This article documents the comparisons with experiments presented in Northwestern University joint TAM/SPREE lecture on Feb. 9, 2022 (see www.youtube.com/watch?v=3aNiTC_igM4www.youtube.com/watch?v=3aNiTC_igM4), which in turn was an extension of an ASME/IMECE plenary lecture on Nov. 1, 2021 (www.youtube.com/watch?v=HvVOobdHUlwwww.youtube.com/watch?v=HvVOobdHUlw).

The coding of model M7, both explicit and implicit, can be downloaded from www.civil.northwestern.edu/people/bazantwww.civil.northwestern.edu/people/bazant. It also includes the extensions to fiber reinforced concrete (M7f), rate effect, and comminution under impact.

Prof. Bazilevs and Dr. Behzadinasab of Brown University are thanked for providing the c++ code ready to be used with Peridigm, which facilitated the peridynamic implementation of the CDPM2 model. Both *PD-Gr* and *PDba-Gr* can be downloaded from github.com/htn403github.com/htn403.

Thanks are due to Dr. Fei and Prof. Choo of The University of Hong Kong for providing the code.

## Acknowledgment

Partial financial support under NSF grant CMMI-202964 and ARO grant W911NF-19-1-003, both to Northwestern University, is gratefully acknowledged. Thanks are due to Yuri Bazilevs, Masoud Behzadinasab, Jinhyun Choo, and Fan Fei for sharing their computer codes, and also to Florin Bobaru, Stewart Silling, Jian-Ying Wu, and Junuthula N. Reddy for their critical comments on the aforementioned ASME-IMECE lecture, which induced us to include further model variants and make the present analysis more detailed. Dat Ha is thanked for performing some preliminary checks of the phase-field models.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

All the data and all the information that support the findings of this article are freely available.^{12}

## Appendices

### Appendix A: The Calibration Process

In view of their diverse features, each model in this study had to be calibrated differently. The finite element discretization (which serves as the geometry for all methods) began by selecting the element size. For CB models applied to experiments lacking size effect data (which are what matter most), the element size was initialized as the maximum aggregate size of concrete, *d*_{a}, and then updated between *d*_{a} and 3 *d*_{a} based on the optimum fitting of existing data. For PF models, a sufficient number of elements need to be included in the band width *w*_{0} to avoid mesh bias and ensure the convergence, which normally range from three to six elements. For PD models, the horizon size was chosen to be 2.5–3 times the element size. For fairness of comparison, the horizon radius was chosen to be 1/2 of the crack band width. This ensured that the same RVE volume would be involved in the energy dissipation and that the PD results would be as close as possible to CB results (note that the PD models behave strangely when the element sizes are spatially graded to vary between a chosen minimum and maximum; e.g., in the unnotched three-point bend specimens in Fig. 7, the crack tends to form where small elements transit to larger elements rather than at midspan).

In the next step, the model parameters have been calibrated to fit best the data for each experiment. In the case of size effect tests, the calibration was based on the peak loads and the load–displacement curve for the smallest specimen, and then adjustments were made for larger sizes. The size effect data, if available, were prioritized over the postpeak data for one size, as they are more unambiguous [75]. Whenever available, data of other tests for the same concrete were also considered; e.g., Hoover et al. [33,74] performed uniaxial compression of cylinders and prisms, Brazilian tests, size effect tests on beams with different notch-to-depth ratios, etc., on the same concrete with a high consistency. For experiments with more limited material data, qualitative characteristics like crack patterns were also taken into account.

One must also note that each model requires different amount of data. The *bPF* model uses only three parameters—*E*, *G*_{f} and Poisson ratio *ν*. Hence, it requires fewer data than the *CB-M7* model, which involves eight free parameters *E*, *ν*, *k*_{1}, … *k*_{6}, although default values are often used for the last 3 free parameters (also, some parameters are more relevant than others, e.g., *ν* is unimportant for the type 2 size effect tests).

Even though the parameters of *CB-Gr*, *PD-Gr*, *PDba-Gr*) were kept the same in most cases, *PD-Gr* and *PDba-Gr* yielded some odd results. While, for example, the same parameter set (p1) resulted in the same nominal strength of the smallest specimens for type 1 size effect, a discrepancy occurred for type 2 (Fig. 23). To make the results comparable, *G*_{f} and *f*_{t} had to be rescaled to give the same strength for the smallest-size specimens (p2).

### Appendix B: Comparisons of Bond-Based Peridynamics to Size Effect Tests

For reasons stated in Sec. 4, the present comparisons used the (ordinary and nonordinary) state-based PD. However, the recent results of the bond-based version are worth mentioning. Hobbs et al. [52] recently studied the bond-based model for concrete and was the first to present (in Nov. 2021) the PD simulations of size effect in fracture, using published test data on similar notched and unnotched three-point-bend specimens of various sizes. These authors compared the PD results with the maximum load data measured in [73], but did not present the standard log-log diagram of size effect. This diagram is here reproduced (Figs. 24(a) and 24(b)) from the maximum load values evident from Hobbs et al.’s plots of load–deflection curves for various specimen sizes, for both type 1 and type 2 size effect laws. As shown, Hobbs et al.’s simulations with the bond-based PD completely miss the experimentally observed transitional size effect, for both types. Further note that, similar to some cases in Fig. 6, type 2 again gives a straight line, which corresponds to a power law and its slope corresponds to an exponent different from −1/2. This value makes the model thermodynamically inadmissible [27,71,78]; see Appendix H.

Hobbs et al.’s type 1 simulations give no size effect. This error cannot be attributed to a lack of Weibull statistics because the specimens are of the three-point-bend type, in which the zone of maximum stress involves only one RVE. To observe the statistical size effect, there must be a large zone in which the crack can originate from many RVEs, at many random locations [27,71]. Such a zone is best created in uniform tension specimens and, in a limited way, also in four-point bend beams with a very long segment of constant bending moment.

### Appendix C: Phase-Field Modifications by Feng and Wu

**u**, φ]

^{T}that minimizes the integral of the energy function (C1) over the entire domain.

### Appendix D: Phase-Field Modifications by Fei and Choo

### Appendix E: The Mathematical Basis of Bond-Associated Peridynamics

**X**and

**X**′ has been proposed:

### Appendix F: The Mathematical Form of CDPM2 Model by Grassl et al.

*g*

_{p}and plastic loading function

*f*

_{p}are different. The growth of compressive and tensile damages is described as follows:

### Appendix G: Basic Features of Microplane Model

The classical plasticity and damage models use the stress and strain tensors with tensorial invariants and their constitutive law is defined by loading surfaces, normality rules, plastic-hardening parameter(s) and damage parameter(s). Such models are called tensorial. The microplane constitutive model also relates the stress and strain tensors. But the constitutive law is vectorial [26,27,31,71,105]. It is defined as a relation of the stress vectors to the strain vectors, the latter representing the projections of the strain tensor onto a generic plane of arbitrary orientation within the material, called the microplane. The use of vectors helps physical insight, intuitively reflecting the microcrack openings, compression splitting cracks, shear slip, and frictional dilatancy. The algorithm is explicit, delivering the stress from strain, but it can be an advantage that on the microplane level it is easily converted to an implicit algorithm delivering the tangential stiffness tensor on the continuum level [53].

M7 is the latest in a series of progressively improved models M1–M7. In M7, as well as M3 and M4, all the inelastic behavior is characterized by the stress–strain boundaries. The return (or drop) to a boundary is in each loading step done at a constant strain (which amounts to a special case of the radial return algorithm for the tensorial models). The drop to the boundary on the microplane level suffices to guarantee nonnegativeness of energy dissipation. There is no need for nonassociated plasticity violating the normality rule. Unequal friction and dilatancy angles on the continuum level are reproduced automatically, with no possibility of negative energy dissipation (the nonassociatedness is a nonissue for microplane model).

The stress tensor is obtained variationally according to the principle of virtual work, by properly weighted integration over all spatial orientations. The integration is carried out numerically over a hemispherical surface according to one of the optimal Gaussian formulas. It amounts to weighted summation over a set of discrete microplanes (for good accuracy, their number per hemisphere must be at least 21 but typically is 37).

The boundaries of negative slope define the evolution of damage, and the horizontal ones define plasticity. There is no need for hardening plasticity, since it is automatically generated by interaction of the microplanes. Also, there is no use for a separate damage parameter per se, and there exists no single damage variable (they are many). There are one normal component and two shear components on each microplane. By using 37 microplanes per hemisphere, one has effectively 37 × 3 = 111 possible damage and plasticity sources. In effect, this is a sort of analog of multisurface plasticity with 111 independent loading surfaces, which are, however, vectorial rather than tensorial. This feature explains why the vertex effect is endemic and why it can occur at any point of the 9D space of stress tensor components. Beginning with Koiter [106,107] and Phillips et al. [108], it has been widely acknowledged that multisurface plasticity is more realistic, but 111 (even 10) plastic loading surfaces in the tensorial, rather than vectorial, space would be virtually intractable, both for model development and computer modeling (an omnipresent vertex effect nevertheless exists also in the tensorial endochronic theory [109] which, however has other limitations).

The fact that all the response within the boundaries is treated as elastic is a significant simplification. Curved rising stress–strain response is nevertheless automatically reproduced, thanks to different boundaries kicking in gradually in subsequent loading steps. The damage is generated on softening stress–strain boundaries at some microplanes but not at others, which is what creates the strong path dependence of the constitutive law at continuum level. The generation of damage on a microplane does not proceed monotonically. The damaged material can even stiffen when, e.g., hydrostatic pressure or transverse compression gets superposed on damage in shear.

During unloading and reloading, different microplanes become active. This is what reproduces the Bauschinger effect and the correct response under load cycles [110]. Even though the plentiful test data on cyclic and fatigue fracture are not covered here, the ability to capture the opening and closing of microcracks and sequential microplane activation in a material experiencing multiple loading cycles is also an important capability of model M7. The M7 crack band model [25,110] has been shown to possess such a capability. It can predict the behavior up to several thousand load cycles. The explicit description of damage on each microplane helps.

One feature that simplifies model formulation is that one-to-one relations between the conjugate stress and strain components on the microplanes are sufficient. There appears to be no need for cross-interactions such as the Poisson effect on the microplanes. The reason is that such cross effects are automatically generated by the interaction of microplanes of different orientations.

For a finite crack front width, the direction of crack growth is significantly affected by the triaxial stress state at the front. This effect plays a role in the curved propagation of a crack band in a compact tension specimen in which a hole has been drilled on the side of the crack extension line. As mentioned in Secs. 3 and 4, both PF and PD [23,111] predict a smoothly curved crack running into the hole. Despite the ruggedness imposed by the mesh, the M7 crack band model, too, predicts a similar curved path into the hole (Fig. 25). Note that, even if the band were not rugged, the path could not be identical. The reason is that the aforementioned triaxial stress effect, influencing the crack path direction, cannot be captured correctly by the PD and PF models. Therefore, and also because of missing the crack-parallel compression, the PD and PF crack paths cannot be quite realistic. This triaxiality effect in the FPZ is best manifested by sideways crack propagation in orthotropic fiber composites, making a sharp angle on the crack path. Such as path is correctly predicted correctly by the crack band model [112] but not by LEFM, PF, and PD models (nor by the cohesive crack model).

What matters for the good performance of CB-M7 is that this model implies three independent material characteristic lengths. The M7 correctly captures the fact [27,71] that the tensile postpeak softening curve begins with a steep softening slope and ends with a very long tail. The area under the initial slope is called the initial fracture energy *G*_{f}. The total area including the whole tail is called the total fracture energy *G*_{F}. The corresponding Irwin’s lengths, which are $lf=EGf/ft2$ and $lF=EGF/ft2$ matter mainly for the length of the FPZ, while the *w*_{0} is a material length characterizing the FPZ width, considered as the crack band width. The ability to capture all the three material characteristic lengths is an advantage of the crack band model compared to others.

Normally, the microplane model for concrete is calibrated by adjusting its three scaling parameters to fit the tensile strength and Young’s modulus. Four more parameters can be easily adjusted to fit the confined compression data. The initial fracture energy *G*_{f} is obtained through the fitting of the maximum loads for scaled notched fracture specimens using the size effect law (if such data exist). If *G*_{f} is specified, the scaling parameters of M7 are best adjusted so as to fit the size effect law (SEL) that corresponds to the *G*_{f} value [27,71]. Alternatively, a given *G*_{f} can be matched as the area under the initial tangent of the postpeak softening load-displacement curve (corrected for dissipation away from the FPZ, if any). The total fracture energy *G*_{F} is obtained as *G*_{f} times the ratio of the total areas under the total load–displacement curve and under its initial tangent.

### Appendix H: Why the Slope of a Straight-Line Size Effect in the Log-Log Plot Cannot Differ From −1/2

*C*of radius

*r*containing the FPZ:

*u*

_{i}∝

*r*

^{λ}. Hence,

*ε*

_{ij}∝

*r*

^{λ−1}and

*σ*

_{ij}∝

*r*

^{λ−1}. Therefore,

*J*∝

*r*

^{2λ−1}. Since

*J*must be finite, the only admissible value of Re

*λ*is 1/2. The scaling law must, therefore, be

*σ*

_{N}∝

*D*

^{−1/2}. Hence, if the log-log plot of structural strength versus size shows the slope of −1/

*η*with

*η*> 2, it means that

*σ*

_{ij}∝

*r*

^{−1/η}and

*J*∝

*r*

^{1/2−1/η}→ 0, which is impossible.