State of art in regularization methods for numerical analysis of structures with softening

Abstract

This paper provides an extensive review of popular regularization methods utilized in numerical models to stabilize the structural response of materials exhibiting significant softening. The necessity for regularization is highlighted in cases of material softening, which is attributed to the loss of ellipticity in the governing differential equations. It discusses the advantages and disadvantages of the regularization methods most commonly employed in the scientific community. Furthermore, the paper highlights recent advancements, particularly in defining internal length within nonlocal models and characteristic element length in fracture energy regularization methods, as alternative solutions to address the limitations inherent in traditional approaches.

Introduction

Quasi-brittle materials, such as concrete, rocks, and composites, exhibit inherent heterogeneity that results in complex mechanical behavior. Even in their unloaded state, these materials contain defects, including micro-cracks and voids. In the initial stages of the fracturing process, both pre-existing and newly formed micro-cracks expand, leading to the formation of a distributed network of micro-cracks as illustrated in Figure 1(a). This network occurs within what is termed the fracture process zone (FPZ), where micro-crack development predominantly takes place. As the loading continues, the propagation and coalescence of these micro-cracks are confined to a progressively narrower band, depicted in Figure 1(b), eventually resulting in the emergence of a visible, macroscopic crack while many other micro-cracks are deactivated during the process of elastic unloading, as shown in Figure 1(c). This evolution of the FPZ leads to a decrease in material stiffness and introduces a softening behavior observable at the macro-scale, which is identified as mechanical damage. Such phenomena can be accounted for by models founded on continuum damage mechanics (CDM).

Figure 1.

A schematic showing evolution of macro-crack from a diffused network of micro-cracks (Sarkar et al., 2019).

Originating from the seminal works of Kachanov (1958) and Rabotnov (1969), damage models incorporate internal variables to characterize the density and orientation of micro-defects. The simplest form is the isotropic damage model, wherein the damaged stiffness tensor is considered a scalar multiple of the initial elastic stiffness tensor (Mazars, 1984; Mazars et al., 1990). Furthermore, the integration of damage models with plasticity or inelastic strain mechanisms, as explored by Grassl and Jirásek (2006) and Nechnech et al. (2002), addresses the inherent brittleness of pure damage models by accounting for permanent strain. More refined theories account for damage anisotropy using vectors or tensors, as detailed in the works of Desmorat et al. (2007), Halm and Dragon (1998), and Krajcinovic and Fonseka (1981), which provide a more accurate depiction of damage by acknowledging its directional properties.

In structures subjected to extreme loading conditions, strain increments often concentrate in narrow zones, whereas the remaining part of the structure undergoes unloading. This phenomenon, known as strain localization, results in damage being confined to the minimal possible volume. As a result, during the softening phase, the governing differential equations lose ellipticity in a static formulation or hyperbolicity in a dynamic formulation, rendering the boundary value problem ill-posed. Consequently, models based on classical continuum theory are unable to objectively predict post-peak results and localized failure patterns accurately (Bažant, 1976; De Borst et al., 1993). These deficiencies manifest themselves in numerical calculations by displaying a pronounced sensitivity to spatial discretization factors such as the orientation and size of the finite-element (FE) mesh during softening (Bažant and Oh, 1983; Cervera et al., 2010b). Specifically, the dissipated energy and force-displacement response of the numerical model are highly dependent on the FE mesh density (Bažant and Jirásek, 2002; Oliver et al., 2012). Furthermore, the direction of strain localization bands tends to align with the FE mesh lines, affecting the propagation patterns and challenging the objective representation of failure (Cervera and Chiumenti, 2006; Cervera et al., 2010a).

In recent years, various regularization techniques have been developed to enhance the objectivity of continuum damage models (Bažant and Jirásek, 2002). The softening phenomenon is modeled as dependent on the energy dissipation across a representative volume, leading to the adoption of nonlocal formulations (Poh and Swaddiwudhipong, 2009b). Nonlocal theories (Kröner, 1967) propose that the state of any material point within a domain is influenced by the entire domain, thereby ensuring a more comprehensive response to damage. To mitigate strain localization into an arbitrarily small volume, those variables that cause strain-softening are regularized as nonlocal through an internal length scale parameter (Bažant and Lin, 1988a, 1988b). This regularization is achieved using either integral (Pijaudier-Cabot and Bažant, 1987) or gradient (Peerlings et al., 1996) forms, each designed to capture the spatial distribution of damage more effectively. Interestingly, the gradient form, derived from the integral form with the Taylor series, demonstrates that both approaches are equivalent when specific weight functions are applied, highlighting a fundamental consistency within nonlocal theories (Jirásek, 2007b).

Nonlocal damage models have been demonstrated to possess outstanding regularization capabilities, offering mesh-size independent solutions and avoiding zero energy dissipation. However, both standard nonlocal integral and gradient-enhanced damage models still face various numerical challenges. Issues such as incorrect predictions of damage initiation and the emergence of spurious damage zones at the final softening stage have been reported in Geers et al. (1998) and Simone et al. (2004). Additionally, the inaccurate treatment of boundary conditions has been shown to lead to incorrect energy dissipation (Bažant et al., 2010; Jirásek et al., 2004). Furthermore, size effects are often misrepresented, contributing to the limitations of models (Jirásek et al., 2004). These issues can often be attributed to the constant internal length parameter used in standard nonlocal damage models. Consequently, research has explored evolving interaction domains, initiated by pioneering work (Geers et al., 1998), to mitigate spurious energetic interactions and ensure localized damage bandwidths during the final failure stages.

In contrast to the significant computational costs associated with the nonlocal averaging process required across the entire domain in nonlocal damage models, the crack band model (Bažant and Oh, 1983) offers a simpler and more easily implemented remedy. This model necessitates estimating the crack bandwidth, correlating with the element size within the localization zone in numerical models. Subsequently, the stress–strain softening relations for these elements are adjusted to preserve the fracture energy dissipation during softening, thereby ensuring an objective response from the global FE model. This concept underlies the term “fracture energy regularization” (Jirásek and Bauer, 2012). However, it is crucial to note that this approach represents a partial regularization technique. Mathematically, the boundary value problem remains ill-posed, leading to some inherent limitations (Lopes et al., 2020).

Viscous regularization introduces a time-dependent viscous term into the constitutive equations, offering an alternative approach to material behavior modeling. This term adds a rate-dependent behavior to the material, ensuring the ellipticity of the incremental equilibrium equations and addressing the issue of strain localization into infinitesimally small regions (Geers et al., 1994; Needleman, 1988). Incorporating viscosity renders the softening behavior dependent on the time step, ensuring that the tangent stiffness matrix stays positive definite for adequately small time increments (Langenfeld et al., 2018). Consequently, viscous regularization allows for the recovery of a well-posed boundary value problem, enhancing solution stability and reducing its sensitivity to the discretization scheme. However, choosing an appropriate viscosity parameter is critical.

Despite the demonstrated effectiveness of previously discussed regularization techniques in addressing ill-posed material models for (quasi)-brittle damage, certain limitations persist. Furthermore, there are few comprehensive review works on this subject. Bažant and Jirásek (2002) provided a detailed review of various regularization techniques, though this work is now over 20 years old. Simone (2007) and Jirásek (2007b), respectively, reviewed gradient and integral nonlocal damage models. Jirásek and Bauer (2012) offered a comprehensive investigation into numerical aspects of fracture energy regularization. All these contributions date back more than 10 years. Recently, Liu et al. (2022) summarized some regularization approaches in dealing with mesh dependency. Moreover, Langenfeld et al. (2022) presented a concise review, offering unified comparisons of these three regularization techniques within a variational framework. However, this review did not cover the latest modifications to these methods, essential for achieving mesh-independent results without spurious phenomena.

Given this context, this work is primarily motivated to offer a detailed review of recent modifications applied to these regularization techniques. Accordingly, the issues of non-objectivity and mesh dependency in standard strain-softening continuum models are initially introduced. Subsequently, commonly used standard regularization techniques, detailing their advantages and disadvantages, are presented. Following that, modifications of conventional regularization techniques are explored, with a special focus on the definition of internal length in nonlocal damage models and characteristic element length in fracture energy regularization. Finally, some concluding remarks regarding future applications of these techniques are provided.

Non-objectivity in the presence of softening

A continuum model for micro-cracking in quasi-brittle materials leads inevitably to strain softening, a phenomenon that the stresses in a strain-softening material reduce gradually with the increasing of strains due to the growing of micro-cracks and finally reduce to zero as the micro-cracks develop into a major visible crack. To capture this process, CDM introduces damage variables that quantify the degradation of material stiffness caused by micro-cracking. For instance, an isotropic damage variable $d$ can be defined to represent material degradation, where $d=0$ indicates an intact material, and $d=1$ signifies a fully cracked material. The simplest elasticity-based damage model (Lemaitre and Chaboche, 1994) can thus be expressed as:

$$\mathit{\sigma }=(1-d)\mathbb{C}:\mathit{\epsilon }$$

(1)

where $\mathit{\sigma }$ and $\mathit{\epsilon }$ are the stress and strain tensors, respectively. $\mathbb{C}$ denotes the elasticity tensor.

The governing equations for a quasi-static problem are given by the equilibrium condition:

$$\mathrm{\nabla }\cdot \mathit{\sigma }+\mathbf{f}=\mathbf{0}$$

(2)

where $\mathbf{f}$ represents the body force per unit volume. Substituting the constitutive relationship into the equilibrium equation yields the following second-order partial differential equation (PDE) in terms of the displacement field $\mathbf{u}$ :

$$(1-d)\mathbb{C}:{\mathrm{\nabla }}^2\mathbf{u}-\mathrm{\nabla }d\cdot (\mathbb{C}:\mathrm{\nabla }\mathbf{u})+\mathbf{f}=\mathbf{0}$$

(3)

For a given damage field $d<1$ , the system governed by this PDE remains elliptic, and the displacement field can be determined uniquely. However, when the material reaches a fully cracked state ( $d=1$ ), both terms involving $d$ in the PDE vanish, causing the system to lose the ellipticity. This degeneration renders the boundary value problem ill-posed (Peerlings et al., 2002), leading to a lack of unique solutions within the softening region.

Jirásek (2004, 2007a) analyzed the pathological characteristics of strain-softening solutions in CDM, highlighting the following issues: (i) an infinitely small softening region; (ii) a load-displacement diagram that consistently exhibits snapback; and (iii) zero total energy dissipation during the softening phase.

The non-objectivity associated with strain softening can be effectively illustrated by the example of a simple bar subjected to uniaxial tension. As depicted in Figure 2(a), a straight bar with a constant cross-sectional area $A$ and total length $L$ undergoes tension due to an applied displacement $u$ . The material is modeled to exhibit linear strain softening post-peak stress $f_{\mathrm{t}}$ , as illustrated in Figure 2(b). The strain at the elastic limit, corresponding to peak stress, is denoted as ${\epsilon }_0$ , while the ultimate strain is represented by ${\epsilon }_{\mathrm{u}}$ .

Figure 2.

(a) A uniform bar under uniaxial tension and (b) constitutive law with linear softening.

Before reaching peak stress, the strain is uniquely determined by the stress. The strain field is uniform along the bar since the stress field is uniform according to the static equilibrium. The response of the bar is linear elastic up to a displacement $u_0=L{\epsilon }_0$ , at which point the reaction force reaches its maximum value of $F_0=A{f}_{\mathrm{t}}$ . Subsequently, the resistance of the bar begins to decrease. Across each cross-section, stress may decrease with increasing strain (indicating softening) or with decreasing strain (indicating elastic unloading). As a result, the uniqueness of the solution is lost. Considering the static equilibrium equation, for any given stress between 0 and $f_{\mathrm{t}}$ , two strain values, ${\epsilon }_{\mathrm{s}}$ and ${\epsilon }_{\mathrm{u}}$ , exist. Consequently, the response of the bar is expressed as follows:

$$u=L_{\mathrm{s}}{\epsilon }_{\mathrm{s}}+L_{\mathrm{u}}{\epsilon }_{\mathrm{u}}\mathrm{with}{L}_{\mathrm{s}}+L_{\mathrm{u}}=L$$

(4)

where $L_{\mathrm{s}}$ and $L_{\mathrm{u}}$ are the cumulative lengths of the softening and unloading regions, respectively.

A primary challenge presented by equation (4) is determining the length $L_{\mathrm{s}}$ , which may vary discontinuously between 0 and $L$ . This uncertainty regarding the softening region suggests that no unique and objective solutions can be obtained after the bar reaches peak load. Consequently, in the context of softening, the boundary value problem may become ill-posed, as the governing equations permit an infinite number of solutions.

In reality, material properties are not perfectly uniform, evidenced by reduced strength in localized regions. When stress reaches this reduced strength, softening initiates within the weakened area, while regions outside it undergo elastic unloading. Consequently, the size of the softening region is determined by the size of these imperfect areas. However, the softening region is often infinitesimally small due to the minor size of the imperfections. As a result, the entire degradation of material properties becomes concentrated within this small volume of imperfections. Such localization results in a perfectly brittle macroscopic fracture behavior, deviating from the anticipated gradual stiffness reduction.

Mesh dependency and sensitivity in FE method-based numerical analyses are clearly illustrated through the given example of a tension bar shown in Figure 3. Consider the bar subdivided into $m$ elements, with a minor defect, such as a reduced cross-section, positioned at the midpoint. The element with this imperfection is likely to reach the strength threshold, $f_{\mathrm{t}}$ , more rapidly than the intact elements, initiating material failure at this point. Subsequently, the stresses in neighboring elements cannot reach the tensile strength and will decrease elastically. Consequently, strain distribution concentrates into a narrow band, the width of which is significantly influenced by the size of the element, as depicted in Figure 3(a). The fracture zone is determined by the spatial resolution of the discretization, which approaches zero as the mesh becomes infinitely fine. Similarly, the energy dissipated by the fracture trends toward zero, a result that is non-physical.

Figure 3.

Response of an imperfect bar: (a) strain localization and (b) load-displacement diagrams due to different meshes.

Since only one element undergoes softening while the others are unloading, the ultimate displacement, corresponding to the zero reaction force, can be expressed as:

