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 Φ(l0η), 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 Φ(l0η). 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 [13], 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 [59], and phase-field models [1014] 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) [1518], for which a realistic damage constitutive model, such as the microplane model M7 for concrete or shale [1922], 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.

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, Ψ¯, 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, ε.

Fig. 1
(a) If free energy density Φ depends on the change of length of infinitesimal line segment only, it must be a function of strain tensor only, (b) if free energy density Φ also depends on the change of displacement gradient u′, it must also depend on the second displacement gradient u″, (c) the change of displacement u″ must be distributed over a finite material length l0 characterizing heterogeneity (the same applies in 2D to duy/dx), and (d) a simple example of spatial periodic strain distribution in an elastic bar
Fig. 1
(a) If free energy density Φ depends on the change of length of infinitesimal line segment only, it must be a function of strain tensor only, (b) if free energy density Φ also depends on the change of displacement gradient u′, it must also depend on the second displacement gradient u″, (c) the change of displacement u″ must be distributed over a finite material length l0 characterizing heterogeneity (the same applies in 2D to duy/dx), and (d) a simple example of spatial periodic strain distribution in an elastic bar
Close modal

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=u2u1 (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′ = l0u″. 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 Δux,y = l0ux,yx of the y-derivative of the x-displacement over length l0 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.

Consequently, as our basic hypothesis for the modeling of distributed softening damage and fracture, we postulate that the total homogenized Helmholtz free energy density of the material, Ψ¯, must be a sum of two terms:
Ψ¯(ε,l0η)=Ψ(ε)+Φ(l0η)
(1)
Here, ε is the strain tensor, and Φ(l0η) is the 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 η is a third-order tensor of Cartesian components ηijk = uk,ij. It represents the third-order tensor of the second gradient of displacement vector uk (subscripts k = 1, 2, 3 label the Cartesian coordinates xk, and the summation rule is implied). For the sake of brevity, η 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 Ψ(ε), Φ(l0η) 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 η behaves differently than the strain-gradient tensor μ, with components μkij=εij,k (where, for small strains, ε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:

  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.

To clarify the difference in a primitive way, consider a spatial strain distribution ε=εstiffness+asin(2πx/l) in an elastic bar (εstiffness,a,l being constant, n = 1, 2, 3, …), see Fig. 1(d). The average of strain energy density based on stiffness homogenization is Ψ(ε¯)=(E/2)εstiffness2 (solid horizontal line in Fig. 1(d). But the averaging of energy density gives a higher value, Ψ¯(ε(x))=(1/nl)0nl(E/2)(εstiffness+asin(2πx/l))2dx=(E/2)ε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:
Φ(lu)=E4a2
(2)

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,4248], 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.

First, we discuss a one-dimensional (1D) case since it is instructive by its simplicity. Let us consider the 1D case of strain localization in a heterogeneous bar in direction x1x. 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 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:
Φ(x)=κ2l0|u(x)|C2
(3)
where C is the “sprain” threshold, and κ is the localization-resisting stiffness (or “sprain” stiffness). Its dimension is that of stress (i.e., Pa); l0 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.

Let us introduce a uniform subdivision of coordinate x1x into line elements of size Δx = h delimited by nodes r = 1, 2, …, N. In terms of nodal displacement ur, the second displacement gradient, i.e., the second derivative at node r, [d2u(x)/dx2]r=u(xr)=ur, is approximated by
urur+12ur+ur1h2
(4)

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.

The discrete approximation of the localization-resisting sprain energy density at node r is:
Φ(xr)=κ2l0|ur12ur+ur+1|h2C2
(5)
Differentiation with respect to the nodal displacements yields the nodal sprain forces that need to be applied to resist the excess energy of damage localization, i.e.:
fr=bthΦ(x)ur,fr1=bthΦ(x)ur1,fr+1=bthΦ(x)ur+1
(6)

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

By carrying out the differentiation in Eq. (6) with Eq. (5), one gets the scaling rule of the nodal forces:
fr=btκl0hf0¯,fr1=fr+1=fr2
(7)
where f0¯ is a dimensionless constant:
f0¯=2l0|ur12ur+ur+1|h2Csgn(u)
(8)

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 fr 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:

foru0:f0¯=max[2(l0uC),0]
(9)
foru<0:f0¯=min[2(l0u+C),0]
(10)
Figure 2 depicts graphically the triplet of axial nodal forces
(fr1,fr,fr+1)=(fr2,fr,fr2)forallr=2,3,4,,N1
(11)
applied at each triplet of nodes. Note that these forces are always self-equilibrated, i.e.,
fr1+fr+fr+1=0
(12)
(and so, in 2D or 3D generalization, according to St. Venant principle, their effect decays exponentially with distance). The adjacent triplets are overlapping, and the forces from the adjacent triplets must be superposed and assembled in a way similar to finite element stiffness matrices. Within intervals of a uniform second gradient u″, the nodal forces sum up to zero. Further, note that lim h→0fr = ∞ in 1D (this changes for 2D if the element sizes in both directions tend to zero, as demonstrated in Sec. 2.3).
Fig. 2
Triplet forces countering the excessive positive second displacement gradient u″ at node r in 1D
Fig. 2
Triplet forces countering the excessive positive second displacement gradient u″ at node r in 1D
Close modal

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

Keeping κ constant is sufficient to handle early postpeak damage, but not full softening. Under excessive normal strain ε=u, the material degrades and κ decreases. So we express κ, the sprain (localization-resisting) stiffness, as follows:
κ(λ)=λκ0
(13)
where κ0 is a material constant and λ is a stiffness parameter (representing an internal variable). The simplest realistic evolution law for λ is expressed as follows:
Δλ=2kd2λεΔε(ε=u)
(14)
where kd is a dimensionless material constant. The 〈..〉 ensures that the damage cannot decrease. For monotonic ε(t) and infinitesimal Δεdε, this integrates to yield the bell-shaped curve:
λ=e(u/kd)2
(15)
which is proportional to the Gaussian probability density function. Here, by having u2 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 kd ≈ 3 ft/E, which gives a rough estimation of the order of magnitude.
The discrete form of Eq. (14) for loading increment Δt is given as follows:
Δλr=2kd2λrur+1ur12hΔur+1Δur12h
(16)
For the boundary nodes, the calculation of λr must be omitted.

In step-by-step integration, a simultaneous calculation of new nodal values of ur 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 x1 = x, x2 = y, the displacement vector u has Cartesian components u1 = ux and u2 = uy. 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 ux(x, y) and uy(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 ux(x, y) and uy (x, y), are vectors (while the mean curvature of each surface is an invariant, equal to one half of the Laplacian).

To generalize Eq. (3) to 2D (or 3D), one must further realize that, for tensorial consistency, threshold C must now be replaced by a vector, Ci, even though, in the case of isotropy, its components must be equal, i.e., C1 = C2 = C = constants = dimensionless thresholds. So, the proper generalization of Eq. (3) giving the sprain (or localization-resisting) energy density as a tensorial invariant is
Φ(x)=κ2l0ui,jjCil0ui,kkCi
(17)
where 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:
Φ(x,y)=κ2[l0(ux,xx+ux,yy)Cx2+l0(uy,xx+uy,yy)Cy2]
(18)
=κ2(l02uxCx2+l02uyCy2)=κ2|ξ|2
(19)
where |ξ|2=ξiξi,ξi=<l02uiCi> and Cx=Cy=C1=C2 in case of isotropy.The metric dimensions of Φ and κ are N/m2 (or J/m3). 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., ux,yy and uy,xx) nor the Laplacian could be expressed in terms of the strain-gradient tensor.
It may be noted that the sprain tensorη, which represents the third-order second-gradient tensor of the displacement vector and is defined as
ηijk=uk,ij
(20)
has five invariants. They have similar forms as the invariants of the strain-gradient tensor μ with components μkij=ε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 hx and hy 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:

ux,xx=(uxr1,s2uxr,s+uxr+1,s)hx2,uy,xx=(uyr1,s2uyr,s+uyr+1,s)hx2
(21)
ux,yy=(uxr,s12uxr,s+uxr,s+1)hy2,uy,yy=(uyr,s12uyr,s+uyr,s+1)hy2
(22)
where coordinate subscripts preceded by a comma denote derivatives, but not for the nodal indices 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:
fxr,s=hxhytΦr,suxr,s,fxr1,s=hxhytΦr,suxr1,s,fxr+1,s=hxhytΦr,suxr+1,s
(23)
fyr,s=hxhytΦr,suyr,s,fyr1,s=hxhytΦr,suyr1,s,fyr+1,s=hxhytΦr,suyr+1,s
(24)
fxr,s1=hxhytΦr,suxr,s1,fxr,s+1=hxhytΦr,suxr,s+1
(25)
fyr,s1=hxhytΦr,suyr,s1,fyr,s+1=hxhytΦr,suyr,s+1
(26)
For the individual nodes, calculation of the derivatives of Laplacian of ux yields,
(2ux)(uxr,s)=2hx22hy2,(2ux)(uxr1,s)=(2ux)(uxr+1,s)=1hx2,(2ux)(uxr,s1)=(2ux)(uxr,s+1)=1hy2
(27)
For uy, the calculation is similar. Introducing dimensionless Laplacians,
X=l0|2ux|,Y=l0|2uy|
(28)
we obtain, after rearrangements, the following localization-resisting in-plane nodal forces in the set of five nodes centered at node (r, s):
fxr,s=2κ(λ)l0t(hyhx+hxhy)XC=fx
(29)
fyr,s=2κ(λ)l0t(hyhx+hxhy)YC=fy
(30)
fxr1,s=fxr+1,s=κ(λ)l0t(hyhx)XC
(31)
fxr,s1=fxr,s+1=κ(λ)l0t(hxhy)XC
(32)
fyr1,s=fyr+1,s=κ(λ)l0t(hyhx)YC
(33)
fyr,s1=fyr,s+1=κ(λ)l0t(hxhy)YC
(34)
Physically, these forces represent the resultants of localization-resisting stresses in directions x and y (see Fig. 3).
Fig. 3
Forces generated by twin triplets with the curvature value at center node (r, s) in a 2D domain with local coordinate system (x′, y′) and global coordinate system (x, y)
Fig. 3
Forces generated by twin triplets with the curvature value at center node (r, s) in a 2D domain with local coordinate system (x′, y′) and global coordinate system (x, y)
Close modal

The scaling rule of nodal forces in 2D is given by Eqs. (29)(34). When both steps hx and hy are scaled in proportion, the nodal forces fx and fy do not change with the step size. When hx is refined while hy does not change, the fxr,s scales, for small hx/hy, asymptotically in proportion to 1/hx, the same as in the one-dimensional case, and asymptotically for large hx/hy in proportion to hx. Similar scaling properties hold for other components.

It remains to generalize Eq. (14) to two dimensions. The sprain (localization-resisting) stiffness parameter must be an invariant of strain εi,j=(ui,j+uj,i)/2. The logical and the simplest is to take the (linearized) area strain εA=(ε11+ε22)/2. The degradation law is again κ(λ) = λκ0 and the 2D evolution law for λ is
Δλ=2kd2λεAΔεAwithεA=12(u1,1+u2,2)
(35)
For monotonic εA(t) and infinitesimal ΔεAdεA, this integrates as follows:
λ=e(εA/kd)2
(36)
where kd is an empirical dimensionless material constant. The discrete form of Eq. (35) in increment Δt is expressed as follows:
Δλr,s=2kd2(λr,s)12(ur1,sur+1,s2hx+ur,s1ur,s+12hy)(Δur1,sΔur+1,s2hx+Δur,s1Δur,s+12hy)
(37)

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

Material orthotropy, a special case of anisotropy, may be characterized by tensor αij. Equation (17) may then be simply generalized as follows:
Φ(x,y)=κ2l0αijuk,ijCkl0αmnuk,mnCk
(38)
where Ck are the thresholds for the orthotropy directions, and the sprain tensor ηijk = uk,ij is defined by Eq. (20) (we skip discussing the conversion of κ to general anisotropy).

Material isotropy is the special case for Ci = 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 Ci) 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.

The differential equation underlying the sCBM can be derived variationally using the Gauss integral theorem. For the sake of brevity, a one-dimensional bar is considered here. We want to minimize the potential energy Π:
Π=bt0L[e=0e=uσ(e)de+12κ(λ)l0uC2pu]dx
(39)
where L,b,andt are bar length, width, and thickness, respectively, −pu is the potential energy of the axial distributed load p applied per unit volume of material, u=du/dx, and u=d2u/dx2. Let us restrict ourselves to the case that u″ <−C/l0 (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:
δΠ=bt0L[σδuκ(λ)l0uCl0δupδu]dx=0
(40)
=bt0L[σδuκ(λ)l02uδupδu]dx+[]0L=0
(41)
=bt0L(σ+κ(λ)l02uIVp)δudx+[]0L=0
(42)
The boundary terms []0L stemming from the integration by parts are omitted since they give the boundary conditions, while we seek here only the differential equation. Since δΠ must vanish for any δu, the expression in parenthesis must vanish, too. Thus, we obtain the differential equation:
κl02d4u(x)dx4=q(x)whereq(x)=dσ(x)dx+p(x)
(43)

Interestingly, Eq. (43) is formally identical to the bending equation for a beam of bending stiffness κ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).

To obtain the resisting nodal forces fr, Eq. (43) may be rewritten as uIV=(u)=(σ+p)/(κl02)=q/(κl02). Then, replacing (·)″ by its finite difference approximation at node r, we have
ur12ur+ur+1h2=qrκ(λ)l02
(44)
where qr is the q value at node r. Now, to relate ur to the nodal force, we take a derivative of Eq. (3) and obtain (for ur < 0, hill):
fr=2btκ(λ)l0hl0urC
(45)
Then, expressing ur and similarly ur−1, ur+1, from this equation, we obtain, upon substitution into Eq. (44) and after rearrangements:
fr12fr+fr+1=2hbtqr
(46)

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:

Π=tΩ(εσij(ε)dεij+κ2l0ui,jjCil0ui,kkCipiui)dxdy
(47)
Calculating the variation δΠ 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, l0, and the strain, ε=u, is uniform across the band. In the present sCBM, l0 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).

