Smooth Crack Band Model—A Computational Paragon Based on Unorthodox Continuum Homogenization

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:

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

2. 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, $η$⁠, 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, $η$⁠, and sprain energy density, Φ. Tensor $η$ 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.

In standard continuum mechanics, the homogenized free energy density, $Ψ¯,$ 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 $Ψ¯$ is a function of, and only of, the strain tensor, $ε$⁠.

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′=du/dx$ from $u1′$ to $u2′$⁠. 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, $Δu′=u2′−u1′$ (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₀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₀u_x,yx of the y-derivative of the x-displacement over length l₀ 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 $ε′$⁠, it does not in two or three dimensions as discussed later.

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:

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

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

The energy-based homogenization is characterized by the third-order tensor, $η$⁠, 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.

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.

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.

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¯$ 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¯$⁠, with u″ given by Eq. (4), may be described as follows:

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.

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.

In a planar sheet of thickness t, with Cartesian coordinates x₁ = x, x₂ = y, the displacement vector $u$ has Cartesian components u₁ = u_x and u₂ = 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 $∇2ux=u1,jj$ and $∇2uy=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).

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:

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_xr,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.

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.

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, $η$⁠. Tensor α_ijη_ijk is a generalization of the vector of Laplacian $∇2uk$⁠.

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.

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.

Interestingly, Eq. (43) is formally identical to the bending equation for a beam of bending stiffness $κl02$⁠, under distributed transverse load q(x).² 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).

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 the original CBM, the single-element width of the crack band is taken equal to the material characteristic length, l₀, and the strain, $ε=u′$⁠, is uniform across the band. In the present sCBM, l₀ 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).

For tension $u′=ε>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.

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₀ is increased linearly with time t and reaches u_0max = 4 × 10⁻⁴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 $ε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 $εc=2Gf/(l0ft)$⁠. The material characteristic length l₀ in Eq. (3) is taken to be the aggregate size l₀ = 50 mm, and then, $εc/ε0=2.5$⁠. The problem is static, but the explicit dynamic relaxation algorithm with artificial mass density ρ = 2.6 × 10⁻⁹ ton/mm³, damping coefficient 10⁻³ ton/s, and time-step $Δt=0.8h/E/ρ$ 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 κ₀ = 10 MPa in Eq. (13), second gradient threshold C = 10⁻³, 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₀ = 1), Fig. 7 for an element size of 6 mm (h/l₀ = 3/25), and Fig. 8 for an element size of 2 mm (l₀/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₀/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₀/ = 1; 7 elements in Fig. 7(b) with h/l₀ = 3/25; and 25 elements in Fig. 8(b) with h/l₀ = 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₀ = 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).

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⁻³. 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₀/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.

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₀ in Eq. (17) is taken to be l₀ = 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⁻⁵. For the simulated fracture to be quasi-static, we ensure with mass density ρ′ = 2.6 × 10⁻⁵ ton/mm³ that the ratio of kinetic energy to strain energy is less than 10⁻³ 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 $ε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 $ε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 $ε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 $ε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₀ 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.

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

2. In continuum mechanics of damage of quasibrittle materials, the total energy density, $Ψ¯$⁠, 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₀. 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.

3. Below a certain threshold C, the sprain (localization-resisting) energy Φ vanishes and, above C, it is a function of the sprain tensor times l₀. During the postpeak softening, the corresponding localization-resisting (or sprain) stiffness, κ, needs to be reduced as the damage, characterized by volumetric strain, increases.

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

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

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

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.

There are no conflicts of interest.

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

Equations (A1) and (A2) use a spread-out change of volumetric strain, $3l0ε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.

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

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.