$$u_{\mathrm{f}}=L_{\mathrm{s}}{\epsilon }_{\mathrm{f}}=\frac{L}{m}{\epsilon }_{\mathrm{f}}$$

(5)

Subsequently, the load-displacement diagram, significantly influenced by the element size, is illustrated in Figure 3(b). As the number of elements increases and their size decreases, the load-displacement curve consistently demonstrates snap-back behavior. For an infinite number of elements, the post-peak curve may coincide with the elastic loading curve. From a numerical perspective, ill-posedness is thus manifested by the pathological sensitivity of numerical results to the FE size, which determines the size of the softening region.

In general, standard CDM often yields mesh-size dependent results, which originates from the loss of ellipticity in the governing PDEs (De Borst et al., 1993). Achieving objective results in FE analyses becomes challenging due to the computed fracture energy and the corresponding structural response being influenced by the mesh size. Upon mesh refinement, a physically unrealistic situation emerges: an infinitesimally small element size leads to negligible energy dissipation and crack formation (Bažant and Jirásek, 2002). Beyond mesh-size dependency, spurious crack trajectories have been observed, indicating that numerical results can vary with the orientation of the FE mesh. This phenomenon, which may be attributed to a direct consequence of spatial discretization inadequacies, is referred to as mesh-bias dependency (Cervera and Chiumenti, 2006; Cervera et al., 2010a; Jirásek and Grassl, 2008).

Additionally, further mesh refinement in simulations designed for multiaxial conditions often leads to convergence issues in equilibrium iterations at peak loads (Jirásek, 2004). This issue arises because refining the mesh increases the number of potential equilibrium states. Consequently, the complexity of determining the correct combination of loading and unloading at individual Gauss points escalates. During the iterative solution process, the algorithm must “select” from among these equilibrium states and tends to choose a different state in each subsequent iteration. This ambiguity ultimately causes the iterative solution process to diverge (Borst, 2001).

Conventional regularization techniques

The mesh sensitivity in local damage models, caused by the loss of ellipticity, stems from an inherent mathematical nature rather than numerical ones. The scientific community has proposed several regularization techniques to restore objectivity in numerical analyses and render the boundary value problem well-posed with reliable outcomes. This section will overview commonly employed approaches, some of which have been implemented in commercial software. The core idea behind these methods is the introduction of a localization limiter to prevent strain localization into a line, thereby avoiding zero energy dissipation.

The first two approaches involve nonlocal models, which are classified into integral and gradient types based on their formulation. The third method, fracture energy regularization, establishes a direct correlation between fracture energy and mesh size for numerical analyses. The final approach is viscous regularization, which incorporates viscosity into the damage evolution.

Integral nonlocal model

The concept of “nonlocal” might have been first proposed in the 1960s by Eringen (1966) and Kröner (1967), introducing the idea that the stress at a point could depend not only on the state variables at that point but also on the distribution of state variables across the entire body. This concept was then proposed as a localization limiter by Bažant et al. (1984) and later applied to CDM by Pijaudier-Cabot and Bažant (1987). Nonlocal versions of a smeared crack model, a model with softening plasticity, and a microplane model were further proposed by Bažant and Lin (1988a, 1988b) and Bažant and Ožbolt (1990) to address the issue of strain localization. In integral nonlocal damage models (Pijaudier-Cabot and Bažant, 1987), the local variable $f(x)$ is replaced by its nonlocal counterpart $\overline{f}(\mathit{x})$ , obtained through spatial averaging over the surrounding domain $V$ with respect to a weight function ${\alpha }_0$ :

$$\overline{f}(\mathit{x})={\int }_V\alpha (\mathit{x},\mathit{\xi })f(\mathit{\xi })\mathrm{d}\mathit{\xi }$$

(6)

where $\alpha (\mathit{x},\mathit{\xi })$ represents a nonlocal weight function that is determined by the distance between the “receiver” point, $\mathit{x}$ , and a “source” point, $\mathit{\xi }$ , within the integration domain $V$ .

The nonlocal weight function, $\alpha (\mathit{x},\mathit{\xi })$ , must be a non-negative function that decreases monotonically as the distance between point $\mathit{x}$ and point $\mathit{\xi }$ increases. A notable challenge arises near the boundaries of the domain: the domain of influence contributing to the nonlocal average at a point close to the boundary is smaller than it is for points located further from the boundary. This discrepancy results in lower nonlocal values at the boundary compared to the interior of the body. To address this issue, the weight function is usually rescaled in the vicinity of a boundary to ensure a uniform field. This can be achieved by setting:

$$\alpha (\mathit{x},\mathit{\xi })=\frac{{\alpha }_0(\mathit{x},\mathit{\xi })}{{\int }_V{\alpha }_0(\mathit{x},\mathit{\xi })\mathrm{d}\mathit{\xi }}$$

(7)

with a requirement known as the reproducibility of order zero (Jirásek, 2007b):

$${\int }_V\alpha (\mathit{x},\mathit{\xi })\mathrm{d}\mathit{\xi }=1$$

(8)

where ${\alpha }_0$ serves as the basic weighting function that illustrates the diminishing effect of nonlocal interactions as the distance $r=|\mathit{x}-\mathit{\xi }|$ increases. The most widely employed nonlocal weight function is the Gaussian distribution function:

$${\alpha }_0(r)=\exp (-\frac{r^2}{2l^2})$$

(9)

where $l$ is an internal parameter signifying the characteristic length scale of nonlocal averaging. This scale can be associated with the FPZ length introduced by Hillerborg et al. (1976), defined as:

$$\$l=\frac{E{G}_{\mathrm{f}}}{f_{\mathrm{t}}^2}\$$$

(10)

where $E$ is the elasticity modulus, $G_{\mathrm{f}}$ is the fracture energy, and $f_{\mathrm{t}}$ is tensile strength.

Investigations by Bažant and Pijaudier-Cabot (1989) reveal that equation (10) typically yields a larger value for $l$ . As a practical alternative, the characteristic length was linked to the material heterogeneity and empirically suggested to be approximately 2.7 times the maximum aggregate size for concrete. Furthermore, Le Bellégo et al. (2003) recommended that the FE size within the process zone should be no larger than one-third of the characteristic length. Consequently, the characteristic length in FE models can be determined as three times the FE size. Additionally, Le Bellégo et al. (2003) emphasized that determining an accurate characteristic length requires calibration to achieve the best fit with experimental data.

In addition, other functional forms for ${\alpha }_0$ include the bell-shaped function (Bažant and Ožbolt, 1990) and the Green distribution function (Grassl et al., 2014), which are expressed as equations (11) and (12), respectively.

$${\alpha }_0(r)={⟨1-\frac{r^2}{R^2}⟩}^2$$

(11)

$${\alpha }_0(r)=\exp (-\frac{r}{l})$$

(12)

where $⟨\cdot ⟩$ is the Macaulay bracket that intends to extract the positive part; $R$ denotes the interaction radius, related to the internal length scale, and represents the maximum distance at which point $\mathit{\xi }$ influences the nonlocal average at point $\mathit{x}$ .

To prevent the concentration of local strain at a central point and to ensure that nonlocal plastic strain is distributed over a larger region during plastic softening, Galavi and Schweiger (2010) introduced a new weight function known as the G&S function (Su et al., 2023). This function is mathematically expressed as follows:

$${\alpha }_0(r)=\frac{r}{l}\exp (-\frac{r^2}{l^2})$$

(13)

Equations (9) to (13), with internal length $l=$ 1, are illustrated in Figure 4(a) for comparison. The bell-shaped weight function exhibits bounded support, signifying that the interaction between two points vanishes when their distance exceeds $R$ , which is $\sqrt{2}l$ in Figure 4(a). In contrast, both the Gaussian and Green weight functions are unbounded, indicating their interaction radius $R=\mathrm{\infty }$ . The G&S function is different from others in that the nonlocal variable at a specific point is independent of the local variable at that point.

Figure 4.

(a) Comparison of different nonlocal weight functions ${\alpha }_0$ and (b) distribution of local and different nonlocal variables.

Although differences exist in the contribution from a neighboring point at the same distance across different weight functions, they all enable the smoothing of the spatial distribution of the local variable, even after localization onset, as exemplified in Figure 4(b). In Figure 4(b), local variables along a bar of length 10 are modeled to follow a Gaussian distribution. Equations (9) and (12) yield the similar regularized results. Though equation (11) presents a slightly different outcome, it can achieve similar results as the other two weight functions by adjusting $R$ . Equation (13) is expected to exhibit different nonlocal performance due to its consideration of a broader nonlocal distribution.

The selection of nonlocal variables is crucial for incorporating nonlocality within damaged domains while allowing for a return to the standard local elastic continuum model in cases where material behavior remains within the elastic range (Jirásek, 1998). Bažant and Pijaudier-Cabot (1988) proposed that a nonlocal variable should be defined either as a weighted average of the damage variable itself or of a strain measure. A comparative study, conducted in Jirásek (1998), reveals only minor differences between results obtained using nonlocal damage energy release rates and nonlocal equivalent strain. However, formulations based on the damage variable can lead to spurious stress-locking effects during later stages of softening, a limitation that hinders the modeling of stress-free macroscopic cracks (Jirásek, 1998, 2007b). An innovative alternative, the averaging of the local displacement field, was first discussed in Rodrıguez-Ferran et al. (2004). However, its practical implementation faces challenges due to the need for specific modifications of the averaging operator, as further detailed in Jirásek and Marfia (2005) and Rodríguez-Ferran et al. (2005). Although the use of nonlocal damage energy release rates is physically appealing, the nonlocal strain variable is more commonly adopted due to its versatility and applicability across a wider range of constitutive laws. Commonly utilized strain variables in the literature include the energy norm-based equivalent strain (Jirásek and Patzák, 2002), Von Mises equivalent strain (Peerlings et al., 2001), modified Von Mises equivalent strain (De Vree et al., 1995), Rankine equivalent strain (Jirásek, 2004), and Mazars equivalent strain (Mazars and Pijaudier-Cabot, 1989). Furthermore, two strain measures can be decomposed to specifically address tension and compression, respectively (Comi, 2001; He et al., 2015).

Some careful treatments are required when applying integral-type nonlocal damage models. Currently, there is no consensus on the appropriate adjustment of averaging near the physical boundaries of a body. A commonly adopted solution involves modifying and normalizing the weight function, as shown in equation (7), to comply with equation (8). However, this arbitrary treatment results in the asymmetry of the weight function, that is, $\alpha (\mathit{x},\mathit{\xi })\ne \alpha (\mathit{\xi },\mathit{x})$ , and consequently leading to the non-symmetry of the tangent operators (Bažant and Pijaudier-Cabot, 1988). Furthermore, the change of nonlocal interactions near the boundary of the solid suggests that the constitutive formulation (Pijaudier-Cabot and Bažant, 1987) may not adhere to thermodynamic principles.

More critical numerical defects include incorrect damage initiation and improper damage propagation behaviors (Krayani et al., 2009; Pijaudier-Cabot and Dufour, 2010; Simone et al., 2004), which are non-physcial. Notably, the development of a spurious damage zone during the final stage of softening prevents the formulation of a macroscopic discontinuous crack. This issue arises from the unexpected influence of highly damaged zones on their surrounding undamaged area. Moreover, non-physical interactions between material points across damaged bands represent another drawback of this regularization technique.

Another significant drawback concerns computational costs. The numerical implementation of the nonlocal integral formulation often requires a global averaging procedure. A frequently employed approach is the traversal algorithm whose time complexity is $O(N^2)$ (Sarkar et al., 2019), with $N$ denoting the number of data points. This quadratic growth underscores the unsuitability of traditional nonlocal models for large-scale problems.

To mitigate the computational challenge, Jirásek (2007b) proposed evaluating weight coefficients once during the initial computation and storing them for subsequent use. However, this approach neglects potential changes in the interaction domain due to structural deformations. Other studies (Summersgill et al., 2017b; Tang et al., 2021) have proposed limiting nonlocal interactions to a circular or spherical radius. While this method reduces the computational domain, the radius selection significantly impacts computational costs, and the process of identifying nonlocal neighbors by evaluating pairwise distances between Gaussian points (GPs) remains a bottleneck. A more efficient approach based on the quadtree (in two-dimensional (2D)) or octree (in three-dimensional (3D)) algorithm was suggested by Jirásek (2004) and later successfully implemented by Lu et al. (2022, 2024). This method reduces the time complexity to $O(\text{log}N)$ , achieving up to a 100 times speedup at the integration point level compared to the traversal algorithm.

Gradient nonlocal model

Instead of addressing spatial interactions through integrals, the gradient-type nonlocal model integrates the gradients of internal variables into the constitutive relations (Bažant, 1984; Lasry and Belytschko, 1988; Schreyer and Chen, 1986). This approach was first introduced in the domain of damage mechanics by Peerlings et al. (1996). A gradient theory emerges from the nonlocal formulation given in equation (6) by employing a Taylor series expansion to approximate the local quantity as follows:

$$\begin{array}{rl}f(\mathit{x}+\mathit{\xi })& =f(\mathit{x})+\mathrm{\nabla }f(\mathit{x})\cdot \frac{\mathrm{\partial }\mathit{\xi }}{\mathrm{\partial }\mathit{x}}+\frac{1}{2!}{\mathrm{\nabla }}^{2}f(\mathit{x})\cdot \frac{{\mathrm{\partial }}^2\mathit{\xi }}{{\mathrm{\partial }}^2\mathit{x}}\\ & +\frac{1}{3!}{\mathrm{\nabla }}^{3}f(\mathit{x})\cdot \frac{{\mathrm{\partial }}^3\mathit{\xi }}{{\mathrm{\partial }}^3\mathit{x}}+\cdots \end{array}$$

(14)

Substituting equation (14) into equation (6) yields the differential form of the nonlocal quantity:

$$\overline{f}(\mathit{x})=f(\mathit{x})+c{\mathrm{\nabla }}^{2}f(\mathit{x})+d{\mathrm{\nabla }}^{4}f(\mathit{x})+\cdots$$

(15)

in which ${\mathrm{\nabla }}^2$ is the Laplacian operator, $c$ and $d$ are parameters associated with internal length. By neglecting higher-order gradient terms in equation (15), the explicit second-order gradient enhancement formulation is derived:

