## Abstract

The crack band model, which was shown to provide a superior computational representation of fracture of quasibrittle materials (in this journal, May 2022), still suffers from three limitations: (1) The material damage is forced to be uniform across a one-element wide band because of unrestricted strain localization instability; (2) the width of the fracture process zone is fixed as the width of a single element; and (3) cracks inclined to rectangular mesh lines are represented by a rough zig-zag damage band. Presented is a generalization that overcomes all three, by enforcing a variable multi-element width of the crack band front controlled by a material characteristic length $l0$. This is achieved by introducing a homogenized localization energy density that increases, after a certain threshold, as a function of an invariant of the third-order tensor of second gradient of the displacement vector, called the *sprain tensor***η**, representing (in isotropic materials) the magnitude of its Laplacian (not expressible as a strain-gradient tensor). The continuum free energy density must be augmented by additional sprain energy $\Phi (l0\eta )$, which affects only the postpeak softening damage. In finite element discretization, the localization resistance is effected by applying triplets of self-equilibrated in-plane nodal forces, which follow as partial derivatives of $\Phi (l0\eta )$. The force triplets enforce a variable multi-element crack band width. The damage distribution across the fracture process zone is non-uniform but smoothed. The standard boundary conditions of the finite element method apply. Numerical simulations document that the crack band propagates through regular rectangular meshes with virtually no directional bias.

## 1 Introduction and Basic Idea

The fracture process zone (FPZ) at the front of propagating crack always has a finite width, and a distinct line crack appears only behind the crack front. The importance of a finite front width was recently highlighted by the gap test [1–3], which revealed an enormous effect of crack-parallel stresses on the material fracture energy. This effect cannot be realistically captured by the line crack models, such as the linear elastic fracture mechanics, the volumetric damage (Mazars) model [4], cohesive crack models [5–9], and phase-field models [10–14] in their current concept. Remedy might be sought in making the fracture energy dependent on crack-parallel stress, but this is insufficient since the dependence is strongly history dependent. It requires modeling fracture as a band of triaxial damage with a front of finite width. This is the hallmark of the crack band model (CBM) [15–18], for which a realistic damage constitutive model, such as the microplane model M7 for concrete or shale [19–22], is required (M7 subroutine can be freely downloaded from Ref. [23]).

An important advantage of the CBM [15,16,21] over other existing fracture models with distributed damage has been that it can describe the boundary and crack face conditions correctly because these conditions are the same as those of the finite element method. A superior performance of CBM compared to peridynamics and phase-field models was recently documented by extensive comparisons to 11 *distinctive* fracture tests important for practical applications [22,24]. The CBM, however, has three limitations:

The main one is that the single-element width of the band requires the damage to be uniform across the band.

The second is that the width of the FPZ cannot vary, being fixed to a single-element width.

The third is that, in a regular rectangular mesh, the single-element width leads to a zig-zag crack band when the propagation direction forms an angle with the mesh lines (although the overall propagation path, whose direction is controlled by the maximization of energy release rate, is still approximately correct in large specimens).

In contrast to the strongly nonlocal integral-type models [25] and to the weakly nonlocal gradient-type models (as well as to peridynamics), the CBM has two advantages:

The main one is that the boundary conditions are clear and physically justified. This is particularly important for the crack faces and the interface of the FPZ.

Another advantage is the tensorial character of material damage in the crack band, which is essential for taking into account the enormous effect of the nonsingular crack-parallel stresses on the material fracture energy and on the FPZ width.

The bias of propagation direction in the CBM is a significant problem only for a regular rectagular mesh. Random meshes, including the Voronoi mesh, and even meshes of squares consisting of four triangles, largely overcome this bias; see the crack band simulation of a curved crack in Fig. 5(b) of [24], in which a crack band with rugged interfaces was shown to approximate the curved, experimentally documented, crack path quite well. The strain-gradient model for softening damage in concrete developed by Cusatis et al. [26,27] has recently been shown to provide a good continuum representation of the lattice discrete particle model (LDPM). The LDPM has been very successful in discrete simulation of fracture in the microstructure of quasibrittle materials such as concrete.

Proposed here is a smooth crack band model (sCBM) which overcomes all the three limitations. Its basic idea is that the crack front width must be controlled by the localization-resisting energy characterized by the third-order tensor, $\eta $, of the second gradient (i.e., gradient of gradient or second derivative) of the displacement vector. This tensor, and the energy density, Φ, associated with it, is a new kind of localization limiter, which needs to be distinguished from the strain-gradient tensor, ** μ**, and from the strain energy density, Ψ. Therefore, we will call it briefly the

*sprain tensor*, $\eta $, and

*sprain energy density*, Φ. Tensor $\eta $ will be used to limit the localization of damage strain and to force it to be distributed over a certain material characteristic length $l0$, just like in medicine where the sprain of a ligament is not a rupture but the damage that is spread over a certain length of the ligament. In the same spirit, the sprain tensor is here employed only for softening damage to ensure that there are several finite elements with varying degrees of damage across the crack band width, making the crack front “smooth.” The sprain tensor is used to limit excessive damage localization only after a certain threshold,

*C*, has been exceeded. The threshold is set so that no appreciable “sprain” could occur during elastic and inelastic-hardening deformations.

### 1.1 The Concept of Homogenization of Damage and Localization-Resisting Energy Governed by “Sprain” Tensor.

In standard continuum mechanics, the homogenized free energy density, $\Psi \xaf,$ is assumed to depend solely on the change of length of infinitesimal line segments of the material, Fig. 1(a). From this, it inevitably follows that $\Psi \xaf$ is a function of, and only of, the strain tensor, $\epsilon $.