Fig. 4
(a) A parabolic displacement profile with a coarse mesh, (b) a parabolic displacement profile with a fine mesh, (c) first derivative of displacement in (a), and (d) first derivative of displacement in (b) showing hill-type localization
Fig. 4
(a) A parabolic displacement profile with a coarse mesh, (b) a parabolic displacement profile with a fine mesh, (c) first derivative of displacement in (a), and (d) first derivative of displacement in (b) showing hill-type localization
Close modal
For the present sCBM, we need first to estimate a suitable value of threshold C (to be later recalibrated by experiments). We introduce a subdivision of l0 into many elements. In that case, the distribution of strain ε=u 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 umax = umax/(l0/2) = C. From this we get the estimate (a crude one):
C2ftE
(48)

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 fc is about ten times higher than the tensile strength ft. 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 u0 is increased linearly with time t and reaches u0max = 4 × 10−4L at t/tmax = 1.0 (Fig. 5(b) where the total time tmax = 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 ft, and the material in the rest of the bar is assumed to be elastic up to a tensile strength ft = 4 MPa. The corresponding strain at the tensile strength limit is ε0=ft/E.

Fig. 5
One-dimensional bar under tension. In one portion (1/20 of the total length, indicated by black color), the material strength is 0.99 of that in the rest of the bar: (a) symmetric constant-rate growth is imposed at both ends, (b) the applied displacement u0 varies linearly with time t, and (c) bilinear stress versus strain response of the material of the bar in (a), ending with σ/ft = 2.5 × 10−16 after complete softening
Fig. 5
One-dimensional bar under tension. In one portion (1/20 of the total length, indicated by black color), the material strength is 0.99 of that in the rest of the bar: (a) symmetric constant-rate growth is imposed at both ends, (b) the applied displacement u0 varies linearly with time t, and (c) bilinear stress versus strain response of the material of the bar in (a), ending with σ/ft = 2.5 × 10−16 after complete softening
Close modal

Simple linear softening, corresponding to fracture energy Gf = 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 l0 in Eq. (3) is taken to be the aggregate size l0 = 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−9 ton/mm3, damping coefficient 10−3 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 κ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 = tmax are shown in Fig. 6 for element sizes of 50 mm (h/l0 = 1), Fig. 7 for an element size of 6 mm (h/l0 = 3/25), and Fig. 8 for an element size of 2 mm (l0/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.

Fig. 6
Results for one-dimensional bar meshed with elements of size 50 mm (h/l0 = 1) at t = tmax in Fig. 5: (a) displacement profile, (b) strain profile, (c) profile of the second gradient of displacement, and (d) profile of superposed triplet forces
Fig. 6
Results for one-dimensional bar meshed with elements of size 50 mm (h/l0 = 1) at t = tmax in Fig. 5: (a) displacement profile, (b) strain profile, (c) profile of the second gradient of displacement, and (d) profile of superposed triplet forces
Close modal
Fig. 7
Results for a one-dimensional bar meshed with elements of size 6 mm (h/l0=3/25) at t = tmax in Fig. 5: (a) displacement profile, (b) strain profile, (c) profile of the second gradient of displacement, and (d) profile of superposed triplet forces
Fig. 7
Results for a one-dimensional bar meshed with elements of size 6 mm (h/l0=3/25) at t = tmax in Fig. 5: (a) displacement profile, (b) strain profile, (c) profile of the second gradient of displacement, and (d) profile of superposed triplet forces
Close modal
Fig. 8
Results for the one-dimensional bar meshed with elements of size 2 mm (h/l0=1/25) at t = tmax in Fig. 5: (a) displacement profile, (b) strain profile, (c) profile of the second gradient of displacement, and (d) profile of superposed triplet forces
Fig. 8
Results for the one-dimensional bar meshed with elements of size 2 mm (h/l0=1/25) at t = tmax in Fig. 5: (a) displacement profile, (b) strain profile, (c) profile of the second gradient of displacement, and (d) profile of superposed triplet forces
Close modal

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 68 show that larger l0/h induces damage distribution over more elements representing the softening “damage zone” (or sprain zone); see the subfigures (a), (b), and (c) in Figs. 68. Across the “sprain zones,” there is 1 element in Fig. 6(b) with h/l0/ = 1; 7 elements in Fig. 7(b) with h/l0 = 3/25; and 25 elements in Fig. 8(b) with h/l0 = 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 l0 = 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. 68). 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 hyhx).