$$\overline{f}(\mathit{x})=f(\mathit{x})+c{\mathrm{\nabla }}^{2}f(\mathit{x})$$

(16)

From a mathematical perspective, the explicit gradient model, which incorporates only the local quantity and its second-order derivatives, remains inherently local. This model retains its local character unless the truncation of higher-order terms in equation (15) is avoided (Peerlings et al., 2001). Furthermore, the numerical implementation of this model within an FE framework poses challenges due to the requirement for $C^1$ -continuity in displacement approximation. This requirement arises from the incorporation of second derivatives of local variables in equation (16), complicating the implementation process. The limitations inherent to the explicit model are effectively addressed by its implicit counterpart. This involves an additional differentiation of equation (16), resulting in:

$${\mathrm{\nabla }}^{2}f(\mathit{x})={\mathrm{\nabla }}^2\overline{f}(\mathit{x})-c{\mathrm{\nabla }}^{4}f(\mathit{x})-d{\mathrm{\nabla }}^{6}f(\mathit{x})+\cdots$$

(17)

substitution of which into equation (15) and neglecting the terms higher than the second order lead to a second-order implicit gradient formulation:

$$\overline{f}(\mathit{x})-c{\mathrm{\nabla }}^2\overline{f}(\mathit{x})=f(\mathit{x})$$

(18)

in which the nonlocal variable is defined as a function of the corresponding local quantity, allowing for the direct application of $C^0$ -continuous displacement interpolation.

To solve the Helmholtz-type differential equation presented in equation (18), a homogeneous Neumann boundary condition is imposed on the nonlocal quantity across the entire physical surface:

$$\mathbf{n}\cdot \mathbf{\nabla }\overline{f}=0$$

(19)

where $\mathbf{n}$ denotes the external normal unit vector to the surface. Equation (19) presents the vanishing normal strain gradient at a free edge. Nevertheless, the physical interpretation of this natural boundary condition remains ambiguous (Jirásek, 2007b).

In contrast to equation (16), the derivatives of the nonlocal variable in equation (18) include an infinite series of higher-order derivatives of the local quantity implicitly (Peerlings et al., 2001). This characteristic signifies that the implicit gradient model is a truly nonlocal model. Additionally, the truncation of higher-order terms in the implicit gradient model exerts a subtle influence on its performance (Askes et al., 2000). Moreover, studies by Peerlings et al. (2001) and Simone (2007) reveal that the implicit gradient model is equivalent to the integral model when the Green-type weight function is utilized. Performance comparisons in Askes et al. (2000) demonstrate that the implicit gradient and integral models possess similar regularization capabilities, with both significantly outperforming the explicit gradient model.

Numerical challenges similar to those encountered in integral nonlocal models also require careful consideration in gradient nonlocal models. First, the boundary condition as defined in equation (19) requires the normal component of the gradient of the nonlocal variable to vanish at the boundary. This constraint causes nonlocal interactions near an axis of symmetry to mirror those near the physical boundary of the solid (Krayani et al., 2009).

Additionally, stress oscillation in the case of low-order elements within gradient models has been reported in Simone et al. (2003). Notably, this oscillation intensifies with mesh refinement, indicating the significance of employing integral nonlocal smoothing techniques in certain scenarios (Negi and Kumar, 2019; Nguyen et al., 2018).

Furthermore, similar numerical deficiencies observed in integral nonlocal models—namely, incorrect damage initiation, spurious damage distribution, and stress-locking effects—are also present in gradient nonlocal models (Geers et al., 1998; Simone et al., 2004; Xue et al., 2024).

Fracture energy regularization

As discussed in the previous section, equation (4) admits an infinite number of solutions, among which strain softening localizes in specific regions while the remainder of the material undergoes unloading. Particularly, some solutions exhibit damage growth localized within an interval of length $h_{\mathrm{b}}$ , which can be interpreted as a crack bandwidth. This bandwidth is influenced by the size and spacing of material heterogeneities, as described by Bažant and Pijaudier-Cabot (1989). In numerical FE simulations, the damaged zone typically localizes to a single layer of FEs, leading to the phenomenon where finer mesh discretization results in negligible energy dissipation (Bažant and Jirásek, 2002; Oliver et al., 2012).

To prevent the issue of diminishing energy dissipation with infinitesimally small element sizes, a fracture energy regularization technique based on the crack band model (Bažant and Oh, 1983; Comi and Perego, 2001) is often employed for practical engineering computations. This technique combines the cohesive crack model and classical continuum mechanics (Kurumatani et al., 2016; Xu and Waas, 2016). The fundamental concept involves adjusting the softening branch of the stress–strain curve relative to the crack bandwidth to ensure consistent energy dissipation per unit area across varying element sizes. In other words, the stress–strain relationship at each material point is no longer a direct representation of material properties but varies according to mesh size. As a result, this approach ensures that global load-displacement outcomes and fracture energy are matched to experimental data.

The adjustment of softening branches is illustrated in Figure 5, which shows the derivation of an effective continuum stress–strain law from the traction–separation law of a crack. Within the cohesive crack model framework (Hillerborg et al., 1976), the fracture energy $G_{\mathrm{f}}$ , defined as the energy dissipation per unit area of a surface during failure, is represented by the area beneath the traction–separation curve. The expression of fracture energy is:

$$G_{\mathrm{f}}={\int }_0^{\mathrm{\infty }}f\mathrm{d}w$$

(20)

where $f$ represents the cohesive force, which is directly related to the crack opening displacement (COD), $w$ .

Figure 5.

Damage zone and fracture energy regularization (He et al., 2019).

The crack band approach assumes that during the failure process, a damage zone of limited width undergoes a softening period, while the remaining areas experience unloading (Bažant and Oh, 1983). Supposing that a single element accommodates a single crack or localization zone, the COD can be expressed as follows:

$$w=\epsilon h-{\epsilon }_{0}h=h(\epsilon -{\epsilon }_0)$$

(21)

where $\epsilon$ represents the measured strain and ${\epsilon }_0$ denotes the elastic limit strain. $h$ signifies the crack bandwidth in which damage is evaluated. This bandwidth corresponds to the size of a single FE in FE analysis. Importantly, $h$ works as a converter of strain $\epsilon$ from damage models and COD $w$ from a cohesive crack model.

In damage models, the cohesive force $f$ acting on the fracture surface is equated to the stress $\sigma$ . By substituting equation (21) into equation (20), the fracture energy can be expressed as:

$$G_{\mathrm{f}}=h{g}_{\mathrm{f}}with{g}_{\mathrm{f}}={\int }_{{\epsilon }_0}^{\mathrm{\infty }}\sigma \mathrm{d}\epsilon$$

(22)

where $g_{\mathrm{f}}$ is the volumetric fracture energy, defined as the energy dissipation per under volume. This approach only requires uniaxial constitutive laws, allowing for the application in multi-dimensional problems by replacing the unidirectional strain with equivalent strain measurement.

As illustrated in Figure 5, the volumetric fracture energy $g_{\mathrm{f}}$ is represented by the shaded area beneath the stress–strain curve. The softening branch gets rescaled through the introduction of $h$ , and the global fracture energy is preserved regardless of mesh size. The determination of $h$ is thus linked to the characteristics of FEs. For instance, Bažant and Oh (1983) initially assumed $h$ to be the square root of the area of quadrilateral elements. Subsequently, Rots et al. (1988) refined this assumption by incorporating considerations of triangular elements. These definitions can be further extended to 3D FEs (Kurumatani et al., 2016), and the resulting calculations for both 2D and 3D elements are summarized as follows:

$$\text{Triangular element}:{\mathit{h}}_{\mathit{e}}=\sqrt{2{\mathit{A}}_{\mathit{e}}}$$

(23a)

$$\text{Tetrahedral element}:{\mathit{h}}_{\mathit{e}}=(12{\mathit{V}}_{\mathit{e}}{)}^{1/3}$$

(23b)

$$\text{Quadrilateral element}:{\mathit{h}}_{\mathit{e}}=\sqrt{{\mathit{A}}_{\mathit{e}}}$$

(23c)

$$\text{Hexahedral element}:{\mathit{h}}_{\mathit{e}}={{\mathit{V}}_{\mathit{e}}}^{1/3}$$

(23d)

where $h_e$ is the characteristic length of an individual element. $A_e$ represents the area of a triangular or quadrilateral element. $V_e$ denotes the volume of a tetrahedral or hexahedral element.

The fracture energy regularization technique offers significant advantages, notably that the formulation remains inherently local and its implementation is simpler than that of previous nonlocal enhancement formulations. Additionally, the boundary conditions established by this approach are both clear and physically justified, contributing to its robustness. This regularization method has been successfully applied in the development of new damage models (Arruda et al., 2022; Bui and Tran, 2022), enabling the achievement of mesh-size independent results.

However, the definition of the characteristic element length is crucial for the effectiveness of the fracture energy regularization technique. The original definition based on the element area in Bažant and Oh (1983) requires square or cubic element mesh refinement (Comi and Perego, 2001). Without meeting this requirement may result in continued mesh sensitivity. Additionally, constitutive-level snapback may arise either from overly large FEs or insufficiently fractured energy (Bažant, 2002; Bažant and Oh, 1983). Both excessively fine and overly coarse meshes can lead to inaccurate crack distribution. The width of the FPZ is fixed due to the premise of accommodating a single crack per element (Červenka et al., 2018). Moreover, this single-element crack band assumption requires uniform damage distribution across the band (Gorgogianni et al., 2020).

From a mathematical standpoint, the governing equations of the fracture energy regularization approach remain ill-posed, offering only a partial regularization of the problem. The mesh-induced directional bias and corresponding incorrect crack trajectories still exist (Jirásek and Grassl, 2008). Furthermore, damage localization into a zero-thickness layer with finer mesh discretization remains an unresolved challenge (Jirásek and Bauer, 2012).

Viscous regularization

Under dynamic loading conditions, quasi-brittle materials exhibit a pronounced rate-dependent behavior, characterized by a significant increase in dynamic strength at high strain rates. This strain-rate effect can be attributed to both viscous effects and the microscopic influences of mass inertia (Häussler-Combe and Panteki, 2016). In the framework of CDM, this effect leads to the retard of damage evolution at high strain rates. To account for the micro-inertia effects, Häußler-Combe and Kitzig (2009) and Wang et al. (2019) introduced the second time derivative of a nonlocal variable into the gradient model. Furthermore, Davaze et al. (2021), Häussler-Combe and Kühn (2012), Häussler-Combe and Panteki (2016), and Rosenbusch et al. (2024) incorporated an additional damping term, corresponding to the first time derivative of the nonlocal variable, which represents micro-viscosity as discussed in Huynh and Abedi (2025). The strain-rate effect can also be included by the introduction of viscosity terms into the constitutive relationships (Pedersen et al., 2008), which can transform the micro-viscosity effects into the macroscopic scale (Häussler-Combe and Kühn, 2012; Häussler-Combe and Panteki, 2016).

The dynamic effect observed during cracking is expected to be influenced by inertial effects, viscosity, rate effects, and damping, which tend to be neglected in quasi-static loading scenarios. Employing these factors in damage models could enable the prediction of nonlocal effects through damage models, presenting a feasible strategy for regularization. Both real and artificial factors may serve as an expedient substitute for nonlocality.

The simplest approach involves directly employing classical Rayleigh damping, which is proportional to specified stiffness and mass, along with a given damping percentage for each material. However, caution is recommended, since it has been demonstrated (Bažant and Jirásek, 2002; de Borst and Duretz, 2020; Needleman, 1988) that incorporating any rate dependence along with inertial effects into the constitutive relation introduces an internal length scale, $l_{\mathrm{c}}$ , implicitly. This length scale is associated with the wave propagation velocity, $v$ , of an FE, which depends on Young’s modulus ( $E$ ), viscosity ( $\eta$ ), and mass density ( $\rho$ ), as shown in equation (24):