But for softening damage, this is insufficient. To explain it, consider first, for simplicity, the uniaxial deformation of a statistically heterogeneous bar shown in Fig. 1(c), with axial coordinate *x*, and consider the change of gradient $u\u2032=du/dx$ from $u1\u2032$ to $u2\u2032$. In the case of damage in heterogeneous material, e.g., a particulate composite such as concrete, the change of in-line gradient characterizing material damage, $\Delta u\u2032=u2\u2032\u2212u1\u2032$ (pictured in Fig. 1(b) as a change of slope), cannot be point-wise. Rather, it must be spread out over a certain material characteristic length $l0$, because the material heterogeneity in a bar under tension causes the microscale damage to be distributed over a certain material characteristic length $l0$, i.e., does not allow the damage to localize immediately into a line across the bar thickness. Consequently, we must set Δ*u*′ = *l*_{0}*u*″. This is similar to ligament *sprain* in medicine, which is the term for a ligament damage distributed over a certain length of the ligament, to be distinguished from ligament rupture. Note also that *u*″ cannot be generalized to 3D as the gradient of strain rather than displacement. Indeed, the change Δ*u*_{x,y} = *l*_{0}*u*_{x,yx} of the *y*-derivative of the *x*-displacement over length *l*_{0} cannot be expressed as the change of shear strain.

It is thus clear that the density of energy, Φ, resisting the localization is not a function of the displacement gradient, *u*′. Instead, it must depend on its change, Δ*u*′. Furthermore, while in one dimension, *u*″ coincides with the strain-gradient $\epsilon \u2032$, it does not in two or three dimensions as discussed later.

*localization-resisting energy density*, which does not exist in the continuum mechanics yet. Its absence is acceptable in all situations except softening damage mechanics (or quasibrittle fracture mechanics). Tensor $\eta $ is a third-order tensor of Cartesian components

*η*

_{ijk}=

*u*

_{k,ij}. It represents the third-order tensor of the second gradient of displacement vector

*u*

_{k}(subscripts

*k*= 1, 2, 3 label the Cartesian coordinates

*x*

_{k}, and the summation rule is implied). For the sake of brevity, $\eta $ will be called the

*sprain*tensor (according to the previous comment on ligament sprain in medicine). This tensor characterizes the gradual change of displacement gradient in any direction. Like $\Psi (\epsilon )$, $\Phi (l0\eta )$ will later have to be expressed in two or three dimensions as a tensorial invariant (i.e., independent of coordinate rotations). We will see that $\eta $ behaves differently than the strain-gradient tensor

**, with components $\mu kij=\epsilon ij,k$ (where, for small strains, $\epsilon ij=(ui,j+uj,i)/2$ = linearized strain tensor).**

*μ*In continuum mechanics of statistically heterogeneous materials (or quasibrittle fracture mechanics), the continuum strain is defined by homogenization. It is now proposed that, in damage mechanics, two different kinds of homogenization must be distinguished:

*Stiffness-based homogenization:*So far, material homogenization has generally been based on the equivalence of material stiffness and has typically been obtained variationally on the basis of the principle of virtual work, enforcing equilibrium while heeding the kinematic compatibility requirements. These homogenization methods include, for example, Voigt [28], Reuss [29] and Hashin–Shtrickman bounds [30]; Eshelby [31], Hashin [32], Hill’s self-consistent [33], Mori–Tanaka [34], Eringen [35,36], Bažant [37,38], Christensen’s composite spheres [39,40], and Dvorak methods [41]. Accordingly, the standard continuum thermodynamics of heterogenous materials is based on the free energy density Ψ that is a function of the strain of the continuum homogenized by stiffness only.*Energy-based homogenization:*As already pointed out, in heterogeneous materials exhibiting distributed damage, the stiffness-based homogenization misses the part of the homogenized Helmholtz free energy density, Φ, which is due to the resistance of damage localization in a heterogeneous microstructure. This part cannot be obtained by stiffness homogenization, even though it may be obtained by minimization of the total Helmholtz free energy of the system. It represents an*energy-based*homogenization, which is here achieved via second derivatives of displacements.