Fig. 9
Zoomed-in (a) Fig. 6(b), (b) Fig. 7(b), and (c) Fig. 8(b); a smoothed hill-type localization of damage is demonstrated in (c)
Fig. 9
Zoomed-in (a) Fig. 6(b), (b) Fig. 7(b), and (c) Fig. 8(b); a smoothed hill-type localization of damage is demonstrated in (c)
Close modal

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 kd = 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 (l0/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).

Fig. 10
(a) Diagram of reaction force versus displacement at the right end of the bar with the boundary condition specified in Fig. 5 and (b) a zoomed-in residual strain profile near the center of the bar in Fig. 11(e)
Fig. 10
(a) Diagram of reaction force versus displacement at the right end of the bar with the boundary condition specified in Fig. 5 and (b) a zoomed-in residual strain profile near the center of the bar in Fig. 11(e)
Close modal
Fig. 11
Profiles at maximum reaction force (first row) denoted by the dot in Fig. 10(a) and profiles at t = tmax after rupture (second row). The element size is 2 mm (h/l0=1/25): (a) displacement profile at maximum load, (b) strain profile ε at maximum load, (c) normalized sprain (localization-resisting) stiffness profile at maximum load, (d) displacement profile at t = tmax, (e) strain profile ε at t = tmax, and (f) normalized sprain (localization-resisting) stiffness profile at t = tmax.
Fig. 11
Profiles at maximum reaction force (first row) denoted by the dot in Fig. 10(a) and profiles at t = tmax after rupture (second row). The element size is 2 mm (h/l0=1/25): (a) displacement profile at maximum load, (b) strain profile ε at maximum load, (c) normalized sprain (localization-resisting) stiffness profile at maximum load, (d) displacement profile at t = tmax, (e) strain profile ε at t = tmax, and (f) normalized sprain (localization-resisting) stiffness profile at t = tmax.
Close modal