$$\{\begin{array}{c}v{l}_{\mathrm{c}}\rho =\eta \\ v=\sqrt{E/\rho }\end{array}\Rightarrow l_{\mathrm{c}}=\frac{\eta }{\sqrt{E/\rho }}$$

(24)

To avoid undesired numerical artifacts, it is advisable to use an artificial value for viscosity, rate dependence, or damping, rather than the material’s actual properties. This approach enables effective regularization that incorporates nonlocal effects without compromising the dynamic response.

Several researchers (Loret and Prevost, 1990; Needleman, 1988; Sluys and de Borst, 1992) have explored utilizing viscosity-based effects as localization limiters. These methods delay damage evolution by introducing a length scale into the damage law (Ladeveze, 1992; Ladevèze et al., 2000). However, as strain rates decrease, viscosity-induced nonlocal effects diminish. Consequently, the effectiveness of viscosity in serving as a stand-in for a nonlocal model is limited to a narrow range of time delays and strain rates, typically not exceeding one order of magnitude (Bažant and Jirásek, 2002).

For an FE with a given characteristic length, $l_{\mathrm{c}}$ , and applying the concept of a wavefront corresponding to this length, one can determine a characteristic time integration, $t_0$ , as shown in equation (25). In fast dynamic scenarios, such as high-speed impact, where the event duration is shorter than $t_0$ , regularization is unnecessary. In such cases, regardless of whether viscosity is considered, there is insufficient time for erroneous localization of strain-softening damage to emerge because of inertia effects.

$$\{\begin{array}{c}v{l}_{\mathrm{c}}\rho =\eta \\ v=\sqrt{E/\rho }\\ v=\frac{l_{\mathrm{c}}}{t_0}\end{array}\Rightarrow t_0=\frac{\eta }{E}$$

(25)

The inclusion of viscosity has been widely applied in visco-plasticity theory (de Borst and Duretz, 2020; Dias da Silva, 2004; Loret and Prevost, 1990) and damage models (Chaboche et al., 2001; Dubé et al., 1996; Faria et al., 1998) for localization problems based on the theory of Perzyna (1966) and Duvaut and Lions (1972). Prezyna and Duvaut–Lions models belong to the overstress theory in which the stress can exceed the rate-independent yield surface, thus not satisfying the consistency condition. Therefore, Wang et al. (1997) introduce a consistency viscoplastic model through a rate-dependent yield surface. Niazi et al. (2013) employed the consistency viscoplastic model to critically evaluate the reason for the existence of this regularization strategy for quasi-static problems.

The simplest method to implement regularization through viscosity for fast damage evolution relies on the Duvaut–Lions model (Duvaut and Lions, 1972), which separates non-viscous and viscous components of plastic strain ( ${\epsilon }_{vp}$ ) and the hardening parameters ( $\varkappa$ ):

$$\begin{array}{rl}\delta {\epsilon }_{vp}& =\frac{1}{\eta }{\left[D_e\right]}^{-1}(\left\{{\sigma }_v\right\}-\{\sigma \})\delta t\\ \delta {\varkappa }_v& =-\frac{1}{\eta }({\varkappa }_v-\varkappa )\delta t\end{array}$$

(26)

Subsequently, various researchers (Maimí et al., 2007) have adapted this formulation to explicitly account for the effects of viscous damage. The key idea is introducing a time dependence to the evolution of the internal variable responsible for ill-posedness (Faria et al., 1998; Geers et al., 1994; Needleman, 1988). Within this regularization framework for damage models, as shown in Figure 6, a viscous damage variable is established as follows:

$$\delta d_v=\frac{1}{\eta }(d-d_v)\times \delta t$$

(27)

where $\delta d_v$ represents the change of the viscous damage variable $d_v$ over time $\delta t$ . $\eta$ is the viscosity parameter which represents the relaxation time of the viscous system. $d$ denotes the damage variable evaluated in the non-viscous damage model.

Figure 6.

Rate-dependent softening model: (a) fracture domain and (b) rheological model referred to Bažant and Li (1997).

A numerical algorithm needs to be implemented for the solution of equation (27). A straightforward approach involves updating the damage variable using the backward Euler method as follows:

$$\delta d_v(t_0+\mathrm{\Delta }t)=\frac{\mathrm{\Delta }t}{\eta +\mathrm{\Delta }t}\delta d(t_0+\mathrm{\Delta }t)+\frac{\eta }{\eta +\mathrm{\Delta }t}\delta d_v(\mathrm{\Delta }t)$$

(28)

where $\mathrm{\Delta }t$ is the size of the time step.

The viscous formulation implicitly incorporates a length scale, presenting the similar regularization effects observed in nonlocal models (Bažant and Jirásek, 2002). This feature helps limit and thus slow down the evolution of damage within specific elements. Additionally, for sufficiently small time increments, this approach ensures that the tangent stiffness matrix of softening materials remains positive and definite, contributing to stable and realistic simulations. Utilizing viscous regularization with a minimal viscosity parameter $\eta$ enhances model convergence in the softening regime without compromising the accuracy of the numerical solution significantly (Lapczyk and Hurtado, 2007; Lopes et al., 2020). Moreover, this technique does not need any additional global discretization and only requires supplementary operations at the local level (Liu et al., 2022), indicating its straightforward implementation.

However, the selection of the viscosity coefficient $\eta$ is critical. An inappropriate value can result in either inadequate regularization, failing to prevent localization, or introducing excessive damping, masking the true brittle nature of the material. Furthermore, the nonlocal effects induced by viscosity diminish over time (Bažant and Jirásek, 2002). This implies that the long-term behavior of the viscous system gradually resembles that of its inviscid counterpart.

Modifications on nonlocal regularization technique

Nonlocal formulations act as efficient localization limiters and provide an objective description of strain localization (Pijaudier-Cabot and Benallal, 1993). Nonetheless, standard nonlocal models, including both integral-type and gradient-type, have been observed to exhibit several numerical artifacts:

Incorrect failure initiation and propagation: Geers et al. (1998) and Simone et al. (2004) have reported that the shift of maximum nonlocal equivalent strain leads to damage wrongly initiated away from the crack tip in a pre-cracked structure. Moreover, the standard nonlocal models fail to predict the shear damage band.

Inability to accurately model marco cracks: The issue of spurious damage distribution, also known as damage widening, characterized by an overly diffused appearance at the end stage of fracture (Giry et al., 2011). This phenomenon, referred to as a spreading effect by Zhang et al. (2021), indicates the unrealistic computational predictions on visible marco cracks. An analytical analysis of gradient-enhanced models in Xue et al. (2024) demonstrates that this damage widening is an inherent defect of these models.

Excessive spurious energy dissipation: The commonly employed scaling procedure near the physical boundary of the body may result in excessive energy dissipation (Jirásek et al., 2004), leading to an undesirable boundary effect. This effect results in nonlocal variables near the boundary being lower than those in the central regions. Consequently, this leads to damage erroneously spreading toward the region below the notch tip (Zhang et al., 2021).

Incorrect shielding effects: Nonlocal interactions are often inaccurately transmitted across notches or cracks, suggesting that material points on opposing sides of a crack should not influence each other despite their small geometric distance (Peerlings et al., 2001). Both existing non-convex boundaries and newly forming crack surfaces can induce shielding effects and thus require careful consideration in computational analyses (Pijaudier-Cabot and Dufour, 2010).

Size effect phenomenon reproduction: The standard nonlocal models are unable to replicate the size effect phenomenon in notched and unnotched concrete beams using a unified set of material parameters (Grégoire et al., 2013; Havlásek et al., 2016; Jirásek et al., 2004).

These nonphysical behaviors observed in conventional nonlocal damage models are largely attributed to the use of a constant interaction domain. Fundamental research by Bažant (1991, 1994) has identified micro-crack interactions as the primary source of nonlocality. Since the size of the active FPZ, where micro-cracking occurs, is evolving with the fracturing process (Jankowski and Styś, 1990; Otsuka and Date, 2000), the nonlocal interaction of the nonlocal continuum models should not remain constant. To address this limitation, the concept of an evolving characteristic length was introduced in Geers et al. (1998). This section presents modifications to standard nonlocal damage models that incorporate varying physical principles, aiming to improve the modeling of the gradual shift from diffuse damage to strain localization. These adaptations seek to bridge CDM and fracture mechanics, enabling a more accurate representation of quasi-brittle material behavior.

Strain-based transient formulation

To mitigate the issue of energy transfer from highly damaged zones to neighboring elastically unloading regions, the introduction of a new nonlocal damage model incorporating a transient internal length was proposed in Geers et al. (1998). This model introduces a novel variable, gradient activity $\zeta$ , which replaces the constant $c$ in equation (18). Gradient activity is designed to selectively activate nonlocal interactions within areas of significant damage. It is directly related to the local equivalent strain ( ${\epsilon }_{\mathrm{eq}}$ ) and can be expressed as follows:

$$\zeta \left({\epsilon }_{\mathrm{eq}}\right)=\{\begin{array}{cc}c{\left(\frac{{\epsilon }_{\mathrm{eq}}}{{\epsilon }_{\mathrm{eq},\zeta }}\right)}^{n_{\zeta }}& \text{if }{\epsilon }_{\mathrm{eq}}\le {\epsilon }_{\mathrm{eq},\zeta }\\ c& \text{if }{\epsilon }_{\mathrm{eq}}>{\epsilon }_{\mathrm{eq},\zeta }\end{array}$$

(29)

where ${\epsilon }_{\mathrm{eq},\zeta }$ and $n_{\zeta }$ are model parameters. The gradient activity $\zeta$ increases with the equivalent strain state, yet it is constrained by the parameter $c$ when the local equivalent strain reaches ${\epsilon }_{\mathrm{eq},\zeta }$ .

This strain-based transient formulation effectively maintains local material behavior around highly damaged areas, preventing the expansion of the damage zone. However, this approach requires additional continuity for the gradient activity, introducing an extra degree of freedom per node and thereby increasing computational costs (Geers et al., 1998). To address these challenges, a simplified version of the gradient model was subsequently proposed in Saroukhani et al. (2013):

$$\frac{{\overline{\epsilon }}_{\mathrm{eq}}}{\zeta }-{\mathrm{\nabla }}^2{\overline{\epsilon }}_{\mathrm{eq}}=\frac{{\epsilon }_{\mathrm{eq}}}{\zeta }$$

(30)

with

$$\zeta \left({\epsilon }_{\mathrm{eq}}\right)=\{\begin{array}{cc}{c}_0+(c-c_0){\left(\frac{{\epsilon }_{\mathrm{eq}}}{{\epsilon }_{\mathrm{eq},\zeta }}\right)}^{n_{\zeta }}& \text{if }{\epsilon }_{\mathrm{eq}}\le {\epsilon }_{\mathrm{eq},\zeta }\\ c& \text{if }{\epsilon }_{\mathrm{eq}}>{\epsilon }_{\mathrm{eq},\zeta }\end{array}$$

(31)

where ${\overline{\epsilon }}_{\mathrm{eq}}$ denotes the nonlocal equivalent strain; $c_0$ is assigned a positive value that is less than the square of the shortest distance between any two integration points. This setting ensures $\zeta \ne 0$ in equation (30) and maintains the equivalence between equations (31) and (29).

Pijaudier-Cabot et al. (2004) found that nonlocal effects intensify as damage progresses due to the porosity effect through the superposition theorem (Kachanov, 1987). Based on this understanding, Pijaudier-Cabot et al. (2004) proposed the application of a strain-based evolving internal length, $l_c({\epsilon }_{\mathrm{eq}})$ , to integral nonlocal models. This approach utilizes an evolving internal length in replace of $l$ in equation (9) to enhance the capacity of integral nonlocal models:

$$l_{\mathrm{c}}({\epsilon }_{\mathrm{eq}})={\alpha }_{0}f({\epsilon }_{\mathrm{eq}})+l_{\mathrm{c}0}$$

(32)

where ${\alpha }_0$ serves as a material parameter that can be set equal to the maximum internal length $l_c$ . $f({\epsilon }_{\mathrm{eq}})$ aligns with the damage evolution function which provides a range from 0 to 1. $l_{\mathrm{c}0}$ is related to the initial defects present in the material.

The above description of evolving internal length is based on the idea that the length scale parameter expands with increasing damage, eventually stabilizing at a predetermined internal length ( $l_{\mathrm{c}}$ ). It might seem counterintuitive at first since the increased length scale parameter indicates the nonlocal interaction domains expand. This mechanism plays an important role in the regularization process. Initially, the progressive nonlocal interaction from a minimal value can effectively prevent diffuse damage distribution. As fracture progresses, these interactions are appropriately confined, effectively avoiding spurious damage growth at higher deformation levels. This is achieved because certain regions surrounding the process zone are unloaded. In these unloaded areas, the gradient activity diminishes and material response is local. Consequently, the interactions between the damaged and unloaded zone disappear (Geers et al., 1998; Pijaudier-Cabot et al., 2004; Saroukhani et al., 2013).

Damage-dependent method

Geers et al. (2000) proposed several purely phenomenological nonlocal gradient models, including one that incorporates a damage-based transient internal length. According to Geers et al. (2000), the gradient parameter $c$ in equation (18) is dependent on the local damage variable $d$ , which can be simply expressed as follows:

$$c(d)=(1-d)c_{max}$$

(33)

in which the intrinsic length diminishes linearly with the increase of damage, approaching zero as the local damage variable nears unity. This reduction effectively eliminates the nonlocal treatment, ensuring the problem remains well-posed even after full crack initiation. However, this simple transient model cannot provide numerically stable and physically meaningful solutions for analyzing crack propagation phenomena (Geers et al., 2000).

Based on the model described by equation (33), Tran et al. (2023) and Tran and Bui (2023) introduced an enhanced gradient-damage model for the analysis of brittle fracture through the application of an energy limiter method:

$$\frac{\overline{d}}{(1-d)+\zeta }-c{\mathrm{\nabla }}^2\overline{d}=\frac{d}{(1-d)+\zeta }$$

(34)

where $\zeta$ , suggested to be set at ${10}^{-10}$ , is a small positive number for numerical stability purpose. $\overline{d}$ is the primary unknown nonlocal damage variable. The internal length refers to a problem-specific parameter associated with the energy limiter and the fracture energy of the material.

Nonlocality should not only diminish in the direction normal to the boundary at the boundary itself but also account for the shielding effect, which assumes that points separated by a crack should not interact. Pijaudier-Cabot and Dufour (2010) further addressed this by adjusting the nonlocal interaction distance between points $\mathit{x}$ and $\mathit{\xi }$ . This adjustment involves scaling the interaction distance by a function, $\gamma (\mathit{\xi })$ , a strategy similar to that employed in Krayani et al. (2009). According to Pijaudier-Cabot and Dufour (2010), the function $\gamma (\mathit{\xi })$ is defined as:

$$\gamma (\mathit{\xi })=\sqrt{1-d(\mathit{\xi })}$$

(35)

where $d(\mathit{\xi })$ is the local value of damage computed at point $\mathit{\xi }$ . As the damage increases, $\gamma (\mathit{\xi })$ approaches zero, which effectively results in the interaction distance becoming infinitely large. Upon the formation of a macro-crack, new boundary surfaces can be implicitly considered by assigning a damage value of 1.0.

By considering the evolving length effect proposed by Pijaudier-Cabot et al. (2004) and the boundary effect discussed by Krayani et al. (2009), Chen et al. (2024) introduced a modified characteristic length scale that incorporates two exponential terms. This length scale is expressed as:

$$l_{\mathrm{c}}=\left[\frac{e^{-max(d(\mathit{x}),d(\mathit{\xi }))}-e^{-1}}{1-e^{-1}}\right][1-e^{-\frac{D(\mathit{x})}{l_{c0}}}]l_{\mathrm{c}0}$$

(36)

where $D(\mathit{x})$ represents the distance of point $\mathit{x}$ from the boundary. In equation (36), the first term adjusts the length scale based on the damage variable, ensuring that $l_{\mathrm{c}}$ prior to damage initiation and increases gradually with the evolution of damage. The second term reduces $l_{\mathrm{c}}$ near the boundary to incorporate the boundary effect, aligning with the concept of the distance-based method, which will be discussed in the following section.

Distance-based method

The concept of distance-based averaging involves adjusting the averaging function at a point according to its distance from the nearest boundary, thereby mitigating spurious energy dissipation. This approach was first proposed by Bolander Jr and Hikosaka (1995), where its application was demonstrated by reducing the characteristic length near the notch tip.

To account for the boundary effect, where interaction stresses are expected to vanish at the boundary of a solid, Krayani et al. (2009) introduced a method that involves remapping the coordinate system utilized in the weighting function during the computation of averages near the boundary. The transformation of the coordinate system in 3D is defined as follows:

$$\Vert \mathit{x}-\mathit{\xi }\Vert =l_{\mathrm{c}}\sqrt{\frac{(x_1-{\xi }_1{)}^2}{a^2}+\frac{(x_2-{\xi }_2{)}^2}{b^2}+\frac{(x_3-{\xi }_3{)}^2}{c^2}}$$

(37)

where $a$ , $b$ , and $c$ represent the minimum values between the internal length and the distance from a given point to the nearest boundary in three directions in a 3D case. As the values of $a$ , $b$ , and $c$ decrease, the effective interaction distance between points $\mathit{x}$ and $\mathit{\xi }$ increases, mitigating spurious nonlocal interactions near the boundary. Moreover, the standard nonlocal formulation can be restored when $a$ , $b$ , and $c$ all equal to the internal length $l_{\mathrm{c}}$ .

As reported by Krayani et al. (2009), the scaling of the interaction domain results in it becoming elliptical and eventually degenerating into a straight segment. Similarly, Grassl et al. (2014) modified the averaging function by associating the basic weight function, ${\alpha }_0$ , with the minimum distance from point $\mathit{x}$ to the specimen boundary. This modification is equivalent to altering the internal length:

$$l_{\mathrm{c}}(\mathit{x},\mathit{\xi })=l_{\mathrm{c}}\cdot \gamma (\mathit{x})$$

(38)

with

$$\gamma (\mathit{x})=\{\begin{array}{cc}1,& D(\mathit{x})⩾tR\\ \frac{1-\beta }{tR}D(\mathit{x})+\beta ,& D(\mathit{x})<tR\end{array}$$

(39)

where $\beta$ and $t$ are model parameters. $D(\mathit{x})$ represents the minimum distance of point $\mathit{x}$ to the specimen boundary. The modification described in equation (38) addresses the issue by maintaining a reduced but circular influence domain, in contrast to the elliptical domain derived from equation (37). Moreover, this approach ensures that material behavior becomes fully local when the point is situated on the boundary. To achieve a more uniform distribution of dissipated energy near the notch tip, Grassl et al. (2014) also proposed a smoothed exponential formulation:

$$\gamma (\mathit{x})=1-(1-\beta )\exp (-\frac{D(\mathit{x})}{tR})$$

(40)

The distance-based method has been shown to accurately model energy dissipation near notches (Grassl et al., 2014). Furthermore, this approach effectively predicts the experimental size effect in both unnotched and notched beams using a single set of input parameters (Havlásek et al., 2016). However, the implementation of this method involves determining two additional input parameters, as outlined in equations (39) and (40).

Local-complement method

Methods that preserve the symmetry of nonlocal weight functions near the boundary, such as that outlined in Borino et al. (2003), are referred to as local-complement approaches. This term stems from these methods incorporating additional terms into the standard nonlocal formulation. The key idea is to compensate for the contribution of material lost beyond the physical boundaries by multiplying the local value at the receiving point $\mathit{x}$ by an appropriate factor. For instance, Grassl et al. (2014) introduced a generalized weight distribution formula:

$$\alpha (\mathit{x},\mathit{\xi })={\alpha }_{\mathrm{\infty }}(\Vert \mathit{x}-\mathit{\xi }\Vert )+v(\mathit{x})\delta (\mathit{x}-\mathit{\xi })$$

(41)

with

$$v(\mathit{x})=1-{\int }_V{\alpha }_{\mathrm{\infty }}(\Vert \mathit{x}-\mathit{\xi }\Vert )\mathrm{d}\mathit{\xi }$$

(42)

where ${\alpha }_{\mathrm{\infty }}$ represents the normalized weight function unaffected by boundary influences. $v(\mathit{x})$ quantifies the relative weight of the volume absent which diminishes to zero for points distanced from the boundary seen from Figure 7(a). $\delta$ is the Dirac distribution. Figure 7(a) reports two distinct contributions from nonlocal and local terms in a one-dimensional (1D) bar with fixed internal length. Following this modification, the computation of the nonlocal variable is defined as:

$$\overline{f}(\mathit{x})={\int }_V{\alpha }_{\mathrm{\infty }}(\mathit{x},\mathit{\xi })f(\mathit{\xi })\mathrm{d}\mathit{\xi }+v(\mathit{x})f(\mathit{\xi })$$

(43)

in which the first term on the right-hand side represents the standard nonlocal variable, which remains unaffected by boundary influences. The second term introduces an additional component that incorporates the local value at the point of interest.

Figure 7.

Diagram showing the contribution of different terms of (a) equation (41) and (b) equation (44) for one-dimensional (1D) bar in a uniform state of damage. Referred to Borino et al. (2003) and Krayani et al. (2009).

To mitigate nonlocal effects at the boundary of a solid, a strategic modification of the weight function is required in proximity to the boundary. According to Krayani et al. (2009), this involves conceiving an extension of the solid beyond its physical boundary by introducing a fictitious solid. This fictitious extension is designed to have a state of strain that is symmetric with respect to the actual boundary. Consequently, the computation of the nonlocal variable for a point $\mathit{x}$ near the boundary incorporates the standard expression, adjusted by subtracting the contribution from its symmetric counterpart, ${\mathit{x}}^{\star }$ , located within the fictitious solid. The formulation of this revised weight function is specified as follows:

$$\begin{array}{rl}{\alpha }^{\star }(\mathit{x},\mathit{\xi })& ={\alpha }_{\mathrm{\infty }}(\Vert \mathit{x}-\mathit{\xi }\Vert )-{\alpha }_{\mathrm{\infty }}(\Vert {\mathit{x}}^{\star }-\mathit{\xi }\Vert )\\ & +v^{\star }(\mathit{x})\delta (\mathit{x}-\mathit{\xi })\end{array}$$

(44)

with

$$v^{\star }(\mathit{x})=1-{\int }_V{\alpha }_{\mathrm{\infty }}(\Vert \mathit{x}-\mathit{\xi }\Vert )\mathrm{d}\mathit{\xi }+{\int }_V{\alpha }_{\mathrm{\infty }}(\Vert {\mathit{x}}^{\star }-\mathit{\xi }\Vert )\mathrm{d}\mathit{\xi }$$

(45)

Figure 7(b) shows contributions from different terms in equation (44) in a 1D bar with uniform damage. For points located far away from the boundary, the contribution from the symmetric point ${\mathit{x}}^{\star }$ vanishes, allowing for the recovery of the standard integral nonlocal damage model. Conversely, at the extremities of the bar, the material response becomes inherently local, as indicated by equation (44) when $\mathit{x}={\mathit{x}}^{\star }$ . This aligns with the intuitive understanding that nonlocal effects should diminish at the boundary of a solid. However, the applicability of this concept is restricted to 1D implementations. The challenge arises in defining the symmetric counterpart for any point within 2D or 3D frameworks, as a given material point may have multiple overlapping symmetric counterparts (Krayani et al., 2009).

Interaction-based method

A novel interaction-based nonlocal formulation, explicitly computing nonlocal interactions through the dilation of circular inclusions centered on each material point, was introduced in Rojas-Solano et al. (2013) and Pijaudier-Cabot and Grégoire (2014). Within the integral formulation framework, the basic weight function in equation (7) is defined as follows:

$${\alpha }_0(\mathit{x},\mathit{\xi })\equiv A^{*}(\mathit{x},\mathit{\xi },a)=\frac{A(\mathit{x},\mathit{\xi },{\mathit{\epsilon }}^{*},a)}{\Vert {\mathit{\epsilon }}^{*}\Vert }$$

(46)

where $A^{*}(\mathit{x},\mathit{\xi },a)$ is the interaction function that is governed by the positions of two points and the radius $a$ of a circular inclusion centered at point $\mathit{\xi }$ ; $A(\mathit{x},\mathit{\xi },{\mathit{\epsilon }}^{*},a)=\sqrt{\sum _{i=1}^3|{\epsilon }_i(\mathit{x}){|}^2}$ is the norm of strain tensor at point $\mathit{x}$ with the dilation ${\mathit{\epsilon }}^{*}$ from point $\mathit{\xi }$ .

To address both boundary and shielding effects, it was proposed that the size of the inclusions should decrease as a function of their distance to the boundary and in relation to the damage level. This approach was suggested by Rojas-Solano et al. (2013) as follows:

$$a(\mathit{x})=\text{min}(a_0\sqrt{1-d(\mathit{x})},D(\mathit{x}))$$

(47)

where $a_0$ represents a model parameter associated with the maximum aggregate size. $d(\mathit{x})$ quantifies the local damage at point $\mathit{x}$ , $D(\mathit{x})$ denotes the minimal distance from point $\mathit{x}$ to the nearest free boundary.

Stress-based formulation

Bažant (1994) introduced an integral nonlocal model characterized by an anisotropic interaction kernel with a constant length scale, where the principal stress directions determine the shape of the nonlocal averaging domain. This concept was further advanced by Giry et al. (2011), where a transient length scale is incorporated to enable the nonlocal interaction to diminish as the principal stress perpendicular to the crack decreases. According to Giry et al. (2011), the nonlocality is based on the source point $\mathit{\xi }$ rather than the receiver point $\mathit{x}$ . The new interaction domain around $\mathit{\xi }$ is established as an ellipsoid, with its principal axes determined by the stress direction at the source point and magnitude specified as follows:

$$a_i={\rho }_i{l}_{\mathrm{c}}$$

(48)

with

$${\rho }_i=\{\begin{array}{cc}\left|\frac{{\sigma }_i}{f_{\mathrm{t}}}\right|& \text{if}-f_{\mathrm{t}}<{\sigma }_i<f_{\mathrm{t}}\\ 1& \text{if}\text{otherwise}\end{array}$$

(49)

where ${\rho }_i$ represents the parameters that characterize the stress state. ${\sigma }_i$ denotes the principal stresses. $f_{\mathrm{t}}$ is the material tensile strength.

The internal length $l_{x\xi }$ between the source point $\mathit{\xi }$ and the receiver point $\mathit{x}$ is defined based on the geometric configuration of an ellipsoid centered at $\mathit{\xi }$ . Specifically, $l_{x\xi }$ corresponds to the radial distance from $\mathit{\xi }$ to $\mathit{x}$ within this ellipsoidal domain, following the direction vector $(\mathit{\xi }-\mathit{x})$ . The expression for calculating $l_{x\xi }$ is given as:

$$l_{x\xi }={\rho }_{x\xi }{l}_{\mathrm{c}}$$

(50)

with

$${\rho }_{x\xi }^2=\frac{1}{[{\sin }^2(\varphi ){\cos }^2(\theta )]/{\rho }_1^2+[{\sin }^2(\varphi ){\cos }^2(\theta )]/{\rho }_2^2+({\cos }^2(\varphi )/{\rho }_3^2)}$$

(51)

where ${\rho }_{x\xi }$ represents the scaled factor, which is derived within the framework of spherical coordinates denoted as $(\rho ,\varphi ,\theta )$ . Consequently, the value of ${\rho }_{x\xi }$ varies from 0, indicating principal stresses tend to 0, to 1, demonstrating principal stresses achieve the maximum.

The intensity and shape of the nonlocal interaction domain are governed by the stress state, allowing nonlocal interactions to be automatically canceled in regions where principal stress diminishes. Consequently, the influence of the boundary on a specific material point is implicitly perceived by the stress state of the source point. An illustration of this principle is provided in Figure 8(a), which depicts the elliptical influence exerted by the point located at $\mathit{\xi }$ on its surroundings. As shown in Figure 8(b), point C is unaffected by the notch, displaying a circular influence. Point B, positioned near the crack tip, experiences a highly oriented stress state, leading to an elliptical influence domain. Point A, shielded by the crack, is subjected to a low-stress state and exerts minimal influence on adjacent points. This approach has been validated for its effectiveness in modeling damage initiation and propagation without explicit boundary treatments. It is extended to integration with rate-dependent damage models, as demonstrated in Pereira et al. (2016).

Figure 8.

Stress-based nonlocal model proposed in Giry et al. (2011): (a) influence of a point and (b) isovalues of the influence of various points in the specimen.

While the stress-based model effectively corrects damage initiation and reduces artificial damage spreading, it exhibits some mesh dependency due to the length scale parameters diminishing as principal stresses approach zero. To address the challenge posed by zero stress in certain directions, a minimum value for $l_{x\xi }$ is set equal to the characteristic size of the element in Giry et al. (2011). Grassl et al. (2014) proposed a similar strategy in which the modification of the nonlocal weight function depends on the effective stress state at the “receiver” point $\mathit{x}$ .

The principal stress-based model was further developed into a gradient-enhanced variant by Bongers (2011), where the nonlocal interaction domain at the receiver point $\mathit{x}$ is adjusted by stretching it along the principal axes. In addressing a 3D problem, the constant gradient parameter $c$ in equation (18) is replaced by a set of stress-based evolving parameters. Consequently, the standard gradient model is reformulated as:

$$\overline{f}(\mathit{x})-c_1\frac{{\mathrm{\partial }}^2\overline{f}(\mathit{x})}{\mathrm{\partial }{x_1}^2}-c_2\frac{{\mathrm{\partial }}^2\overline{f}(\mathit{x})}{\mathrm{\partial }{x_2}^2}-c_3\frac{{\mathrm{\partial }}^2\overline{f}(\mathit{x})}{\mathrm{\partial }{x_3}^2}=f(\mathit{x})$$

(52)

in which $[x_1,x_2,x_3]$ denote the principal coordinates. The gradient parameters $c_1$ , $c_2$ , and $c_3$ are determined by taking the dependence of the characteristic length on the stress state into account, with their definitions inspired by equation (49):

$$c_{1,2,3}=\frac{{\sigma }_{1,2,3}^2}{f_t^2}c$$

(53)

Furthermore, Thai et al. (2016) developed a higher-order gradient-enhanced model within the stress-based framework of Bongers (2011). This model was specifically designed to mitigate the truncation effects associated with higher-order gradient terms, thereby improving the accuracy and stability of the formulation.

Nguyen et al. (2018) highlighted that equation (53) does not effectively prevent spurious damage growth due to the incorrect calculation of the nonlocal equivalent strain. To address this issue, Nguyen et al. (2018) proposed a modified diffusive interaction domain incorporating the local equivalent strain, formulated as follows:

$$c_{1,2,3}={\left(\frac{{\kappa }_0}{max({\kappa }_0,{\epsilon }_{\mathrm{eq}})}\right)}^{\eta }\frac{{\overline{\sigma }}_{1,2,3}^2}{f_{\mathrm{t}}^2}c$$

(54)

where ${\kappa }_0$ represents the damage threshold. ${\overline{\sigma }}_{1,2,3}^2$ denote the smoothed principal stresses proposed in Nguyen et al. (2018) to mitigate stress oscillation effects. $\eta$ is defined as a coefficient potentially related to the rate of decrease in nonlocal interaction. As the equivalent strain ${\epsilon }_{\mathrm{eq}}$ increases, there is a corresponding reduction in the anisotropic gradient coefficients $c_{1,2,3}$ , which in turn leads to a diminished nonlocal interaction. This concept aligns with observations by Poh and Sun (2017), which assert that nonlocal interactions decrease with deformation after damage initiation.

Moreover, Nguyen et al. (2018) reformulated the implicit gradient-enhanced model described in equation (18) in terms of the principal coordinates $(x_1,x_2,x_3)$ , which is given as

$$\overline{f}(\mathit{x})-\mathrm{\nabla }\cdot (\mathbf{c}\mathrm{\nabla }\overline{f}(\mathit{x}))=f(\mathit{x})$$

(55)

with

$$\mathbf{c}=\left[\begin{array}{ccc}{c}_1& 0& 0\\ 0& c_2& 0\\ 0& 0& c_3\end{array}\right]$$

(56)

where the gradient parameters $c_{1,2,3}$ are determined according to equation (54).

Meanwhile, Vandoren and Simone (2018) introduced a general format within a Cartesian coordinate system, where the anisotropic gradient matrix $\mathbf{c}$ is expressed as follows:

$$\mathbf{c}={\mathbf{c}}_{\sigma }=\mathbf{R}\left[\begin{array}{ccc}{c}_1& 0& 0\\ 0& c_2& 0\\ 0& 0& c_3\end{array}\right]{\mathbf{R}}^{⊺}$$

(57)

where $\mathbf{R}$ represents the rotation matrix that realizes the transformation between the principal axis system, $(x_1,x_2,x_3)$ , and the Cartesian coordinate system, $(x,y,z)$ . The gradient parameters $c_{1,2,3}$ are derived based on the principal stresses using equation (53).

The anisotropic gradient presented in equation (57) does not differentiate between compressive and tensile stresses. To address this limitation, Vandoren and Simone (2018) also proposed a modified formulation:

$$\mathbf{c}={\mathbf{c}}_{\overline{\sigma }}=\frac{{\overline{\sigma }}_{\mathrm{eq}}{f}_{\mathrm{t}}}{max({\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2)}{\mathbf{c}}_{\sigma }$$

(58)

where ${\overline{\sigma }}_{\mathrm{eq}}$ is the nonlocal principal stress and calculated by ${\overline{\sigma }}_{\mathrm{eq}}=(1-w)E{\overline{\epsilon }}_{\mathrm{eq}}$ ; $w$ is the damage evolution.

The modified interaction domain defined by equation (58) expands with increasing deformation in intact material. However, it diminishes with deformation after damage initiation, an observation that aligns with the description provided in equation (54). A small modification was made to equation (58) by Lale et al. (2023) to accelerate the decrease of nonlocal interactions in the direction perpendicular to the maximum stress.

The stress-based approach, while not explicitly accounting for boundary influences, effectively “feels” them through variations in the stress field. In such models, an anisotropic interaction kernel is utilized to control both the intensity and the geometry of the nonlocal interaction domain. Given that this kernel functions as a stress-dependent tensor, nonlocal interactions are automatically nullified in areas experiencing a decrease in stress values.

Localizing gradient damage model

Based on the generalized micromorphic theory as detailed by Forest (2009), a thermodynamically consistent localizing gradient damage model was initially introduced by Desmorat and Gatuingt (2007b) and Poh and Sun (2017). This model incorporates a micro-continuum at each macroscopic point within the domain to characterize the microscopic deformation. According to Poh and Sun (2017), the free energy density function within this framework is expressed as follows:

$$\psi =\frac{1}{2}(1-\omega )\mathit{\epsilon }:\mathbb{C}:\mathit{\epsilon }+\frac{1}{2}h(e-\tilde{e}{)}^2+\frac{1}{2}gh{l}^2\mathrm{\nabla }\tilde{e}\cdot \mathrm{\nabla }\tilde{e}$$

(59)

In equation (59), the first term corresponds to the conventional continuum damage model, in which $\omega$ is the damage variable. The second term is responsible for micro–macro interaction with $h$ as the coupling modulus. Here, $e$ and $\tilde{e}$ denote the equivalent strains at macroscopic and microscopic scales, respectively. The last term addresses linearized interactions among micro-cracks via $\mathrm{\nabla }\tilde{e}$ . Notably, Poh and Sun (2017) proposed that the interaction function $g$ in equation (59) diminishes as damage progresses, formulated as:

$$g=\frac{(1-R)\exp (-\eta \omega )+R-\exp (-\eta )}{1-\exp (-\eta )}$$

(60)

where $R$ and $\eta$ serve as constants that govern the degradation of the interaction domain as the damage evolves at a material point. The impact of these constants on numerical performance was extensively analyzed by Sarkar et al. (2019). According to equation (60), the interaction domain diminishes with increasing damage within the free energy statement, thereby mitigating spurious energetic interactions across cracks.

The governing equations of the localizing gradient model consist of the following equilibrium equation (without body force) and microforce balance equation:

$$\mathrm{\nabla }\cdot \mathit{\sigma }=\mathbf{0}$$

(61a)

$$\tilde{\mathit{\sigma }}=\mathrm{\nabla }\cdot \tilde{\mathit{\xi }}$$

(61b)

where $\tilde{\mathit{\sigma }}$ denotes the coupling stress tensor and $\tilde{\mathit{\xi }}$ is the moment stress tensor. Employing the Coleman–Noll procedure allows for the derivation of constitutive relations as follows:

$$\mathit{\sigma }=\frac{\mathrm{\partial }\psi }{\mathrm{\partial }\mathit{\epsilon }}=(1-\omega )\mathbb{C}:\mathit{\epsilon }+h(e-\tilde{e})\frac{\mathrm{\partial }e}{\mathrm{\partial }\mathit{\epsilon }}$$

(62a)

$$\tilde{\mathit{\sigma }}=\frac{\mathrm{\partial }\psi }{\mathrm{\partial }\tilde{e}}=h(\tilde{e}-e)$$

(62b)

$$\tilde{\mathit{\xi }}=\frac{\mathrm{\partial }\psi }{\mathrm{\partial }\mathrm{\nabla }\tilde{e}}=gh{l}^2\mathrm{\nabla }\tilde{e}$$

(62c)

where equation (62a) represents the classical stress tensor, including the coupling effects between the micro and macro levels. Equation (62b) signifies the equivalent coupling stress arising from interactions between micro- and macro-cracks. Equation (62c) describes the moment stress resulting from interactions among micro-cracks. By integrating equations (62b) and (62c) into the microforce balance equation presented in equation (61b), the localizing gradient model are established as follows:

$$\tilde{e}-e=\mathrm{\nabla }\cdot (g{l}^2\mathrm{\nabla }\tilde{e})$$

(63)

In the framework of the localizing gradient model, a “pseudo-crack” can be described as the damage variable approaches unity, indicating a natural damage-to-fracture transition. This feature enables the model to effectively circumvent the formation of spurious damage zones during the final stages of the softening process (Huang et al., 2022; Sarkar et al., 2020; Wang et al., 2022; Zhang et al., 2021). Additionally, its adaptability has been validated through implementation in commercial FE software (Sarkar et al., 2022, 2019; Zhang et al., 2022). For instance, user subroutines for implementing both conventional and localizing gradient-enhanced models in ABAQUS are made available in Sarkar et al. (2019) and Zhang et al. (2022), while supplemental material for gradient model implementation in COMSOL is provided in Sarkar et al. (2022). Beyond capturing tensile fracture phenomena, the model also accurately reflects compressive behavior (Wang et al., 2023) and mixed failure modes (Shedbale et al., 2021) with designed techniques.

However, Negi and Kumar (2019) observed a minor deviation in the maximum damage localization from the crack tip with the original localizing model introduced by Poh and Sun (2017), which utilized an assumption of an isotropic interaction function domain. To address this, an anisotropic gradient tensor, inspired by the stress-based approach in Vandoren and Simone (2018), was incorporated into equation (63) by Negi and Kumar (2019). This tensor adjusts the orientation and size of the interaction domain based on the principal stress direction and the damage state at specific material points, ensuring that maximum damage values occur directly in front of the crack tip during the initial loading stages. Furthermore, Negi et al. (2020) refined the model by substituting the micromorphic equivalent strain with the micromorphic strain tensor to more accurately represent nonlocal interactions, owing to its immunity to oscillation. This enhancement avoids the use of higher-order elements or separate smoothing procedures previously adopted by Negi and Kumar (2019). Moreover, Negi et al. (2020) highlighted that the generalized anisotropic gradient-enhanced model described in equation (55) can be derived from the micromorphic framework, with equation (63) representing a special case of equation (55) under the assumption of an isotropic interaction domain. More applications of the modified localizing gradient model can be found in Negi et al. (2021, 2022).

Eikonal nonlocal formulation

To account for the influence of cracks and damaged zones between points, Desmorat and Gatuingt (2007a, 2007b) introduced a new nonlocal formulation that retains the integral framework of nonlocal theories. Instead of utilizing the classical distance metric $r=\Vert \mathit{x}-\mathit{\xi }\Vert$ in weight functions between point $\mathit{x}$ and $\mathit{\xi }$ , this model adopts the information propagation time ${\tau }_{x\xi }$ as a novel measure. Consequently, the integral formulation of this approach is expressed as follows:

$$\begin{array}{rl}& \overline{f}(\mathit{x})=\frac{1}{V_0(\mathit{x})}{\int }_V{\alpha }_0\left(\frac{{\tau }_{x\xi }}{{\tau }_{\mathrm{c}}}\right)f(\mathit{\xi })\mathrm{d}\mathit{\xi }\\ & \text{with}{V}_0(\mathit{x})={\int }_V{\alpha }_0\left(\frac{{\tau }_{x\xi }}{{\tau }_{\mathrm{c}}}\right)\mathrm{d}\mathit{\xi }\end{array}$$

(64)

where the conventional internal length $l_{\mathrm{c}}$ in weight functions is substituted by a characteristic or internal time ${\tau }_c$ . This parameter ${\tau }_{\mathrm{c}}$ is defined as the time necessary for an elastic wave to cover a distance $r$ within intact material. Correspondingly, ${\tau }_{x\xi }$ denotes the actual time required for wave propagation from point $\mathit{x}$ to $\mathit{\xi }$ , a metric that naturally increases with the progression of material damage.

The principle of this methodology is that the wave celerity, ${\tilde{c}}_0$ , depending on the stiffness of the material, enables the adjustment of nonlocal interactions as damage progresses. In a homogeneous medium with a scalar damage $d$ , the wave celerity is described as follows:

$${\tilde{c}}_0=\sqrt{\frac{(1-d)E}{\rho }}=c_0\sqrt{1-d}$$

(65)

where $\rho$ represents the density of the material and $c_0$ denotes the wave celerity without damage. Consequently, the internal time required for wave propagation between two points is formulated as:

$${\tau }_{x\xi }=\frac{\Vert \mathit{x}-\mathit{\xi }\Vert }{{\tilde{c}}_0}$$

(66)

By considering the relation $l_c=c_0{\tau }_{\mathrm{c}}$ between the internal length and internal time, an effective or “dynamic” distance $\tilde{r}$ between points $\mathit{x}$ and $\mathit{\xi }$ , which increases as damage progresses, thus is derived as follows:

$$\tilde{r}=c_0{\tau }_{x\xi }=\frac{\Vert \mathit{x}-\mathit{\xi }\Vert }{\sqrt{1-d}}$$

(67)

The “dynamic” distance between two points in equation (67) escalates toward infinity as damage approaches unity, resulting in the vanish of all interactions. This effectively addresses the issue of incorrect interactions across highly damaged zones in standard nonlocal approaches. Desmorat and Gatuingt (2007a, 2007b) illustrated that real cracks and highly damaged zones exhibit equivalence in 1D simulations within the internal time nonlocal framework. However, this model requires the computation of the propagation time ${\tau }_{x\xi }$ at each material point and for each time step, posing significant computational challenges for extending the approach to 2D and 3D simulations.

A novel eikonal nonlocal (ENL) formulation, inspired by the concept of internal time, was introduced by Desmorat et al. (2015). In this formulation, the traditional concept of information propagation time is innovatively replaced with the effective distance between two points. Desmorat et al. (2015) define this effective distance as the solution to an eikonal equation in the context of wave propagation with a Wentzel–Kramers–Brillouin approximation. The unified eikonal equation, accommodating both isotropic and anisotropic damage models, is formulated as follows:

$$\mathrm{\nabla }\tilde{r}\cdot {\mathbf{g}}^{-1}\cdot \mathrm{\nabla }\tilde{r}=1$$

(68)

where ${\mathbf{g}}^{-1}=\mathbf{1}-\mathbf{d}$ is defined as a damage-dependent tensor, alternatively referred to as the Riemannian metric.

Mathematically, the effective interaction distance between material points is determined by solving equation (68). The physical interpretation of the eikonal function is a continuous shortest-path problem. From the perspective of differential geometry, multiplying the information time by the wave velocity $c_0$ identifies the shortest path within a Riemannian space that is distorted by damage, as illustrated in Figure 9. Consequently, Desmorat et al. (2015) established a classical geodesic function via the Hamilton–Jacobi standard theory, which is expressed as follows:

$$\mathrm{d}{\tilde{r}}^2=\mathrm{d}\mathbf{x}\cdot \mathbf{g}\cdot \mathrm{d}\mathbf{x}$$

(69)

Figure 9.

The interaction distance between material points on a Riemannian space with a curvature caused by damage (Nogueira et al., 2024).

Both the eikonal equation (equation (68)) and the geodesic equation (equation (69)) determine the effective distance ${\tilde{r}}_{x\xi }$ between points $\mathit{x}$ and $\mathit{\xi }$ as the shortest distance in a curved spaced which should be a geodesic curve. As illustrated from the top left to the bottom right in Figure 9, the geodesic distance between two points, indicated by a red line, progressively extends toward infinity with increasing levels of damage.

With the effective distances determined, the resultant nonlocal integral formulation is expressed as follows:

$$\begin{array}{rl}& \overline{f}(\mathit{x})=\frac{1}{V_0(\mathit{x})}{\int }_V{\alpha }_0\left(\frac{\tilde{r}}{l_{\mathrm{c}}}\right)f(\mathit{\xi })\mathrm{d}\mathit{\xi }\\ & \text{with}{V}_0(\mathit{x})={\int }_V{\alpha }_0\left(\frac{\tilde{r}}{l_{\mathrm{c}}}\right)\mathrm{d}\mathit{\xi }\end{array}$$

(70)

Eikonal models are particularly effective for addressing newly created boundaries, as they inherently ensure that the effective interaction distance between highly damaged elements and their neighbors tends toward infinity, thereby stopping nonlocal interactions. This feature enables straightforward numerical integrations in 1D cases, as demonstrated in Desmorat et al. (2015), Jirásek and Desmorat (2019), and Thierry et al. (2020). Furthermore, Rastiello et al. (2018) applied the fast marching method to numerically solve the eikonal equation (equation (68)) to obtain effective interaction distances in 2D cases.

Moreover, Desmorat et al. (2015) introduced a gradient-enhanced version of the ENL model, which leads to the modified Helmholtz problem:

$$\tilde{f}(\mathit{x})-c\frac{1}{\sqrt{\text{det}\mathbf{g}}}\mathrm{\nabla }\cdot (\sqrt{\text{det}\mathbf{g}}{\mathbf{g}}^{-1}\cdot \mathrm{\nabla }\tilde{f}(\mathit{x})=f(\mathit{x})$$

(71)

The gradient version of the ENL model was first evaluated by Nogueira et al. (2022) in the context of a 1D spalling test. This study compared the performance of ENL models in both integral and gradient forms, reporting an unexpected diffusion of damage upon localization. This issue was addressed by limiting the damage to a critical value (Geers et al., 1998). Additional applications of the gradient ENL model can be found in Marconi (2022) and Nogueira et al. (2023). Furthermore, Nogueira et al. (2024) derived the eikonal implicit gradient formulation from thermodynamic principles based on a geometric extension of the micromorphic approach (Forest, 2009). This derivation revealed substantial similarities with the previously discussed localizing gradient damage model.

While the eikonal nonlocal damage model outperforms standard nonlocal damage models in simulating realistic crack paths, it also tends to produce load-displacement diagrams that feature pronounced snapbacks and overly brittle responses, diverging from experimental observations. This discrepancy stems from the increasing tendency of the model toward local behavior as damage levels rise. To mitigate this side effect, it is suggested to either explore alternative constitutive models, as recommended in Jirásek and Desmorat (2019) or to refine the calibration of model parameters, a strategy suggested in Wang et al. (2023).

Over-nonlocal formulation

The over-nonlocal model initially conceptualized for nonlocal plasticity by Strömberg and Ristinmaa (1996) and Vermeer and Brinkgreve (1994), was later extended to encompass nonlocal damage models integrated with plasticity by Grassl and Jirásek (2006). This method has proven to effectively mitigate spurious localization both in integral nonlocal models (Di Luzio and Bažant, 2005; Grassl and Jirásek, 2006; Mohamad-Hussein and Shao, 2007) and gradient models (Al-Rub and Voyiadjis, 2009; Poh and Swaddiwudhipong, 2009a, 2009b).

The fundamental principle of the over-nonlocal approach is adopting a weighted average of local and nonlocal variables to characterize the damage process. The formulation of the over-nonlocal approach is described as follows:

$$\hat{f}(\mathit{x})=m\overline{f}(\mathit{x})+(1-m)f(\mathit{x})$$

(72)

where $\hat{f}(\mathit{x})$ represents the effective variable. $m$ is a parameter that adjusts the weight between nonlocal and local components. By setting $m$ to 0, the model recovers to a purely local framework; conversely, setting $m$ to 1 retrieves the conventional nonlocal model.

The previous study in Lu et al. (2009), concerning simple softening plastic models, confirmed that the localization zone is finite only if $m>1$ , a strategy that underlines the concept referred to as “over-nonlocal.” Di Luzio (2007) suggests a range for $m$ of $1<m\le 1.1$ . Summersgill et al. (2017a) found that a high value of $\alpha$ and fine mesh can lead to instability in the over-nonlocal method. It is worth noting that equation (13), as proposed by Galavi and Schweiger (2010), was specifically designed to circumvent the “trial-and-error” approach typically required to determine an appropriate $m$ value in over-nonlocal models. This is because both equations (13) and (72) present a low weight at the center (Chen et al., 2023; Summersgill et al., 2017a).

In fact, equation (72) is similar to the formulation from a local-complement method, as both approaches simultaneously account for local and nonlocal variables. Bui (2010) and Nguyen (2011) have applied equation (72) to regularize pure damage models, which is termed as a mixed local and nonlocal formulation. Notably, the parameter $m$ evolves in response to the fracturing process, inspired by the concept of evolving interactions (Geers et al., 1998). Additionally, in Nguyen (2011), the evolving function $m$ , with a maximum value below 1, effectively captures a finite localization zone in Nguyen (2011).

Crack bandwidth estimation in fracture energy regularization

Determining the crack bandwidth or characteristic element length is essential in fracture energy regularization techniques to maintain the objectivity of energy dissipation. In principle, this measurement should reflect the evolving width of the FPZ throughout the fracture process (Comi and Perego, 2001). Methods that rely solely on the constant calculation of element area or volume hence fall short (Mosalam and Paulino, 1997). This recognition has led to a broader approach, applicable to elements of any shape, where the effective characteristic element length is measured in a direction normal to the local crack plane (Bažant and Planas, 2019; Comi and Perego, 2001). A straightforward estimation method involves aligning this measurement with the principal strain direction, which should be perpendicular to the crack band.

According to Jirásek and Bauer (2012) and Rots et al. (1988), the characteristic element length ( $h$ ) is determined by a number of factors, including the shape and size of the element, the interpolation function used, the integration scheme, and the orientation of the mesh line. This section summarizes the proposed methods for estimating $h$ , taking into account both the classical approaches reviewed by Jirásek and Bauer (2012) and more recent developments, particularly those concerning higher-order elements.

Oliver’s method

Oliver (1989) introduced two methods for estimating the crack bandwidth, specifically designed to address the challenges posed by irregular meshes. These methods have been extensively reviewed by Jirásek and Bauer (2012). The first approach employs an auxiliary function $\varphi =\sum _{i=1}^{n_c}{N}_i{\varphi }_i$ , which assigns values of 1 or 0 to corner nodes based on their position relative to the crack band, with the values being interpolated across the crack band using shape functions. The mathematical formulation of this method is expressed as follows:

$$h_{O1}(\mathbf{x})={(\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\mathbf{x}}\cdot \mathbf{n}(\mathbf{x}))}^{-1}={\left(\sum _{i=1}^{n_c}[\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }\mathbf{x}}\cdot \mathbf{n}(\mathbf{x})]{\varphi }_i\right)}^{-1}$$

(73)

where $n_c$ denotes the number of corner nodes, and $N_i$ is the standard shape function corresponding to corner node $i$ . The symbol $\mathbf{x}$ specifies the target point, such as a GP, where $\mathbf{n}(\mathbf{x})$ is defined as the unit vector oriented in the direction of the major principal strain at point $\mathbf{x}$ . For the purposes of this model, a crack is assumed to be perpendicular to $\mathbf{n}(\mathbf{x})$ . The parameter ${\varphi }_i$ serves as a crack indicator, which is set to 1 if the corner node $i$ is ahead of the crack front and 0 otherwise.

Oliver’s first method geometrically represents the distance measured between the boundary edges of an element, oriented in the normal direction to the crack. Importantly, this method computes the crack bandwidth individually for each Gauss point. To enhance this process, Jirásek and Bauer (2012) proposed a modification that a singular unit vector ${\mathbf{n}}_{\mathrm{c}}$ , representing the major principal strain direction at the center of the element, is employed instead of using the $\mathbf{n}(\mathbf{x})$ at each point $\mathbf{x}$ . This adjustment ensures a consistent crack bandwidth for all Gauss points within the same element.

The second method from Oliver (1989) was based on the concept of equivalent crack length. This length is defined as the distance between projections of the midpoints of the two-element edges intersected by the crack. Consequently, the crack bandwidth is determined by dividing the element area by the equivalent crack length. The expression for Oliver’s second method is expressed as follows:

$$h_{O2}=\frac{A_e}{{\mathit{n}}_{\mathrm{c}}\cdot {\int }_{{\mathrm{\Gamma }}_e}\varphi {\mathbf{n}}_e\mathrm{d}s}$$

(74)

where $A_e$ represents the element area, while ${\mathrm{\Gamma }}_e$ denotes the boundary of the element. ${\mathbf{n}}_e$ is defined as a unit vector perpendicular to the edge of the element and oriented outwards. Given that this method relies on the local information from the center point of an element, the calculated value of $h_{O2}$ is consistent for all GPs within the same element.

Jirásek and Bauer (2012) provided a thorough comparison of Oliver’s methods, highlighting their capacity to account for element shape, size, and crack direction. However, one challenge identified is that the estimated crack bandwidth may become infinite under certain circumstances (Jirásek and Bauer, 2012). Despite their comprehensive approach to incorporating key geometric and crack orientation factors, it is important to note that Oliver’s methods are limited to 2D problems.

Govindjee’s method

Govindjee et al. (1995) found that applying Oliver’s method to 3D problems with common brick elements leads to a discontinuous length function. This issue arises from the discontinuous nature of crack indicator ${\varphi }_i$ . To address this challenge, Govindjee et al. (1995) introduced a continuous definition for the crack indicator, enhancing its applicability to 3D problems:

$$h_G(\mathbf{x})={(\sum _{i=1}^{n_c}\left[\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }\mathbf{x}}{\tau }_i\right]\cdot \mathbf{n}(\mathbf{x}))}^{-1}$$

(75)

with

$$\begin{array}{rl}{\tau }_i& =\frac{({\mathbf{x}}_i-{\mathbf{x}}_c)\cdot \mathbf{n}(\mathbf{x})-{\tau }^{min}}{{\tau }^{max}-{\tau }^{min}}\\ {\tau }^{min}& =\underset{i=1,N_{\text{nodes}}}{min}\{({\mathbf{x}}_i-{\mathbf{x}}_c)\cdot \mathbf{n}(\mathbf{x})\}\\ {\tau }^{max}& =\underset{i=1,N_{\text{nodes}}}{max}\{({\mathbf{x}}_i-{\mathbf{x}}_c)\cdot \mathbf{n}(\mathbf{x})\}\end{array}$$

(76)

where ${\mathbf{x}}_i$ and ${\mathbf{x}}_c$ denote the corner nodal points and the central point of an FE, respectively. Similar to the modifications applied to Oliver’s first method in Jirásek and Bauer (2012), Govindjee’s method is applied using the central point of each element. This ensures that all Gauss points within a single element are assigned the same bandwidth.

Govindjee’s method geometrically determines the crack bandwidth as the maximum distance between the corner nodes, projected onto the direction normal to the crack. This approach accounts for the element shape, element dimensions, and crack direction. It aligns with the projection method introduced in Cervenka et al. (1995), which also expresses these geometric explanations. For addressing compression softening, a similar concept of a crushing band was proposed in Cervenka and Cervenka (2010) and Červenka et al. (2018), relying on the projection of each FE toward the direction of minimal compressive stress. To mitigate the bias introduced by mesh orientation, an orientation factor ranging from 1 to 1.5 was suggested in Cervenka et al. (1995) and Červenka et al. (2018). However, Jirásek and Bauer (2012) pointed out there is no need for this additional orientation factor.

Slobbe’s method

Oliver’s and Govindjee’s methods are only applicable to linear FEs and standard Gauss quadrature. Jirásek and Bauer (2012) identified a potential issue where strain softening might not fully localize across an entire element when quadratic interpolation functions are used in 1D elements, potentially leading to an overestimated crack bandwidth. Consequently, the use of higher-order elements was discouraged. This partial localization issue is attributed to the numerical integration scheme and interpolation function, a phenomenon that was also reported by Slobbe et al. (2013) in the context of 2D elements. To address these challenges in higher-order elements and incorporate mesh orientation considerations, Slobbe et al. (2013) introduced two parameters to refine Govindjee’s method:

$$h_S(\mathbf{x})=\alpha \cdot \gamma \cdot h_G(\mathbf{x})$$

(77)

where $\alpha$ acts as a factor to account for the effect of potential strain localization, with its values for various element types being specified in Table 1. These values are determined by calculating the ratio of softening GPs to the total number of GPs in a single element. $\gamma$ serves as an alignment factor, represented by a Fourier series function containing some undetermined coefficients. The calibration of these Fourier coefficients based on energy dissipation can be found in Slobbe et al. (2013).

Table 1.

$\alpha$ Values for different element type (Slobbe et al., 2013).

Types	Triangle	Quadrilateral	Quadrilateral
	full integration	reduced integration	full integration
Linear	1	1	–
Quadratic	2/3	1/2	13/18

This method comprehensively accounts for the integration scheme and element order, in addition to the element dimension, shape, and crack orientation. However, it is important to note that the validation of this method was confined to 2D problems (Slobbe et al., 2013).

He’s method

To achieve results that are more independent of mesh size and account for mesh irregularities, He et al. (2019) introduced two separate crack bandwidths for tensile and compressive behaviors based on Govindjee’s method. Geometrically, these bandwidths are determined by projecting the midpoint of element edges in directions normal and tangential to the crack, respectively. The mathematical expressions for these calculations are as follows:

$$h_H^i(\mathbf{x})=|l_i^{max}-l_i^{min}|(i=t\text{or}c)$$

(78)

with

$$\begin{array}{rl}{l}_t^{max}& =max[{\mathbf{x}}_N\cdot \mathbf{n}(\mathbf{x})]\\ l_t^{min}& =min[{\mathbf{x}}_N\cdot \mathbf{n}(\mathbf{x})]\\ l_c^{max}& =max[{\mathbf{x}}_N\cdot \mathbf{t}(\mathbf{x})]\\ l_c^{min}& =min[{\mathbf{x}}_N\cdot \mathbf{t}(\mathbf{x})]\end{array}$$

(79)

where $i=t\text{or}c$ serves as the index representing whether the condition is tensile ( $t$ ) or compressive ( $c$ ). $l_i^{max}$ and $l_i^{min}$ represent the maximum and minimum projection values, respectively. ${\mathbf{x}}_N$ denotes the midpoint of the edge for linear elements. $\mathbf{t}(\mathbf{x})$ is defined as the unit vector in the direction tangential to the crack at the point $\mathbf{x}$ .

The validation conducted by He et al. (2019) demonstrates that He’s method effectively achieves objective and mesh size-independent outcomes. It is particularly adept at solving problems characterized by irregular meshes and arbitrary crack orientations, due to its comprehensive consideration of element shape, size, and orientation with respect to crack directions and spatial positions. Despite its proven efficacy in 2D problems, it is important to note that the method was specifically developed for 2D linear elements. Consequently, its applicability to 3D problems and higher-order elements has yet to be validated.

Method for higher-order elements

In recent years, higher-order beam theories based on Carrera unified formulation (CUF) have arisen significant interests in various engineering fields (Arruda et al., 2018; Kaleel et al., 2017; Shen et al., 2022; Trombini et al., 2024). Within the CUF framework, the 3D displacement field of a beam model can be derived through various cross-sectional expansions, which is formulated as equation (80). This approach enables beam models to achieve 3D accuracy while preserving computational efficiency .

$$\mathbf{u}(x,y,z)=F_{\tau }(x,z){\mathbf{u}}_{\tau }(y),\tau =1,2,\dots ,M$$

(80)

where $F_{\tau }(x,z)$ represents the expansion functions on the $x-z$ plane, while ${\mathbf{u}}_{\tau }(y)$ denotes the unknown displacement vectors along the longitudinal axis, $y$ . The symbol $\tau$ serves as a summation index and $M$ indicates the total number of terms in the expansion functions.

Another significant advantage of CUF-based beam models is that there is no need for ad hoc assumptions commonly found in beam theory. Classical beam theories, such as Euler–Bernoulli and Timoshenko, can be integrated as special cases within the CUF framework. This integration enhances a comprehensive approach to beam modeling, where FE approximation can be utilized to derive solutions:

$$\mathbf{u}(x,y,z)=F_{\tau }(x,z)N_i(y){\mathbf{u}}_{\tau i},i=1,\dots ,N_{\mathrm{NE}}$$

(81)

where $N_i$ represents the shape function corresponding to beam node $i$ . $u_{\tau i}$ is the nodal displacement vector to be determined. $N_{\mathrm{NE}}$ denotes the total number of nodes within each beam element.

Equation (81) demonstrates that the choice of expansion functions and the order of beam elements in CUF-based beam models are independent variables, serving as configurable inputs. As depicted in Figure 10, a quadratic finite beam element can be expanded using quadratic Lagrange-like polynomials to approximate the 3D displacement field. This approximation process incorporates the beam shape function $N_i$ to describe displacement along the longitudinal axis of the beam, while the expansion function $F_{\tau }$ is employed to model kinematics across the cross-section. These operations are based on the nodal displacement vector ${\mathbf{u}}_{\tau i}$ at each Lagrange point. The physical meaning of ${\mathbf{u}}_{\tau i}$ is clearly illustrated in Figure 10. The utilization of Lagrange expansion offers the flexibility to accommodate any arbitrary cross-section. More details of the CUF-based model can be found in Carrera et al. (2014).

Figure 10.

Illustration of a quadratic Carrera unified formulation (CUF)-based finite beam element with Lagrange expansion function.

The utilization of CUF-based beam models with FE solutions for simulating strain-softening behavior in quasi-brittle materials inherently faces the challenge of mesh size dependency. To address this issue with fracture energy regularization, especially in the context of higher-order elements, various methods have been proposed to estimate crack bandwidths accurately. The concept of cubic root of taking the cubic root of a volume was introduced in Kaleel et al. (2018) and Nagaraj et al. (2020). This approach is inspired by Bažant and Oh (1983), where the method of utilizing the square root of an element area was suggested. The mathematical expression for this volume-based estimation is presented as follows:

$$h^{*}=(V_{\mathrm{GP}}{)}^{1/3}$$

(82)

where $V_{\mathrm{GP}}$ represents the volume attributed to each GP. It is calculated by evenly distributing the total volume of an element among all its GPs.

Building on the same concept of utilizing the cubic root of a volume, Shen et al. (2023a) proposed a modification for beam FEs that incorporate Lagrange expansions within the CUF:

$$h_1^{*}={}^{\mathrm{a}}\sqrt{\frac{A_e\times L_e}{(\sqrt{M}-1{)}^2\times (N_{\mathrm{NE}}-1)}}$$

(83)

with

$$\frac{\sqrt{A_e}}{L_e}\approx \frac{\sqrt{M}-1}{N_{\mathrm{NE}}-1}$$

(84)

where $L_e$ represents the length of the beam element, and $A_e$ denotes the area of a discretized mesh on the cross-section of the beam.

Figure 11 illustrates the concept of equation (83), based on the example from Figure 10. In this illustration, the total volume is divided into small volumes according to the order of beam element and expansion function. Each of these smaller volumes may contain one or more GPs. The characteristic element length for GPs is computed by the cubic root of the volume containing the target GPs. A key condition stated in equation (84) is that the geometry of the volume containing the target GP should approximate a cubic form. While equations (82) and (83) are based on the same concept, equation (82) tends to yield a higher estimation of the characteristic element length in comparison to equation (83).

Figure 11.

Illustration of calculation for $h_1^{\star }$ .

The constraints imposed by equation (84) often result in overly dense mesh requirements, diminishing the inherent efficiency and flexibility of CUF-based models. Moreover, the crack bandwidth may evolve through the softening period, requiring a more flexible approach. To address this, a novel method inspired by earlier projection techniques was introduced in Shen et al. (2023b). This approach enables the determination of an appropriate crack bandwidth that accounts for mesh characteristics and crack orientation, free from the restriction of equation (84). The formulation of this method, presented in equation (85), similar to the one from He et al. (2019), yet it is expressly adapted for use with higher-order beam elements.

$$h_2^{*}(\mathbf{x})=|max[{\mathbf{x}}_n\cdot \mathbf{n}(\mathbf{x})]-min[{\mathbf{x}}_n\cdot \mathbf{n}(\mathbf{x})]|$$

(85)

where ${\mathbf{x}}_n$ is the middle point of the edge in a 3D solid volume that contains the targeted GP, as illustrated in Figure 12. Notably, this 3D volume, determined as the same method from Figure 11, is a portion of the volume assembled by a beam element with Lagrange expansion. This method allows for the consideration of reduced integration by implicitly adjusting the volume, showcasing its adaptability to various critical computational factors such as element order, integration scheme, mesh size, and crack direction. While some mesh size independent results have been observed with this method, mesh sensitivity that higher-order element might cause less brittle response has also been reported in Shen et al. (2023b).

Figure 12.

Illustration of calculation for $h_2^{\star }$ .

The observation of a more ductile response in some models can be attributed to the fact that the method in Shen et al. (2023b) tends to underestimate the crack bandwidth when employing certain high-order elements. To mitigate the resulting inappropriate mesh sensitivity associated with these higher-order elements, a strategy similar to those discussed in Jirásek and Bauer (2012) and Slobbe et al. (2013) has been integrated into the projection method in Shen et al. (2024). This approach, designed to enhance accuracy and reduce sensitivity, is detailed as follows:

$$h_3^{*}(\mathbf{x})=\alpha \times h_G(\mathbf{x})$$

(86)

where $\alpha$ accounts for the element order, defined as the ratio of the combined weights of the softening GPs to the overall weight of all GPs present in a beam element.

The geometrical explanation behind equation (86) is depicted in Figure 13, utilizing the same example illustrated in Figure 10. Initially, the projection method is applied, followed by determining the parameter $\alpha$ for each GP, which depends on the status of surrounding GPs. For instance, employing a full integration scheme in Figure 13 results in 27 GPs, organized into nine groups along the beam axis. Within each group, three possible scenarios can occur: (i) all points are undergoing softening, resulting in $\alpha =1.0$ ; (ii) softening affects two-thirds of the points, yielding $\alpha =13/18$ from Slobbe et al. (2013); and (iii) all points are either unloading or in elastic loading, leading to $\alpha =0$ . The key advantage of this method lies in its ability to adapt the crack bandwidth according to the specific conditions of the GPs.

Figure 13.

Illustration of calculation for $h_3^{\star }$ .

Comparison of different estimation methods

A comparative analysis of the methods previously discussed for nonsquare quadrilateral elements is illustrated in Figure 14. Figure 14(a) visually compares these methods by projecting the element onto a central point, with superscripts $c$ and $t$ denoting crushing and tensile crack bands, respectively. Oliver’s second method is different by first calculating the equivalent crack length ( $l_{O2}$ ) and then obtaining $h_{O2}^t$ by dividing this length by the area of the element. Interestingly, the determined equivalent crack length ( $l_{O2}$ ) aligns with the crushing band length ( $h_H^c$ ) proposed by He et al. (2019).

Figure 14.

Comparison of each estimation method for a rectangular quadrilateral element with 1 mm width and 1.5 mm height: (a) geometrical meanings of each method and (b) obtained tensile characteristic element length.

Figure 14(b) showcases the estimated lengths of each model across various crack directions. Notably, Oliver’s first method exhibits discontinuity due to a binary crack indicator change when the crack aligns diagonally from the top left to bottom right within the element, whereas the other methods display continuous curves. Oliver’s second method does not consistently yield the same results as his first method. Compared to other projection methods, He’s method generally yields lower estimates.

The comparative analysis of various methods for estimating the characteristic element length in higher-order beam elements is depicted in Figure 15. Specifically, Figure 15(a) illustrates a 3D volume derived from a quadratic beam element utilizing quadratic Lagrange expansion, where the beam element has a length of 2 mm and a square cross-section with 1 mm sides. Among the methods reviewed, only Govindjee’s method is identified as appropriate for 3D applications, with its geometric interpretation depicted in Figure 15(a). Values of $h_3^{\star }$ in Figure 15(b) are derived from $h_G$ by multiplying a factor of 13/18, as presented in equation (86). In a 3D context, the crack direction is influenced by both azimuthal and elevation angles. Comparative analysis, as shown in Figure 15(b), reveals that $h_3^{\star }$ consistently yields higher estimations than $h_2^{\star }$ , offering an effective means to overcome the ductile behavior observed in models using $h_2^{\star }$ .

Figure 15.

Comparison of each estimation method for a beam element with Lagrange expansion: (a) geometrical meanings of each method and (b) obtained tensile characteristic element length.

Final remarks

The strain softening behavior observed in quasi-brittle materials leads to significant computational challenges in standard CDM, such as the loss of ellipticity and an undesirable mesh dependency. These issues further lead to the snapback behavior and zero fracture energy dissipation. Such critical problems are attributed to the mathematical defect rather than the numerical one. The mesh dependency to the discretization thus occurs for any discretization method, including mesh-free method (Pamin et al., 2003) and boundary element method (Lin et al., 2002; Sladek et al., 2003). Therefore, implementing regularization techniques is urgent to enhance the stability and reliability of numerical models dealing with materials that exhibit softening.

Nonlocal models incorporate an internal material length parameter to prevent the FPZ from concentrating at a singular point or line, maintaining a more realistic spread of damage. Fracture energy regularization adjusts the stress–strain softening curve based on spatial discretization, ensuring the total energy dissipated in numerical models aligns with the actual physical fracture energy of the material. Viscous regularization introduces an artificial viscosity term to the model, implicitly accounting for the rate-dependent behavior of materials. This addition helps to control and decelerate the damage evolution within specific elements. These approaches are crucial for ensuring that mathematical and numerical models accurately capture the energy dissipated in the fracture zone of quasi-brittle materials.

While regularization techniques have significantly mitigated the mesh-dependent response in numerical simulations of quasi-brittle materials, certain numerical defects remain challenging. In light of this, this work reviews recent advancements in regularization techniques, with a particular emphasis on refining the definitions of internal length in nonlocal models and characteristic element length in fracture energy regularization. Such advancements are crucial for accurately modeling the response of quasi-brittle materials, ensuring that simulations closely mirror physical behaviors.

In nonlocal models, adopting an evolving internal length concept has proven to be effective in addressing challenges related to inaccurate damage representation and the distribution of spurious damage. This work acknowledges the complexity of categorizing these approaches. For instance, the evolution of internal length within the eikonal nonlocal formulation and localizing gradient models is related to the damage state, a characteristic that could also align with damage-based models. Identifying the most effective regularization technique is challenging, as the efficacy of each model is influenced by multiple factors. While certain methods may be phenomenologically compelling and effective, establishing a well-defined thermodynamic background is crucial to ensure the thermodynamic consistency of the models.

The determination of characteristic element length in fracture energy regularization is crucial for achieving a mesh-independent global structural response. It involves integrating multiple aspects related to the FE and crack direction. This work presents several estimation methods tailored for CUF-based models, enhancing their regularization capabilities. However, these methods have been developed and validated specifically for higher-order beam elements. Their effectiveness for other types of elements, such as plates, shells, or standard 3D brick elements, may be an area for future exploration.

While the refinement of characteristic element length represents a significant step forward in improving regularization capacity, recent advancements have introduced the smooth crack band model, as proposed by Zhang and Bažant (2023). This innovative approach may offer an alternative way to address the challenges associated with fracture energy regularization in numerical modeling.

Viscous regularization is rarely applied within pure damage models. More commonly, a viscous term is introduced and integrated with other regularization methods to obtain a more stable solution with a faster convergence speed.