*n*= 1, 2, 3, …), see Fig. 1(d). The average of strain energy density based on stiffness homogenization is $\Psi (\epsilon \xaf)=(E/2)\epsilon stiffness2$ (solid horizontal line in Fig. 1(d). But the averaging of energy density gives a higher value, $\Psi \xaf(\epsilon (x))=(1/nl)\u222b0nl(E/2)(\epsilon stiffness+asin(2\pi x/l))2dx=(E/2)\epsilon stiffness2+(E/4)a2,$ which is shown by the dashed line in Fig. 1(d). The extra energy density due to energy homogenization is given as follows:

The energy-based homogenization is characterized by the third-order tensor, $\eta $, of the second gradient of displacement vector, the *sprain tensor*. This tensor is necessary only for describing the excess energy contribution of damage localizations that are caused by heterogeneity and are limited by it. In contrast to the classical strain-gradient formulations initiated in 1909 by the Cosserats [35,37,42–48], the strain-gradient cannot, in the case of softening damage, serve as a kinematic variable for the standard homogenized continuum, i.e., for a continuum homogenized in terms of stiffness only.

The present energy-based homogenization is different from the highly successful strain-gradient theory of Gao et al. [47]. Their theory correctly captures the size effect caused by geometrically necessary dislocations in plastic-hardening metals as shown by Huang et al. [48] but is not intended to control softening damage localization.

## 2 Forces Generated by the Localization-Resisting Energy

### 2.1 Softening Damage, “Sprain,” and Rupture in One Dimension—Simple Illustration.

*x*

_{1}≡

*x*. Let

*u*(

*x*) be the axial displacement, and again let Φ be the excess energy density of the localization-resisting energy, which we call, for brevity, the “sprain” energy density. Physically, Φ represents the energy of localized microcracking and microslips in excess of what is captured by the continuum homogenized by stiffness. This energy begins to develop only after the magnitude of $l0u\u2033$ exceeds a certain dimensionless threshold, denoted as

*C*, which must be high enough not to be significantly exceeded during elastic (or inelastic-hardening) deformations. Beyond this threshold, we assume the resistance to increase with the change of displacement gradient linearly, as if the behavior were elastic. This means that the associated (isothermal) Helmholtz free energy density, the sprain energy, should increase beyond the threshold quadratically. Hence, at a point of coordinate

*x*, the sprain (or localization-resisting) energy density is:

*C*is the “sprain” threshold, and

*κ*is the localization-resisting stiffness (or “sprain” stiffness). Its dimension is that of stress (i.e., Pa);

*l*

_{0}is the material characteristic length (characterizing in 3D the effective width of fracture process zone, different from Irwin’s characteristic length); and 〈

*X*〉 = max(

*X*, 0) = Macauley brackets, which help to introduce the assumption of energy equivalence of positive and negative

*u*″ (although here we focus on tensile rupture, Eq. (3) is also valid for localization in compression if threshold

*C*has a different value).

What is the meaning of the force variable, *f*, that must be associated with the excess localization-resisting energy density Φ? One might be tempted to consider the derivative *γ* = ∂Φ/∂*u*″, which has the dimension of surface tension, N/m, or fracture energy. But neither surface tension nor fracture energy has any physical meaning in the present context since no surface, nor crack, has yet formed. So, to define *f*, Φ needs to be discretized first.

*x*

_{1}≡

*x*into line elements of size Δ

*x*=

*h*delimited by nodes

*r*= 1, 2, …,

*N*. In terms of nodal displacement

*u*

_{r}, the second displacement gradient, i.e., the second derivative at node

*r*, $[d2u(x)/dx2]r=u\u2033(xr)=ur\u2033$, is approximated by

Note that we reserve subscripts *r*, *s* for the numbering of mesh nodes, while subscripts *i*, *j*, *k, m, n* will be used as the subscripts of Cartesian coordinates.

*r*is:

Here, *b* and *t* are the width and thickness, respectively, of the axially loaded bar subjected to stress in axial direction *x*.

Here, sgn(*u*″) = sign function = 1 if *u*″ is positive and −1 if it is negative (threshold *C* is considered as always positive). Noting that $f0\xaf$ is a constant if *u*″ is a constant, we see that the nodal forces *f*_{r} must be increased as 1/*h* when the steps *h* are refined. Calculation of $f0\xaf$, with *u*″ given by Eq. (4), may be described as follows:

*u*″, the nodal forces sum up to zero. Further, note that lim

_{h→0}

*f*

_{r}= ∞ in 1D (this changes for 2D if the element sizes in both directions tend to zero, as demonstrated in Sec. 2.3).

No triplet can be centered at the end nodes *r* = 1 and *r* = *N*, since it would protrude beyond the boundary. This ensures the physically correct form of boundary conditions.

### 2.2 Degradation of Sprain (Localization-Resisting) Stiffness *κ*.

*κ*constant is sufficient to handle early postpeak damage, but not full softening. Under excessive normal strain $\epsilon =u\u2032$, the material degrades and

*κ*decreases. So we express

*κ*, the sprain (localization-resisting) stiffness, as follows:

*κ*

_{0}is a material constant and

*λ*is a stiffness parameter (representing an internal variable). The simplest realistic evolution law for

*λ*is expressed as follows:

*k*

_{d}is a dimensionless material constant. The 〈..〉 ensures that the damage cannot decrease. For monotonic $\epsilon (t)$ and infinitesimal $\Delta \epsilon \u2192d\epsilon $, this integrates to yield the bell-shaped curve:

*u*′

^{2}in the exponent, we ensure that positive and negative displacement gradients (strain) cause equal damage. However, the 1D bar considered is under tension, and thus, a compressive damage does not appear in this example. For concrete,

*C*would need to be set much larger in compression than in tension. We may begin with

*k*

_{d}≈ 3

*f*

_{t}/

*E*, which gives a rough estimation of the order of magnitude.

*t*is given as follows:

*λ*

_{r}must be omitted.

In step-by-step integration, a simultaneous calculation of new nodal values of *u*_{r} and *λ*_{r} in each loading (or time) step would lead to a nonlinear problem, even for dynamic relaxation. Therefore, the damage is better delayed by one step, Δ*t*. For diminishing Δ*t*, the delay error should converge to zero. We evaluate the nodal values *λ*_{r} at the end of the previous step from Eq. (14), and then use them as given values in the current step. This way we obtain a system of linear equations in each loading step.

### 2.3 Quasibrittle Fracture Propagation in Two or Three Dimensions.

In a planar sheet of thickness *t*, with Cartesian coordinates *x*_{1} = *x*, *x*_{2} = *y*, the displacement vector $u$ has Cartesian components *u*_{1} = *u*_{x} and *u*_{2} = *u*_{y}. To generalize the expression (3), we must replace *u*″ by a two-dimensional tensor that represents the changes of the gradient (or slopes) of the displacements as functions of *x* and *y*.

So, limiting attention to isotropic materials, we must generalize *u*″ by invariant measures of the curvature of the surfaces of the in-plane displacement vector components *u*_{x}(*x*, *y*) and *u*_{y}(*x*, *y*). In the case of material isotropy, such measures are the Laplacians. The Laplacian of displacement vector $u$ is again a vector and thus not an invariant. The Laplacians $\u22072ux=u1,jj$ and $\u22072uy=u2,jj$, which represent the double of the mean curvatures of functions *u*_{x}(*x*, *y*) and *u*_{y} (*x*, *y*), are vectors (while the mean curvature of each surface is an invariant, equal to one half of the Laplacian).

*C*must now be replaced by a vector,

*C*

_{i}, even though, in the case of isotropy, its components must be equal, i.e.,

*C*

_{1}=

*C*

_{2}=

*C*= constants = dimensionless thresholds. So, the proper generalization of Eq. (3) giving the sprain (or localization-resisting) energy density as a tensorial invariant is

*i*,

*j*,

*k*= 1, 2, 3,

**= (**

*x**x*,

*y*,

*z*) in three dimensions, this equation is valid only for the typical hill-type damage distribution for which the Laplacians $ui,jj$ are negative for all $i$; for limitations of this equation and generalizations (see Comment C1 at the end of this section regarding certain invariance questions and limitations), and

*i*,

*j*,

*k*= 1, 2,

**= (**

*x**x*,

*y*) in two dimensions. Specifically for two dimensions, Eq. (17) for the sprain energy can be written as follows:

*κ*are N/m

^{2}(or J/m

^{3}). Note again that, despite symmetry with respect to

*x*and

*y*, Φ(

*x*,

*y*) is

*not*a function of the strain-gradient tensor. Further note that neither the mixed components (i.e.,

*u*

_{x,yy}and

*u*

_{y,xx}) nor the Laplacian could be expressed in terms of the strain-gradient tensor.

*sprain tensor*$\eta $, which represents the third-order second-gradient tensor of the displacement vector and is defined as

**with components $\mu kij=\epsilon ij,k=$$(ui,j+uj,i),k/2$, identified in [44,50] (see also Eq. (2.5) in Ref. [46] and Ref. [51]). Here, Eq. (18) is an analog of one of them. Analogs of the four others are not used for capturing the effect of the change of gradient of the displacement vector schematically portrayed in Fig. 1(b).**

*μ*For discretization, we consider a planar sheet of thickness *t* to be subdivided by a rectangular mesh with uniform steps *h*_{x} and *h*_{y} in coordinate directions *x* and *y*, respectively, the nodes being numbered as *r*, *s* (*r* = 1, 2, 3, …, *s* = 1, 2, 3, …) in the *x* and *y* directions. To obtain the nodal forces, we introduce in Eq. (18) the second-order approximations of the second partial derivatives:

*r*and

*s*. The

*x*and

*y*components of the localization-resisting forces at the five nodes centered around the node (

*r*,

*s*) are obtained as follows:

*u*

_{x}yields,

*u*

_{y}, the calculation is similar. Introducing dimensionless Laplacians,

*r*,

*s*):

*x*and

*y*(see Fig. 3).

The *scaling rule* of nodal forces in 2D is given by Eqs. (29)–(34). When both steps *h*_{x} and *h*_{y} are scaled in proportion, the nodal forces *f*_{x} and *f*_{y} do not change with the step size. When *h*_{x} is refined while *h*_{y} does not change, the *f*_{x}_{r,s} scales, for small *h*_{x}/*h*_{y}, asymptotically in proportion to 1/*h*_{x}, the same as in the one-dimensional case, and asymptotically for large *h*_{x}/*h*_{y} in proportion to *h*_{x}. Similar scaling properties hold for other components.

*κ*(

*λ*) =

*λκ*

_{0}and the 2D evolution law for

*λ*is

*k*

_{d}is an empirical dimensionless material constant. The discrete form of Eq. (35) in increment Δ

*t*is expressed as follows:

Note that, unlike the compressive volumetric strain, the in-plane compressive area strain is appropriate for characterizing damage in the case of plane stress or low confinement in the *z*-direction because it allows free expansion in the third direction.

In the numerical algorithm, the area strain is again delayed by one step. This preserves the linearity of the equation system for nodal values *u*_{r,s} to be solved in each step Δ*t*.

In three-dimensional fracture problems, a realistic damage constitutive law such as the microplane model is required, and the foregoing equations need to be generalized to three dimensions. Equation (17) is applicable in three dimensions when *i*, *j*, *k* = 1, 2, 3.

*Comment C1: Limitations and Generalizations.* Applying the Macauley brackets to the two Laplacians is convenient but infringes, of course, on the tensorial invariance, if both terms of the product in Eq. (17) are not positive. Yet this seems to have little effect on the present simulations, probably because the Laplacian of $uy$ in the direction of propagation is generally much smaller than that of $ux$, and because both terms of the products in Eq. (17) are positive for the hill-type localization of damage profile across the crack band. To limit the positive curvature (or valley-type localization) at the edges of the band would require introduce a second sprain energy with signs different from Eq. (17), but practically this does not seem important. Also note that, in-the case of unloading and reloading, stiffness κ and threshold $C$ would have to be generalized to a more complex form, possibly involving convex programming aspects akin to Melan's loading–unloading (or Karush–Kuhn-Tucker) criteria in the theory of plasticity (see, e.g., Ref. [52] pp. 204 and 244). Note also that the case of anticlastic curvatures, in which $ui,xx$ and $ui,yy$ have opposite signs, is here dismissed as likely to be unimportant in practical situations.

### 2.4 Anisotropic Generalization and Sprain Tensor.

*α*

_{ij}. Equation (17) may then be simply generalized as follows:

*C*

_{k}are the thresholds for the orthotropy directions, and the sprain tensor

*η*

_{ijk}=

*u*

_{k,ij}is defined by Eq. (20) (we skip discussing the conversion of

*κ*to general anisotropy).

Material isotropy is the special case for *C*_{i} = *C* and for *α*_{ij} = *δ*_{ij} = Kronecker delta (unit tensor). The associated localization-resisting energy density, Φ(*x*, *y*), may again be succinctly called the *sprain energy*.

The strain-gradient tensor may be regarded as a symmetrization of the sprain tensor, $\eta $. Tensor *α*_{ij}*η*_{ijk} is a generalization of the vector of Laplacian $\u22072uk$.

It must be emphasized that the need for the sprain tensor *η*_{ijk}, along with its associated sprain energy Φ, is limited to postpeak softening damage. Threshold *C* (or *C*_{i}) must be set such that the sprain energy does not affect appreciably the elastic and inelastic-hardening deformations.

### 2.5 Effective Numerical Treatment of Boundary Conditions and Element Size Changes.

At the mesh boundary, the last force triplet for the boundary node must be omitted since it may not be allowed to protrude beyond the boundary, as this would cause the in-plane corrective forces at the nodes of the boundary element to be unbalanced and produce spurious damage. The standard boundary conditions of the finite element method apply. Alternatively, it is possible to superpose at the boundary node the doublet of in-plane corrective nodal forces (*f*, −*f*), representing one half of the triplet. This doublet is also self-balanced.

Often the damage does not extend to the boundary. Then one may single out a region should to be free of damage—a region whose boundary is expected to be in the elastic domain, and the original CBM or regular finite elements analysis is used in this region. In this region, the finite elements have the full size for a crack band of single-element width, and the damage development is blocked. This may be simply achieved by setting the threshold in this region to be large enough, or sprain (localization-resisting) stiffness *κ* to be small enough, or both.

The formulation would become more complicated if the element sizes were changed within the damage zone. Therefore, it is preferable to use a uniform mesh over the entire damage zone and make element size changes only in the elastic or hardening inelastic zone, i.e., outside the damage zone. Note that excessive damage may appear when small finite elements are used at boundary supports of small contact areas. At such supports (e.g., in simulations of the three-point bend test), larger finite elements are preferable.

### 2.6 Differential Equation of sCBM and Finite Difference Approach.

*pu*is the potential energy of the axial distributed load

*p*applied per unit volume of material, $u\u2032=du/dx$, and $u\u2033=d2u/dx2$. Let us restrict ourselves to the case that

*u*″ <−

*C*/

*l*

_{0}(i.e., to a hill-type profile of

*u*(

*x*)). Consider first that the sprain (localization-resisting) stiffness,

*κ*(

*λ*), is fixed. We set the variation of Eq. (39) to zero, so as to ensure equilibrium, and then integrate by parts twice, so as to eliminate the derivatives of

*δu*. We obtain the following:

*δ*Π must vanish for any

*δu*, the expression in parenthesis must vanish, too. Thus, we obtain the differential equation:

Interestingly, Eq. (43) is formally identical to the bending equation for a beam of bending stiffness $\kappa l02$, under distributed transverse load *q*(*x*).^{2} This is not surprising since bending stiffness is what prevents localization of curvature. Here, we consider no applied body forces, and so load *q*(*x*) = *σ*′(*x*).

*f*

_{r}, Eq. (43) may be rewritten as $uIV=(u\u2033)\u2033=(\sigma \u2032+p)/(\kappa l02)=q/(\kappa l02)$. Then, replacing (·)″ by its finite difference approximation at node

*r*, we have

*q*

_{r}is the

*q*value at node

*r*. Now, to relate

*u*

_{r}to the nodal force, we take a derivative of Eq. (3) and obtain (for

*u*″

_{r}< 0, hill):

*u*″

_{r}and similarly

*u*″

_{r−1},

*u*″

_{r+1}, from this equation, we obtain, upon substitution into Eq. (44) and after rearrangements:

This agrees with the scaling rule for the triplet of nodal forces, which we previously derived more directly by differentiation of the free energy density with respect to nodal displacements *u* instead of *u*″. The scaling rule for nodal forces is given by Eq. (7).

Equations (45) and (46) vary in time and from node to node. Thus, the fourth-order (flexure) differential equation would have to be solved repeatedly in each subsequent load or time-step Δ*t*, taking the *κ* value from the previous step.

The foregoing analysis can be extended to 2D as well as 3D. The total free energy of a structure of constant thickness *t* with 2D coordinates *x*, *y* may be written as follows:

*δ*Π in 2D and using the Gauss integral theorem, one obtains a fourth-order partial differential equation that is mathematically analogous to the plate bending equation.

## 3 Numerical Simulations

In the original CBM, the single-element width of the crack band is taken equal to the material characteristic length, *l*_{0}, and the strain, $\epsilon =u\u2032$, is uniform across the band. In the present sCBM, *l*_{0} is subdivided into several uniform-strain elements of width *h*. Figure 4 illustrates schematically the parabolas of displacement and the corresponding first gradient of displacement profiles meshed with coarse elements, Figs. 4(a) and 4(c), and fine elements, Figs. 4(b) and 4(d).

*C*(to be later recalibrated by experiments). We introduce a subdivision of

*l*

_{0}into many elements. In that case, the distribution of strain $\epsilon =u\u2032$ across the band width is smoothed to approach the bell-shaped function depicted in Fig. 4(d), which can be approximated by four parabolic arcs. Based on the properties of a parabola, the maximum slope of this bell-shaped function should be approximately equal to threshold

*C*, which means that

*u*″

_{max}=

*u*′

_{max}/(

*l*

_{0}/2) =

*C*. From this we get the estimate (a crude one):

For tension $u\u2032=\epsilon >0$, threshold *C* may be considered as constant. For compression, though, the realistic *C*-value might be different. For example, in concrete, the uniaxial compressive strength *f*_{c} is about ten times higher than the tensile strength *f*_{t}. Then the *C*-value for uniaxial compression would have to be much greater than for tension. With lateral confinement, or under triaxial compression, *C* may have to be much higher still. However, if the *C*-value is set too high, the softening behavior and fracture in compression might be excluded from modeling.

### 3.1 One-Dimensional Tension Test With a Constant Sprain (Localization-Resisting) Stiffness.

Consider now a one-dimensional bar of length *L* = 300 mm subjected to symmetric displacement boundary conditions in Fig. 5(a). The applied displacement *u*_{0} is increased linearly with time *t* and reaches *u*_{0max} = 4 × 10^{−4}*L* at *t*/*t*_{max} = 1.0 (Fig. 5(b) where the total time *t*_{max} = 1.0 s). The material follows a simple bilinear stress versus strain relation for an elastic-softening response (Fig. 5(c)). The elastic strain is given by Hooke’s law with Young’s modulus *E* = 25 GPa. To initiate localization in the desired place, the material in the segment *x* ∈ [150, 165] mm (1/20 of the total length) is assumed to have a slightly reduced tensile strength of 0.99 *f*_{t}, and the material in the rest of the bar is assumed to be elastic up to a tensile strength *f*_{t} = 4 MPa. The corresponding strain at the tensile strength limit is $\epsilon 0=ft/E$.

Simple linear softening, corresponding to fracture energy *G*_{f} = 0.04 N/mm, is assumed after the stress reaches the tensile strength. The strain after total softening is taken as $\epsilon c=2Gf/(l0ft)$. The material characteristic length *l*_{0} in Eq. (3) is taken to be the aggregate size *l*_{0} = 50 mm, and then, $\epsilon c/\epsilon 0=2.5$. The problem is static, but the explicit dynamic relaxation algorithm with artificial mass density *ρ* = 2.6 × 10^{−9} ton/mm^{3}, damping coefficient 10^{−3} ton/s, and time-step $\Delta t=0.8h/E/\rho $ is used for time integration (in which case the ratio of kinetic energy to strain energy is negligibly small, as required for quasi-static solution). The parameters for triplet forces in Eqs. (7) and (8) are material constants consisting of sprain (localization-resisting) stiffness *κ*_{0} = 10 MPa in Eq. (13), second gradient threshold *C* = 10^{−3}, bar width *b* = 1 mm, and bar thickness *t* = 1 mm.

First, we consider the case of constant sprain (localization-resisting) stiffness *κ* with *λ* = 1.0 in Eq. (13). We divide the bar into 6, 50, and 150 elements, which correspond to element sizes $h=50mm,6mm,and2mm$. The profiles of displacement, strain, the superposed triplet forces, and the second derivative of displacement of the bar at *t* = *t*_{max} are shown in Fig. 6 for element sizes of 50 mm (*h*/*l*_{0} = 1), Fig. 7 for an element size of 6 mm (*h*/*l*_{0} = 3/25), and Fig. 8 for an element size of 2 mm (*l*_{0}/*h* = 1/25), respectively. Note that within bar segments of uniform curvature, the superposed triplet forces cancel each other, and the damage then concentrates in one element, as shown in Figs. 6(a) and 6(b). Figure 6(c) shows the second derivative of displacement, which has both a convex and a concave characters. Figure 6(d) shows the superposed triplet forces obtained from *u*″ in Fig. 6(c) as an input, Eqs. (3) and (6). Figure 6(d) illustrates that superposition of two force triplets of opposite directions generates a force of larger magnitude. The deformation at the neighboring nodes is caused by the triplet forces due to the large *u*″ magnitudes at these nodes.

Figure 7 shows that, in contrast to Fig. 6, more elements are deformed to develop a “failure zone,” which represents the length over which the element strains are larger than elsewhere, while the element size is relatively small (Figs. 7(a) and 7(b)). The profile of the second derivative of displacement *u*″ in Fig. 7(c) is similar to that in Fig. 6(c) because more nodes lie in the transitional region between the nodes of maximum convex and maximum concave curvatures of *u*. Because of the symmetry of the deformation with respect to the center element, the superposed triplet forces are significant only at the places near the nodes with the largest *u*″ values. At the remaining nodes, the superposed forces between those of the largest *u*″ magnitude are significant, while the superposed triplet forces at nodes regarded as the “strain boundary” have negligible values, Fig. 6(d).

Figures 6–8 show that larger *l*_{0}/*h* induces damage distribution over more elements representing the softening “damage zone” (or sprain zone); see the subfigures (*a*), (*b*), and (*c*) in Figs. 6–8. Across the “sprain zones,” there is 1 element in Fig. 6(b) with *h*/*l*_{0}/ = 1; 7 elements in Fig. 7(b) with *h*/*l*_{0} = 3/25; and 25 elements in Fig. 8(b) with *h*/*l*_{0} = 1/25. The corresponding “sprain zone” lengths are 50 mm, 42 mm, and 50 mm for the bar meshed with the element sizes of 50 mm, 6 mm, and 2 mm, respectively. These “sprain zone” lengths approximate the characteristic length *l*_{0} = 50 mm.

Zoomed-in figures of the corresponding deformation profiles are plotted in Fig. 9 to show the distributions clearly. The profiles document that, for various element sizes, the “sprain zone” length and damage profiles for the same material parameters (the localization threshold and the resisting stiffness) are very similar (see Figs. 6–8). However, when the element size is decreased, the maximum superposed triplet force values increase, as expected (see Figs. 6(d), 7(d), and 8(d)). According to Eq. (7), when the element size tends to zero, the triplet forces in one dimension tend to infinity (not in 2D, though, if *h*_{y} ∝ *h*_{x}).

### 3.2 One-Dimensional Tension Test of a Bar With a Degrading Sprain (Localization-Resisting) Stiffness.

The resistance to localization should decrease with increasing damage and vanish when the material is totally broken. If a constant sprain stiffness *κ* is used for postpeak damage, there will be residual stresses in the bar; see Fig. 17 in Appendix B. Therefore, a constant sprain (localization-resisting stiffness) *κ* is usable only for near postpeak damage. Beyond that, *κ* must be decreased to cause the triplet forces to diminish as the damage progresses. Here, we show further simulations when the triplet force is diminished according to the stiffness parameter, *κ*(*λ*), in the evolution equation, Eq. (14).

We repeat the previous calculations except for applying the dimensionless material constant *k*_{d} = 2.0 × 10^{−3}. Figure 10(a) shows the reaction force *P* at the right end of the bar versus the displacement *u* at the right end with an element size 2 mm (*l*_{0}/*h* = 25). The first row of Fig. 11 shows the profiles of displacement, strain, and normalized sprain (localization-resisting) stiffness at the maximum reaction force, marked by a dot in Fig. 10(a). The dashed line in Fig. 11(b) represents the elastic strain corresponding to the tensile strength. Thus, only the elements at the center develop damage, while the remaining elements are still in the elastic range. The deformations are smaller near the central localization region due to force reduction in this region at the instant of maximum load. The sprain (localization-resisting) stiffness *κ* decreases with excessive strain, according to the evolution equation, Eq. (16), for the stiffness parameter *λ*; see Fig. 11(c).

The second row of Fig. 11 shows the corresponding profiles at *t* = *t*_{max} when the bar almost breaks. The deformations eventually localize signaling rupture (Figs. 11(d) and 11(e)). However, a small residual deformation still exists around the broken element, as seen in Fig. 11(e) and the zoomed-in picture in Fig. 10(b). The sprain (localization-resisting) stiffness *κ* of the nodes connecting the broken element at the center becomes zero, while the nearby nodes still exhibit small *κ* values due to triplet forces induced by excessive *u*″, while the *κ* values in the remaining elements are unchanged (Fig. 11(f)). Since the *κ* values of the nodes near the broken elements are diminished, both the triplet forces and the corresponding residual stresses have negligible values at *t*/*t*_{max} = 1 (see Fig. 18 in Appendix B). The results for the mesh with element size 6 mm are similar, showing only a slight quantitative difference.

### 3.3 Simulation of Three-Point Bend Fracture Test in 2D.

*D*, width 2

*L*+

*W*, unit thickness

*t*and in an in-plane stress state (Fig. 12). The aspect ratio

*L*/

*D*= 2.5, and

*W*denotes the notch width and

*a*

_{0}the initial notch length. For the explicit integration, the rate boundary conditions are as follows:

*a*)

*b*)