The second row of Fig. 11 shows the corresponding profiles at t = tmax 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/tmax = 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.

An effective way to measure the fracture properties of quasibrittle materials, which may be easily generalized to measure the crack-parallel stress effect [1,2], is the three-point bend test. The specimen is a notched beam of depth D, width 2L + 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 a0 the initial notch length. For the explicit integration, the rate boundary conditions are as follows:
u˙y=u˙atx=L+W/2,y=D
(49a)
u˙x=0,u˙y=0atx=0,y=0
(49b)
u˙y=0atx=2L+W,y=0
(49c)
On the rest of the surface, they are traction free.
Fig. 12
Sketch of three-point bend specimen (with unit thickness). End supports not shown.
Fig. 12
Sketch of three-point bend specimen (with unit thickness). End supports not shown.
Close modal

The material parameters are Young’s modulus E = 30 GPa, fracture energy Gf = 0.05 N/mm, and tensile strength ft = 4 MPa. The characteristic length l0 in Eq. (17) is taken to be l0 = 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/mm3 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 ε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.

Fig. 13
Contours of strain component εxx normal to notch obtained with original CBM for specimens with various mesh inclination angles: (a) α = 0 deg, (b) α = 15 deg, (c) α = 30 deg, and (d) α = 45 deg. (e), (f), (g), and (h) are zoomed-in plots of (a), (b), (c), and (d), respectively
Fig. 13
Contours of strain component εxx normal to notch obtained with original CBM for specimens with various mesh inclination angles: (a) α = 0 deg, (b) α = 15 deg, (c) α = 30 deg, and (d) α = 45 deg. (e), (f), (g), and (h) are zoomed-in plots of (a), (b), (c), and (d), respectively
Close modal

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

