Critical Comparison of Phase-Field, Peridynamics, and Crack Band Model M7 in Light of Gap Test and Classical Fracture Tests

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

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³ [26,27]. Our study focuses mainly on concrete, not only because of the dominance of this material in national economies and its huge CO₂ 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₀ 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₀ from stabilized postpeak softening may be ambiguous [33]. A very rough estimate of l₀ is one to three times the inhomogeneity size of the material.

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.

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 $φ∝e±x/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₀ 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₀ spanning over a sufficient number of finite elements. Thus, w₀ 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.

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.

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:

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

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

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

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

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

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.⁴ For scaled specimens, the size of elements in the damage zone, which represents a material characteristic length, l₀, 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].⁵

• – 

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

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⁷), 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.⁸

• – 

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

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

• – 

  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).¹¹ 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₀ 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:

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

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

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

4. In terms of nation-wide usage, financial outlays, and the associated CO₂ emissions, concrete outstrips every other material by far (in terms of CO₂ 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.

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.

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.

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

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.

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 $ε=−0.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.

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.

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.

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.

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 $ε>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].

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′$⁠. 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.

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.

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

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². 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₀ 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₀ 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₀/h (arbitrarily taken between 3 and 6). This is unobjective. For instance, in the bPF prediction of an unnotched 3PB specimen strength with w₀ = 10 mm (or w₀/h ≈ 3), the peak load in Fig. 6 is different from that in Fig. 22(a), in which it is predicted with w₀ = 20 mm (or w₀/h ≈ 6).

The sensitivity to the choice of w₀ 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₀ [102]. However, if σ_xx ≠ 0, then PF-Wu, too, is sensitive to w₀ 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₀ 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.

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.

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.

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

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

3. Prestressed containment or pressure vessel failures

4. Safety of dams

5. Penetration of projectiles into concrete walls, exit speed

6. Safety of anchors in rock or concrete

7. Failure of composite beams

8. Fracking, especially with poromechanical stress transfer to solid

9. Hydraulic fracturing for geothermal energy or CO₂ sequestration

10. Slow subcritical growth of crack systems in geology

11. Rupture and growth of earthquake faults

12. Sea ice sheet pushing on a fixed structure, icebreaking

13. Load-bearing capacity of floating sea ice

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

15. Tunnel collapse, rock burst in mine stopes and deep boreholes

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

17. Stiff soils, rigid foams, printed materials, etc.

1. Shear failure of RC beams and slabs

2. Cracks in prestressed concrete

3. Longitudinal crack in pressurized aircraft fuselage

4. Crack in the casing of solid-propellant rockets

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

6. Fracking, especially with poromechanical stress transfer to solid

7. Sea ice sheet pushing on a fixed structure

8. Pressure vessel fractures

9. Fracture in arch dam or arch bridge abutments, in footings

10. Crush cans for automobile crashworthiness

11. Cracks caused by projectile impact

12. Cracks in inflatable shells

13. Most thermal cracks

14. Cracks in geology, in seismic events

15. Pullout fracture of anchors from rock or concrete

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

17. Shear crack in fiber composite wind turbine leafs

18. Particle comminution in projectile impact

19. Fracture of fatigued plastic-hardening polycrystalline metals

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

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

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:

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

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

3. Inelastic behavior outside the crack line;

4. The lack of objectivity with respect to the choice of the phase-field band width w₀. 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:

1. enhance the microstructure with particle rotations;

2. introduce dilatant frictional inter-particle shear;

3. avoid a particle-skipping force potential;

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

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

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.

There are no conflicts of interest.

All the data and all the information that support the findings of this article are freely available.¹²

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₀ 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₁, … k₆, 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).

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.

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