*c*)

The material parameters are Young’s modulus *E* = 30 GPa, fracture energy *G*_{f} = 0.05 N/mm, and tensile strength *f*_{t} = 4 MPa. The characteristic length *l*_{0} in Eq. (17) is taken to be *l*_{0} = 25 mm. The parameters for triplet forces in Eqs. (29)–(34) are the sprain (localization-resisting) stiffness *κ* = 200 MPa and the localization threshold *C* = 10^{−5}. For the simulated fracture to be quasi-static, we ensure with mass density *ρ*′ = 2.6 × 10^{−5} ton/mm^{3} that the ratio of kinetic energy to strain energy is less than 10^{−3} before and during crack propagation. The damping coefficient and the formula for correct time-step are the same as in the one-dimensional example. The finite elements are quadrilateral, with full Gaussian quadrature.

With the original CBM, the tests terminate when the postpeak load at the loading point drops to 0.8 of the maximum load value of the test. Figure 13 shows the contour of strain component $\epsilon xx$ calculated with original CBM using finite element meshes rotated by various inclination angles, while keeping the same element size within and near the expected damage area. The figure corresponds to postpeak load drop to 0.8 of the maximum load value.

Figures 13(a)–13(d) show the full-field results, and Figs. 13(e)–13(h) show the corresponding zoomed-in plots. The element sizes within and near of damage zone are the same—5 mm for all four meshes with inclination angles *α* = 0 deg, 15 deg, 30 deg, and 45 deg. Figures 13(e)–13(h) clearly illustrate the dependence of propagation path on the mesh inclination. Specifically, the crack propagates along the vertical mesh line when *α* = 0 deg in Fig. 13(e), and along the tilted mesh line when *α* = 15 deg in Fig. 13(f). In Fig. 13(g) for *α* = 30 deg, the crack propagates vertically, leans slightly to the right and the boundary of the crack band has a zig-zag shape. In Fig. 13(h) for *α* = 45 deg, the crack propagates along a vertical line because of the symmetry of the mesh.