Fig. 14
Contours of strain component εxx normal to notch obtained with sCBM for specimens with various mesh inclination angles: (a) α = 0 deg, (b) α = 15 deg, (c) α = 30 deg, and (d) α = 45 deg. (e), (f), (g), and (h) are zoomed-in plots of (a), (b), (c), and (d), respectively. Note negligible mesh inclination effect.
Fig. 14
Contours of strain component εxx normal to notch obtained with sCBM for specimens with various mesh inclination angles: (a) α = 0 deg, (b) α = 15 deg, (c) α = 30 deg, and (d) α = 45 deg. (e), (f), (g), and (h) are zoomed-in plots of (a), (b), (c), and (d), respectively. Note negligible mesh inclination effect.
Close modal

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

Fig. 15
Contours of stress component σxx normal to notch obtained with sCBM for specimens with various mesh inclination angles: (a) α = 0 deg, (b) α = 15 deg, (c) α = 30 deg, and (d) α = 45 deg. (e), (f), (g), and (h) are zoomed-in plots of (a), (b), (c), and (d), respectively. Note negligible mesh inclination effect.
Fig. 15
Contours of stress component σxx normal to notch obtained with sCBM for specimens with various mesh inclination angles: (a) α = 0 deg, (b) α = 15 deg, (c) α = 30 deg, and (d) α = 45 deg. (e), (f), (g), and (h) are zoomed-in plots of (a), (b), (c), and (d), respectively. Note negligible mesh inclination effect.
Close modal

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.

Fig. 16
Normalized load P/(ftD2) versus displacement u/D of the loading point at the top of the beam in Fig. 12: (a) obtained with original CBM and (b) obtained with sCBM
Fig. 16
Normalized load P/(ftD2) versus displacement u/D of the loading point at the top of the beam in Fig. 12: (a) obtained with original CBM and (b) obtained with sCBM
Close modal

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

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

2

For a solution by finite difference method, one would need to introduce the approximation uIV ≈ (ur−2 − 4ur−1 + 6ur − 4ur+1 + ur+2)/h4, 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