The contours of stress component *σ*_{xx} are similar to the contours of strain component $\epsilon xx$ and are shown in Appendix B. Note that the *x* and *y* directions are defined by the fixed Cartesian system and do not follow the mesh lines when the mesh is rotated.

The deviation of the crack propagation path from the straight-ahead *y* direction is physically wrong. It demonstrates the bias of the original CBM in the case of a regular rectangular mesh. It is one of the three aforementioned arguments in support of sCBM. Note, however, that the directional bias of CBM essentially disappears in the case of random mesh, at the penalty of scatter.

With sCBM, the tests terminate when the prescribed displacement at loading point reaches *u* = 0.12 mm when *t* < *t*_{max}. Figure 14 shows the corresponding contours of strain component $\epsilon xx$ calculated with sCBM when *u* = 0.12 mm. The meshes used in Fig. 14 with the sCBM are the same as the meshes used in Fig. 13 with the original CBM. The zoomed-in plots in Figs. 14(e)–14(h) clearly illustrate that the strain distribution is independent of the mesh inclination, i.e., the contours of $\epsilon xx$ are similar while the finite element meshes used in the calculations have various inclination angles *α* = 0 deg, 15 deg, 30 deg, and 45 deg in Figs. 14(e), 14(f), 14(g), and 14(h), respectively.

In contrast to the one-element wide crack band generated with the original CBM, the sCBM generates a crack band with multiple elements across the band. The width of the crack band is proportional to the characteristic length *l*_{0} and can be numerically controlled by the values of sprain (localization-resisting) stiffness *κ* and threshold *C* for the triplet forces. The contours of stress component *σ*_{xx} normal to notch calculated for various mesh inclination angles are also similar and demonstrate the independence of stress distribution of the mesh inclination. Full-field contours are shown in Figs. 15(a)–15(d), and zoomed-in plots are shown in Figs. 15(e)–15(h).

Figure 16 shows the curves of load versus displacement of the loading point at the top of the beam. Figure 16(a) shows the responses obtained from calculations using original CBM with contours in Figs. 13 and 19. The four curves are indistinguishable before approaching the maximum load and differ at maximum load appreciably. The curve with inclination angle *α* = 0 deg has the lowest maximum load value.

As the mesh inclination angle increases from 0 deg to 30 deg, the maximum load value increases, while the crack propagation path changes from a smooth band to a zig-zag band, as the latter requires more energy to be dissipated within the same vertical crack length. The maximum load when *α* = 30 deg is larger than for the angle *α* = 45 deg for which the mesh becomes symmetric again. The increase in maximum load from *α* = 15 deg to *α* = 30 deg is more noticeable than the increase from *α* = 0 deg to *α* = 15 deg, and from *α* = 45 deg to *α* = 30 deg.