The physical motivation for Φ, used to setup Eq. (17), is that the 3D generalization of curvature is a Laplacian. But it is helpful to clarify why other invariants of tensor η are not appropriate. This tensor has five independent invariants listed in Eq. 2.5 of [46] (and in Eq. 1 of [51]). Equation (17) for the sprain energy corresponds, in its form, to the third of these five, which is ηiikηjjk (also [51]). To explain why the other four independent invariants of η are inapplicable, consider the energy densities Φa, …, Φd corresponding to the other four invariants of tensor η:
ηijjηikkΦa3l0εV,iCi3l0εV,iCi
(A1)
ηiikηkjjΦbl02ukCk3l0εV,kCk
(A2)
ηijkηijkΦcl0uk,ijCijkl0uk,ijCijk
(A3)
ηijkηkjiΦdl0uk,ijCkijl0ui,kjCikj
(A4)
where εV,i=(uj,j),i/3=uj,ji/3=uj,ij/3, which denotes the gradient of the volumetric (or hydrostatic) strain and is different from the Laplacian 2ui=ui,jj.

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 Cijk is dubious (on the other hand, a vector threshold Ci in Eq. (17) is physically necessary in the case of orthotropy). The components such as (u1,2),3 (u2,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 ux,xy, ux,yx, uy,xy, and uy,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 = tmax. 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.

Fig. 17
Profiles of residual stresses in the 1D bar with elements of size (a) 50 mm (h/l0=1), (b) 6 mm (h/l0=3/25), and (c) 2 mm (h/l0=1/25) if a constant sprain (localization-resisting) stiffness κ is used for postpeak damage
Fig. 17
Profiles of residual stresses in the 1D bar with elements of size (a) 50 mm (h/l0=1), (b) 6 mm (h/l0=3/25), and (c) 2 mm (h/l0=1/25) if a constant sprain (localization-resisting) stiffness κ is used for postpeak damage
Close modal
Fig. 18
Profiles at t = tmax after rupture. Note that the residual stress in (c) is negligibly small. The element size is 2 mm (h/l0=1/25): (a) profile of superposed triplet forces (zoomed-in plot), (b) second derivative of displacement, and (c) stress.
Fig. 18
Profiles at t = tmax after rupture. Note that the residual stress in (c) is negligibly small. The element size is 2 mm (h/l0=1/25): (a) profile of superposed triplet forces (zoomed-in plot), (b) second derivative of displacement, and (c) stress.
Close modal
Fig. 19
Contours of stress component σxx normal to notch obtained with original CBM for specimens with various mesh inclination angles: (a) α = 0 deg, (b) α = 15 deg, (c) α = 30 deg, and (d) α = 45 deg. (e), (f), (g), and (h) are zoomed-in pictures of (a), (b), (c), and (d), respectively
Fig. 19
Contours of stress component σxx normal to notch obtained with original CBM for specimens with various mesh inclination angles: (a) α = 0 deg, (b) α = 15 deg, (c) α = 30 deg, and (d) α = 45 deg. (e), (f), (g), and (h) are zoomed-in pictures of (a), (b), (c), and (d), respectively
Close modal

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 ui,j which, in turn, causes dependence on the second gradient tensor, ui,jk. Therefore, there is no need to complicate the formulation by introducing the tensors of material rotation gradient.

One might also think of using
Φ=κ2(l02ui,jjui,kkC2)
(C1)
where κ = constant. But then the initial tangential stiffness, κtan, right after exceeding the threshold C would vary with ui,j rather than being a constant, and to counter this variability would complicate modeling. This is clear by noting that, in 1D,
κtan=u[κ2(l02u2C2)]=κl02uconstant
(C2)

Simple Bilinear Elastic-Softening Constitutive Model

For the two-dimensional case, we consider a strain vector ε=[εxx,εyy,γxy]T. The principal direction is given by θp=(1/2)arctan(γxy/(εxxεyy)), obtained by setting the principal shear strain to 0, and the principal strain vector is given by εp=Qε, where Q is the transformation matrix for rotation of coordinates. For the vector of principal strains εp=[εx,εy,0]T, we consider only the positive parts by defining εp+(i)=max(εp(i),0), i = 1, 2. Then the effective principal strain is given by εeff=εp+Tεp+. Further, we define εpre as the maximum effective strain that has occurred previously up to the current time during the loading history. Also we define ξ=max(εeff,εpre). Then the damage parameter ω is given by
ω={1ξ<ε0[εc(ξε0)][(εcε0)ξ]ε0ξ<εthr[εc(εthrε0)][(εcε0)εthr]ξεthr
(D1)
The stress vector is calculated as σ=(1ω)Dε, where D 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.

References

1.
Nguyen
,
H. T.
,
Pathirage
,
M.
,
Rezaei
,
M.
,
Issa
,
M.
,
Cusatis
,
G.
, and
Bažant
,
Z. P.
,
2020
, “
New Perspective of Fracture Mechanics Inspired by Gap Test With Crack-Parallel Compression
,”
Proc. Natl. Acad. Sci. USA
,
117
(
25
), pp.
14015
14020
.
2.
Nguyen
,
H. T.
,
Pathirage
,
M.
,
Cusatis
,
G.
, and
Bažant
,
Z. P.
,
2020
, “
Gap Test of Crack-Parallel Stress Effect on Quasibrittle Fracture and Its Consequences
,”
ASME J. Appl. Mech.
,
87
(
7
), p.
071012
.
3.
Nguyen
,
H. T.
,
Dönmez
,
A. A.
, and
Bažant
,
Z. P.
,
2021
, “
Structural Strength Scaling Law for Fracture of Plastic-Hardening Metals and Testing of Fracture Properties
,”
Extreme Mech. Lett.
,
43
, p.
101141
.
4.
Mazars
,
J.
,
1981
, “
Mechanical Damage and Fracture of Concrete Structures
,”
Adv. Fracture Res.
,
4
, pp.
1499
1506
.
5.
Barenblatt
,
G. I.
,
1959
, “
The Formation of Equilibrium Cracks During Brittle Fracture. General Ideas and Hypotheses. Axially-Symmetric Cracks
,”
J. Appl. Math. Mech.
,
23
(
3
), pp.
622
636
.
6.
Barenblatt
,
G. I.
,
1962
, “
The Mathematical Theory of Equilibrium Cracks in Brittle Fracture
,”
Adv. Appl. Mech.
,
7
, pp.
55
129
.
7.
Planas
,
J.
, and
Elices
,
M.
,
1992
, “
Asymptotic Analysis of a Cohesive Crack: 1. Theoretical Background
,”
Int. J. Fracture
,
55
(
2
), pp.
153
177
.
8.
Planas
,
J.
, and
Elices
,
M.
,
1993
, “
Asymptotic Analysis of a Cohesive Crack: 2. Influence of the Softening Curve
,”
Int. J. Fracture
,
64
(
3
), pp.
221
237
.
9.
Bažant
,
Z. P.
, and
Planas
,
J.
,
1998
,
Fracture and Size Effect in Concrete and Other Quasibrittle Materials
,
CRC Press
,
Boca Raton, FL
.
10.
Francfort
,
G. A.
, and
Marigo
,
J. -J.
,
1998
, “
Revisiting Brittle Fracture as an Energy Minimization Problem
,”
J. Mech. Phys. Solids.
,
46
(
8
), pp.
1319
1342
.
11.
Bourdin
,
B.
,
Francfort
,
G. A.
, and
Marigo
,
J.-J.
,
2000
, “
Numerical Experiments in Revisited Brittle Fracture
,”
J. Mech. Phys. Solids.
,
48
(
4
), pp.
797
826
.
12.
Bourdin
,
B.
,
Francfort
,
G. A.
, and
Marigo
,
J.-J.
,
2008
, “
The Variational Approach to Fracture
,”
J. Elast.
,
91
(
1–3
), pp.
5
148
.
13.
Wu
,
J.-Y.
,
2017
, “
A Unified Phase-Field Theory for the Mechanics of Damage and Quasi-Brittle Failure
,”
J. Mech. Phys. Solids.
,
103
, pp.
72
99
.
14.
Wu
,
J.-Y.
,
Huang
,
Y.
,
Zhou
,
H.
, and
Nguyen
,
V. P.
,
2021
, “
Three-Dimensional Phase-Field Modeling of Mode I+ II/III Failure in Solids
,”
Comput. Methods. Appl. Mech. Eng.
,
373
, p.
113537
.
15.
Bažant
,
Z. P.
,
1982
, “
Crack Band Model for Fracture of Geomaterials
,”
4th International Conference on Numerical Methods in Geomech
,
University of Alberta, Edmonton
.
16.
Bažant
,
Z. P.
, and
Oh
,
B. H.
,
1983
, “
Crack Band Theory for Fracture of Concrete
,”
Matériaux et construction
,
16
(
3
), pp.
155
177
.
17.
Červenka
,
J.
,
Bažant
,
Z. P.
, and
Wierer
,
M.
,
2005
, “
Equivalent Localization Element for Crack Band Approach to Mesh-Sensitivity in Microplane Model
,”
Int. J. Numer. Methods Eng.
,
62
(
5
), pp.
700
726
.
18.
Le
,
J.-L.
, and
Eliáš
,
J.
,
2016
, “
A Probabilistic Crack Band Model for Quasibrittle Fracture
,”
ASME J. Appl. Mech.
,
83
(
5
), p.
051005
.
19.
Caner
,
F. C.
,
Bažant
,
Z. P.
, and
ASCE
,
H. M.
,
2013
, “
Microplane Model M7 for Plain Concrete. I: Formulation
,”
Mechanics
,
139
(
12
), pp.
1714
1723
.
20.
Caner
,
F. C.
, and
Bažant
,
Z. P.
,
2013
, “
Microplane Model M7 for Plain Concrete. II: Calibration and Verification
,”
J. Eng. Mech.
,
139
(
12
), pp.
1724
1735
.
21.
Caner
,
F. C.
,
Bažant
,
Z. P.
, and
Wendner
,
R.
,
2013
, “
Microplane Model M7F for Fiber Reinforced Concrete
,”
Eng. Fract. Mech.
,
105
, pp.
41
57
.
22.
Bažant
,
Z. P.
,
Nguyen
,
H. T.
, and
Abdullah Dönmez
,
A.
,
2022
, “
Critical Comparison of Phase-Field, Peridynamics, and Crack Band Model M7 in Light of Gap Test and Classical Fracture Tests
,”
ASME J. Appl. Mech.
,
89
(
6
), p.
061008
.
23.
M7 subroutine
,
2013
, http://www.civil.northwestern.edu/people/bazant, Accessed December 2, 2022.
24.
Bažant
,
Z. P.
, and
Nguyen
,
H. T.
,
2022
, “
Proposal of a Model Index, MI, for Experimental Comparison of Fracture and Damage Models
,”
J. Eng. Mech.
,
149
(
11
).
25.
Bažant
,
Z. P.
, and
Jirásek
,
M.
,
2002
, “
Nonlocal Integral Formulations of Plasticity and Damage: Survey of Progress
,”
J. Eng. Mech.
,
128
(
11
), pp.
1119
1149
.
26.
Cusatis
,
G.
,
Bažant
,
Z. P.
, and
Cedolin
,
L.
,
2003
, “
Confinement-Shear Lattice Model for Concrete Damage in Tension and Compression: I. Theory
,”
J. Eng. Mech.
,
129
(
12
), pp.
1439
1448
.
27.
Cusatis
,
G.
,
Pelessone
,
D.
, and
Mencarelli
,
A.
,
2011
, “
Lattice Discrete Particle Model (LDPM) for Failure Behavior of Concrete. I: Theory
,”
Cement Concrete Compos.
,
33
(
9
), pp.
881
890
.
28.
Voigt
,
W.
,
1889
, “
Ueber Die Beziehung Zwischen den Beiden Elasticitätsconstanten Isotroper Körper
,”
Annalen der Physik
,
274
(
12
), pp.
573
587
.
29.
Reuss
,
A.
,
1929
, “
Berechnung der Fließgrenze Von Mischkristallen Auf Grund der Plastizitätsbedingung für Einkristalle
,”
ZAMM-J. Appl. Math. Mech./Zeitschrift für Angewandte Mathematik und Mechanik
,
9
(
1
), pp.
49
58
.
30.
Hashin
,
Z.
, and
Shtrikman
,
S.
,
1963
, “
A Variational Approach to the Theory of the Elastic Behaviour of Multiphase Materials
,”
J. Mech. Phys. Solids.
,
11
(
2
), pp.
127
140
.
31.
Eshelby
,
J. D.
,
1957
, “
The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems
,”
Proc. R. Soc. London., A.
,
241
(
1226
), pp.
376
396
.
32.
Hashin
,
Z.
,
1988
, “
The Differential Scheme and Its Application to Cracked Materials
,”
J. Mech. Phys. Solids.
,
36
(
6
), pp.
719
734
.
33.
Hill
,
R.
,
1965
, “
A Self-Consistent Mechanics of Composite Materials
,”
J. Mech. Phys. Solids.
,
13
(
4
), pp.
213
222
.
34.
Mori
,
T.
, and
Tanaka
,
K.
,
1973
, “
Average Stress in Matrix and Average Elastic Energy of Materials With Misfitting Inclusions
,”
Acta. Metall.
,
21
(
5
), pp.
571
574
.
35.
Eringen
,
A. C.
,
1966
, “
Linear Theory of Micropolar Elasticity
,”
J. Math. Mech.
,
15
(
6
), pp.
909
923
.
36.
Eringen
,
A. C.
,
1999
,
Microcontinuum Field Theories
,
Springer
,
New York
, pp.
101
248
.
37.
Bažant
,
Z.
, and
Christensen
,
M.
,
1972
, “
Analogy Between Micropolar Continuum and Grid Frameworks Under Initial Stress
,”
Int. J. Solids. Struct.
,
8
(
3
), pp.
327
346
.
38.
Bažant
,
Z. P.
, and
Cedolin
,
L.
,
1991
,
Stability of Structures: Elastic, Inelastic, Fracture, and Damage Theories
,
Oxford University Press
,
New York
.
39.
Christensen
,
R.
,
1969
, “
Viscoelastic Properties of Heterogeneous Media
,”
J. Mech. Phys. Solids.
,
17
(
1
), pp.
23
41
.
40.
Christensen
,
R. M.
,
1991
,
Mechanics of Composite Materials
,
Krieger Pub. Co
,
Malabar
.
41.
Dvorak
,
G.
,
2012
,
Micromechanics of Composite Materials
, Vol.
186
,
Springer Science & Business Media
,
Dordrecht
.
42.
Cosserat
,
E.
, and
Cosserat
,
F.
,
1909
,
Théorie des corps déformables
,
Librairie Scientifique A. Hermann et Fils
,
Paris
.
43.
Mindlin
,
R.D.
, and
Tiersten
,
H.F.
,
1962
, “
Effects of Couple-Stresses in Linear Elasticity
,”
Arch. Ration. Mech. Anal.
,
11
, pp.
415
448
.
44.
Toupin
,
R.
,
1962
, “
Elastic Materials With Couple-Stresses
,”
Arch. Rat. Mech. Anal.
,
11
(
1
), pp.
385
414
.
45.
Fleck
,
N.
, and
Hutchinson
,
J.
,
1993
, “
A Phenomenological Theory for Strain Gradient Effects in Plasticity
,”
J. Mech. Phys. Solids.
,
41
(
12
), pp.
1825
1857
.
46.
Fleck
,
N.
, and
Hutchinson
,
J.
,
1997
, “
Strain Gradient Plasticity
,”
Adv. Appl. Mech.
,
33
, pp.
295
361
.
47.
Gao
,
H.
,
Huang
,
Y.
,
Nix
,
W.
, and
Hutchinson
,
J.
,
1999
, “
Mechanism-Based Strain Gradient Plasticity–I. Theory
,”
J. Mech. Phys. Solids.
,
47
(
6
), pp.
1239
1263
.
48.
Huang
,
Y.
,
Gao
,
H.
,
Nix
,
W.
, and
Hutchinson
,
J.
,
2000
, “
Mechanism-Based Strain Gradient Plasticity—II. Analysis
,”
J. Mech. Phys. Solids.
,
48
(
1
), pp.
99
128
.
49.
Nye
,
J. F.
,
1953
, “
Some Geometrical Relations in Dislocated Crystals
,”
Acta. Metall.
,
1
(
2
), pp.
153
162
.
50.
Mindlin
,
R. D.
,
1965
, “
Second Gradient of Strain and Surface-Tension in Linear Elasticity
,”
Int. J. Solids. Struct.
,
1
(
4
), pp.
417
438
.
51.
Bažant
,
Z. P.
, and
Guo
,
Z.
,
2002
, “
Size Effect and Asymptotic Matching Approximations in Strain-Gradient Theories of Micro-Scale Plasticity
,”
Int. J. Solids. Struct.
,
39
(
21–22
), pp.
5633
5657
.
52.
Jirásek
,
M.
, and
Bažant
,
Z. P.
,
2002
, “
Inelastic Analysis of Structures
,”
John Wiley & Sons
,
New York
.
53.
Fleck
,
N.
,
Muller
,
G.
,
Ashby
,
M. F.
, and
Hutchinson
,
J. W.
,
1994
, “
Strain Gradient Plasticity: Theory and Experiment
,”
Acta. Metall. Mater.
,
42
(
2
), pp.
475
487
.