Figure 16(b) shows the responses calculated using the sCBM. The four curves are very close, with only a slight difference at the maximum load value. As expected, the difference between *α* = 0 deg and *α* = 45 deg is the least since both meshes are symmetric with respect to the centerline of the beam. For *α* = 15 deg and *α* = 30 deg, the maximum values are only slightly smaller. The maxima of the load versus displacement curves obtained with sCBM and original CBM do not differ by much, but the corresponding deflection is appreciably larger for sCBM.

## 4 Conclusions

The homogenization of heterogenous materials has so far been conducted on the basis of stiffness (with equilibrium typically imposed via virtual work). This is exemplified, for example, by Voigt [28], Reuss [29] and Hashin–Shtrickman [30] bounds; Eshelby [31], Hashin [32], Hill’s self-consistent [33], Mori–Tanaka [34], Eringen [35,36], Bažant [37,38], Christensen’s composite spheres [39,40], and Dvorak methods [41]. However, the extra energy density of damage strain localization was not captured by these methods. It requires an energy-based homogenization and matters only for softening damage.

In continuum mechanics of damage of quasibrittle materials, the total energy density, $\Psi \xaf$, ought to be the sum of (1) the usual energy density, Ψ, which is a function of the strain tensor and (2) an additional localization resisting (or

*sprain*) energy density Φ that is a function of the third-order tensor of the second gradient of displacement vector (for brevity called the*sprain tensor*) multiplied by the material characteristic length,*l*_{0}. If the material is isotropic, Φ becomes a function of the magnitude of the Laplacian of the displacement vector in excess of a threshold and cannot be expressed as a function of the strain-gradient tensor.Below a certain threshold

*C*, the*sprain*(localization-resisting) energy Φ vanishes and, above*C*, it is a function of the*sprain tensor*times*l*_{0}. During the postpeak softening, the corresponding localization-resisting (or*sprain*) stiffness,*κ*, needs to be reduced as the damage, characterized by volumetric strain, increases.Differentiating the total energy density of a finite element system with respect to the nodal displacements yields self-equilibrated nodal forces resisting localization during postpeak softening. In one dimension, these are self-equlibrated force triplets. In two dimensions, these are self-equlibrated couples of crossing in-plane force triplets.

Finite element calculations document that the sCBM can simulate propagation of crack band through a regular square mesh with no noticeable directional bias when the mesh is inclined.

The crack band front has a multi-element width and thus can represent a

*smoothed*damage profile. The band width depends on the characteristic length*l*_{0},*sprain*(or localization-resisting) stiffness*κ*, and threshold*C*. The band width varies with the overall stress state of a structure. This is important for the effect of crack-parallel stresses, recently revealed by the gap test.

## Footnote

For a solution by finite difference method, one would need to introduce the approximation *u*^{IV} ≈ (*u*_{r−2} − 4*u*_{r−1} + 6*u*_{r} − 4*u*_{r+1} + *u*_{r+2})/*h*^{4}, but this does not seem to be an effective approach.

## Acknowledgment

Partial financial support under NSF (Grant No. CMMI-202964) and ARO (Grant No. W911NF-19-1-003), both to Northwestern University, is gratefully acknowledged. Hoang T. Nguyen of Brown University (formerly Northwestern) is thanked for valuable discussions of computer implementation.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

### Why Other Invariants of Sprain Tensor Are Inapplicable

*η*

_{iik}

*η*

_{jjk}(also [51]). To explain why the other four independent invariants of $\eta $ are inapplicable, consider the energy densities Φ

_{a}, …, Φ

_{d}corresponding to the other four invariants of tensor $\eta $:

Equations (A1) and (A2) use a spread-out change of volumetric strain, $3l0\epsilon V,i$, which cannot capture the spread-out change of damage in Fig. 1(b). The lack of symmetry of Φ_{b} is, for energy, unacceptable. In Eqs. (A3) and (A4), the physical meaning of the third-order threshold tensor *C*_{ijk} is dubious (on the other hand, a vector threshold *C*_{i} in Eq. (17) is physically necessary in the case of orthotropy). The components such as (*u*_{1,2})_{,3} (*u*_{2,1})_{,3} in Φ_{d} have questionable interpretation in terms of work. For these reasons, no finite element models corresponding to these other four invariants have been attempted.

*Comment on Strain Gradient*: It has not been checked in detail whether another one of the invariants of strain gradient or their combination could serve as the energy potential whose derivatives with respect to the nodal displacements would yield forces resisting excessive localization. But this seems to lead to a much more complicated nodal force system. A situation where it could make a difference is the hypothetical case of saddle (anticlastic) surfaces of displacement components with zero Laplacians, in which case one would have opposite *u*-curvatures of equal magnitudes in two orthogonal directions. This would, however, suggest a tensile crack band crossing a compression fracture band at the tip, which hardly looks as a realistic case of interest. Also note that the strain gradient tensor cannot, e.g., capture the curvatures *u*_{x,xy}, *u*_{x,yx}, *u*_{y,xy}, and *u*_{y,yx}.

### Additional Numerical Results for 1D and 2D

Some numerical results discussed in previous sections are shown here. Figure 17 shows the profiles of residual stresses in the 1D bar if a constant sprain (localization-resisting) stiffness *κ* is used for postpeak damage. Figure 18 shows that both the triplet forces and the corresponding residual stresses have negligible values when the sprain (localization-resisting) stiffness *κ* values of the nodes near the broken elements are diminished at *t* = *t*_{max}. Figure 19 shows the contour of strain component $\epsilon xx$ calculated with original CBM using finite element meshes with various inclination angles but with the same element size in the vicinity of the potential damage area when the postpeak load drops to 0.8 of the maximum load value. Figure 19 illustrates the dependence of crack propagation path on the finite element mesh inclination when the original CBM is used. The dependence of crack propagation path on mesh inclination is not a physical phenomenon and thus requires the sCBM to address.

### Problems With Alternative Formulations

It may be noted that the extra sprain (localization-resisting) energy Φ, introduced in Eq. (3) for 1D or Eq. (17) for 2D or 3D, bears partial similarity with the extra energy of geometrically necessary dislocations [49] in flexure and torsion [53], which, aside from the strain-gradient tensor invariant, was also considered as a function of the invariant of the material rotation gradient tensor in 3D. In 2D, though, the material rotation vectors are coaxial and the variation of their magnitude within the plane causes merely a mismatch of displacement gradient *u*_{i,j} which, in turn, causes dependence on the second gradient tensor, *u*_{i,jk}. Therefore, there is no need to complicate the formulation by introducing the tensors of material rotation gradient.

*κ*= constant. But then the initial tangential stiffness,

*κ*

_{tan}, right after exceeding the threshold

*C*would vary with

*u*

_{i,j}rather than being a constant, and to counter this variability would complicate modeling. This is clear by noting that, in 1D,

### Simple Bilinear Elastic-Softening Constitutive Model

*Q*is the transformation matrix for rotation of coordinates. For the vector of principal strains $\epsilon p=[\epsilon x\u2032,\epsilon y\u2032,0]T$, we consider only the positive parts by defining $\epsilon p+(i)=max(\epsilon p(i),0)$,

*i*= 1, 2. Then the effective principal strain is given by $\epsilon eff=\epsilon p+T\u22c5\epsilon p+$. Further, we define $\epsilon pre$ as the maximum effective strain that has occurred previously up to the current time during the loading history. Also we define $\xi =max(\epsilon eff,\epsilon pre)$. Then the damage parameter

*ω*is given by

**is the elastic moduli matrix for 2D plane stress condition. Equation (D1) gets simplified for the one-dimensional case where only one strain component is considered.**

*D*