Smooth Lagrangian Crack Band Model Based on Spress-Sprain Relation and Lagrange Multiplier Constraint of Displacement Gradient

In 1962, Clough [1–3] pioneered finite element analysis (FEA) of continuum bodies. Using his variationally derived constant-strain finite element (FE) of an elastic continuum [3], he launched computational fracture mechanics by his finite element solution of the overall stress state in the cracked Norfolk Dam (Fig. 1), which mesmerized the structural engineering community of that time. Yet, after the lapse of 61 years, no completely satisfactory computational model for fracture exists today. A model in the preceding 2023 paper [4], which introduced the concept of sprain, came close except for severe programming and computational obstacles due to the fact that sprain forces had to be applied on the nodes adjacent to a finite element. In this paper, an improved model, namely Smooth Lagrangian Crack Band Model (slCBM) overcomes these challenges by introducing the new concept of spress as a force variable work-conjugate to the sprain, and by constraining the actual displacement gradient calculated from the nodal displacement vectors to the approximate gradient through the use of Lagrange multipliers. By doing so, finite elements with linear shape functions for both the displacement vector and for the approximate displacement gradient tensor can be used, thus avoiding the need of applying sprain forces on adjacent nodes.

By extensive comparisons with numerous experiments, a recent study in 2022 [5] demonstrated deficiencies of the recently popular computational fracture models, particularly the peridynamics and phase field. Relying on 11 different distinctive classical fracture experiments, including the results of the recently developed gap test [6–8], that study [5] identified serious inherent limitations of peridynamics and of all line crack models including the linear elastic fracture mechanics (LEFM), cohesive crack model, extended finite element method (XFEM), and phase-field model. It further highlighted the misleading nature of model “verifications” by “nondistinctive” tests, i.e., tests that can be fitted closely by very different models [9], commonly abused in many computational papers.

The adoption of displacement vector curvature (or Hessian) [4] as a localization limiter was inspired by the recent discovery of the gap test [6,7] and by comparisons of the classical 1983 crack band model (CBM) with material characteristic length [10–12] to the phase-field models and peridynamics [5,9]. The present theory can also be regarded as a generalization of the classical strain-gradient models [13–21], which differ mainly by missing the material rotation gradient that is important for shear fracture. These theories were developed in the works of Bažant and Jirásek and Gao et al. [13,22], and culminated in the very effective mechanism-based strain-gradient theory of Huang et al. [23], whose objective was hardening plasticity and the effect of geometrically necessary dislocations, rather than damage and fracture. A strain-gradient theory for damage in concrete was formulated by Cusatis et al. [14,15]. For comparisons with experiments and for applications in structural design, it is essential to use a realistic damage constitutive relation, for which the microplane models are essential [24–27]. Achieving superior comparisons with distinctive tests of concrete in Ref. [5] was enabled by the use of not only crack band model [10–12] but also the microplane model M7. The present use of displacement vector curvature also implies significant changes for continuum mechanics of damage in heterogeneous materials, including the homogenization theories [22,23,28–35].

A more detailed review of the subject, including the history, was presented in 2023 in the same journal [4], along with a broader historical review from the same perspective as here. Therefore, it would be superfluous to duplicate it.

The computations in the previous study [4] used, instead of s_ijk, the sprain force f_n = b A_n ∂Φ/∂u, which is the derivative of the sprain energy with respect to the nodal displacement vector u multiplied by the tributary material volume b A_n of node n in a regular finite element mesh (b = body thickness, A_n = tributary area). However, this approach required applying self-equilibrated forces on nodes adjacent to a given finite element, which turned out to be very cumbersome for programming, especially in commercial programs such as abaqus, and extended the running time by almost two orders of magnitude. For computations, it is advantageous to introduce a second-order tensor field, ζ, of components ζ_ij, to serve as an approximation of displacement gradient field u_i,j (dimensionless). Considering the fields of u_i and ζ_ij as independent, each of them may be now represented by C₀ finite elements with linear shape functions, which will greatly simplify FE analysis. To introduce the idea, let us first give a simple example.

Physical Motivation: Looking from a mechanics point of view, one could often find a physical motivation for the Lagrange multipliers [41,42]. Note the similarity of the last integral in Eq. (11) and right-hand side of Eq. (12). In the latter, f(x) · u represents the external work done on displacement u by a constant gravity force f (independent of u) per unit volume. So, in the former, λ (dimension Pa) can be seen as a generalized moment per unit volume doing external work on the gradient difference, like some constant pseudo-gravity (independent of the gradient difference). This external work may be likened to the applied force term in the potential energy, i.e., $−∫Vf(x)⋅udV$⁠. For u, the minimizing condition is ∂W/∂u = 0. By analogy of pseudo-gravity, this explains physically the last extremum condition ∂W/∂λ = 0. The fact that the work of λ is external, not part of the internal energy density, is important. This work is only added to the definition of the internal energy expression in Eq. (11) as a constraint condition in a weak sense to couple the strain energy and the sprain energy as explained earlier.

Why not a quadratic penalty? Indeed, one could impose it instead of using Lagrange multiplier. However, a penalty (with stiffness K_p) would have to consist of a quadratic form $Kp(∇u−ζ)2/2$⁠, which would have the physical meaning of adding fictitious additional strain energy to the system. If the quadratic penalty is not minimized exactly to zero, one would end up with a system containing some non-existent excessive strain energy. The “beauty” of the Lagrange multiplier method is that it represents a fictitious externally applied moment density (or pseudo-gravity) and adds nothing to the internal energy content of the structure.

The spress-sprain relation is based on the sprain energy density whose only purpose is to resist the curvature of the displacement field. This relation is a material property but it is not a constitutive relation in the standard sense because it deals with the second gradient of displacement. This was done in a rather hypothetical way in the preceding study [4], in which it was assumed that what is here named the spress is controlled through a constant threshold C, after which each component of the spress tensor grows from zero with a sprain stiffness κ that decreases with the inelastic volumetric strain $εV$⁠. Better understanding of the spress-sprain relation will necessitate discrete particle simulation of sprain.

The threshold C and the sprain stiffness κ are here considered as constants, although it is suspected that the sprain stiffness κ should be considered as a function of either the total volumetric strain, $εV$ or of its inelastic part $εV″$ (as in Ref. [4]), or |ξ|. The expression for s_ijk in Eq. (14) is smooth, which might help the convergence of the iterative scheme. A high exponent such as n = 10 may have to be used so that the sprain would have a negligible effect in the elastic regime.

After further modeling experience, which showed that decreasing the sprain stiffness κ requires further investigation (and that a sharp nonzero threshold C might be unnecessary and might even slow the convergence rate), we consider the relation shown in Eq. (13) with C = 0 and a constant κ, i.e., s_ijk = κξ_ijk. Taking C = 0 might be a reasonable simplification because, as argued here, the displacement curvatures in the elastic range are negligible and the sprain energy is not big enough to have a substantial influence on the strain energy and to change the elastic behavior of the structure. However, this may not hold if the structure contains defects that cause stress concentration.

This paper is focused on the computational aspects of the sprain theory, while the physical aspects with the comparisons to the results of various types of fracture test are relegated to a subsequent article. Computer simulation results for the uniaxial tension of a 1D bar and the three-point bend fracture test are presented.

Consider a one-dimensional bar of length L = 100 mm using a material following a simple bilinear elastic softening stress–strain relation with a Young’s modulus of E = 30 GPa and a tensile strength of f_t = 4 MPa (Fig. 4). The corresponding elastic strain at the tensile strength limit is $ε0=ft/E$ and the strain after complete softening, $εc$⁠, is chosen such that $εc/ε0=3$⁠. The characteristic length l₀ is 20 mm. In the original CBM, the crack bandwidth is equal to the characteristic length l₀ and the strain, $ε=u′$⁠, is uniform across the band. Here, l₀ is subdivided into several uniform-strain elements and the width of the crack band is roughly equal to the characteristic length. Specifically, the bar is discretized into 5, 21, 51, 101, 201, and 501 elements, and the characteristic length l₀ is thus divided into 1, 5, 11, 21, 41, and 101 elements, respectively. The middle element of the bar will have a slightly reduced tensile strength of 0.99 f_t, in order to induce strain localization at the middle of the bar. The spress-sprain relation used is based on the Eq. (13) with the threshold C = 0 and the sprain stiffness κ = 127.5 GPa.

This is a static problem, and an implicit static algorithm is used. To capture the snapback behavior, a simple algorithm inspired by the crack mouth opening displacement control is used. First, since it is known that the middle (weakest) element is always stretched after each increment, a pair of displacement increments equal in magnitude but opposite in direction is prescribed at the two nodes of the middle element. After that, the nodal displacements at the two boundaries (along with the rest of the DOF) will be obtained using standard Newton–Raphson algorithm with the conditions that the reaction forces at the two nodes of the middle element equal to zero.

The displacement field, strain field, and the approximate sprain field which is the gradient of ζ of the bar are discretized into different numbers of elements, as shown in Fig. 6. The original CBM is used when l₀ is equal to the width of the element (Figs. 6(a), 6(b), and 6(c)). As the bar is discretized into more elements, through the use of the localization resisting sprain energy, the damage gets spread over the characteristic length l₀ with a smooth bell-type distribution using a proper chosen sprain stiffness κ. Figure 5 shows the force versus displacement curves as the bar is discretized into different numbers of elements. Because the strain field is uniform during elastic loading, the elastic response is the same with or without using the sprain energy; therefore, setting the threshold C to be zero is justified. Note that the localization resistance should decrease as the damage increases and vanish when the material is totaly broken, since there will be residual stresses in the bar if a constant sprain stiffness κ is used for post-peak damage [4]. The simulation using a constant sprain stiffness κ is thus terminated near post-peak damage before the residual stresses become large enough to affect the force versus displacement curve.

The strain field (Figs. 6(b), 6(e), 6(h), 6(k), and 6(n)) and the approximate sprain field (Figs. 6(c), 6(f), 6(i), 6(l), and 6(o)) smooth and converge as the number of finite elements increases. The convergence can also be shown by the plot of the L₂ norm of the error of the reaction force of the system, i.e., L₂ = ‖F − F_asympt‖₂. Here F_asympt would be the exact loading force for the given load-point displacement if the exact value were known a priori. Since it is not known, we determine the value F_asympt so that the computed value in the log-log plot of Fig. 10(a) falls into a straight line, and a straight line in the log-log scale represents a power law. A linear regression fitting line shows that the rate of convergence is at least quadratic.

The specimen geometry for the three-point bend test simulation is similar to that used in Ref. [4] but the material model is the M7 microplane model for concrete. Note that different parameter combinations for M7 microplane model correspond to different material properties; it is an important process for objectively describing a quasi-brittle material, which is also detailed in Ref. [43]. Here the elastic (Young’s) modulus is E = 38.3 GPa, the Poisson ratio is ν = 0.18. The calibrated M7 (implicit) parameters were k₁ = 0.00015, k₂ = 110, k₃ = 30, k₄ = 100. The threshold C = 0 and the sprain stiffness κ = 1.9 GPa. The span-to-depth ratio is L/D = 3.8 and the notch depth is a/D = 0.197.

Figure 7 shows the simulation results when a fine square mesh in a zone surrounding the growing crack is deliberately inclined from the notch line at various angles equal to $0deg,5deg,15deg$⁠, and $30deg$⁠. The mesh size in the near-tip zone is l₀/20. Each element is a four-node finite element with linear shape functions for u, ζ and the Lagrange multiplier λ, as already described. Figure 7 compares the widths and shapes of the fracture process zone and demonstrates that the contours of constant transverse strains $εxx$ of different magnitudes are virtually identical. Similarly, the difference between the values of the center-span equilibrium loading force as a function of the vertical load-point displacement for different mesh inclinations is negligible (Fig. 8).

Figure 9 shows the effect of changing mesh size near the crack tip, considering the same identical three-point bend fracture specimens and simulation setup. The square element sizes are chosen as 5 mm, 2.5 mm, and 1 mm, respectively. The results are not only stable but also almost independent of the element size, overcoming the mesh bias. Figure 9 is a practical demonstration of convergence of the present method with linear shape functions and Lagrange multipliers. Figure 10(b) shows the plot of the L₂ norm of the error of the reaction force of the system. Similar to the 1D simulation, here F_asympt is determined as the force value for which the L₂ norm of the error of the reaction force falls into the straight line in the log-log plot, representing a power law. An approximate power law must be expected because of the nearly perfect similarity of the three meshes of different refinements.

1. Computer implementation of the previously proposed concept of limiting damage localization by applying sprain forces on adjacent nodes of a finite element [4] runs into formidable obstacles. It is tedious to program and extends the running time on commercial FE codes such as abaqus by almost two orders of magnitude.

2. What helps to eliminate this problem is to consider both the displacement vector u with components u and v and the approximate gradient tensor ζ with components ζ^u,x, ζ^u,y, ζ^v,x, ζ^v,y as independent nodal variables.

3. The approximate gradient components are constrained to the actual gradient components u_x, u_y, v_x, v_y by means of Lagrangian multipliers λ^u,x, λ^u,y, λ^v,x, λ^v,y applied to the difference between the actual and approximate gradients. Zeroing of the derivatives of the Helmholtz free energy of the structure with respect to the Lagrange multipliers of all the nodes gives the optimum constraint yielding the best approximation of the actual displacement gradient tensor.

4. The gradient tensor of the approximate gradient tensor represents the approximate curvature of the displacement vector. Their limitation by means of the spress-sprain relation based on the sprain energy prevents spurious damage localization at, and behind, the fracture front of finite width. In the generalized stiffness matrix, this is effected by the matrix components representing spress-based resistance to sprain.

5. The Lagrange multipliers can be considered as element properties or the nodal properties. The former is more efficient by reducing the total number of DOF, while treating the Lagrange multipliers as nodal properties aid the programming in commercial FEA software such as abaqus.

6. Thanks to adopting approximate displacement gradients as independent nodal variables and constraining them by Lagrangian multipliers, there is no need for localization control by sprain forces applied on the nodes of the adjacent elements.

7. Physically, the Lagrange multipliers represent applied generalized moments per unit volume working on the gradient differences. It is important that their work is not part of the internal Helmholtz free energy stored in the material but is external, like the work of gravity forces.

8. In two-dimensional fracture modeling, there are ten unknown degrees-of-freedom per node: two nodal displacement components, four approximate displacement gradients, and four Lagrange multipliers. The increased number of variables requires introducing into a user’s element of abaqus a larger and more complex stiffness matrix.

9. A severe test in computer simulations is the mesh bias and convergence when the material characteristic length is subdivided into more elements. FE simulations with abaqus show both tests to be verified. If the material characteristic length l₀ is subdivided into five elements or more, the directional bias of a square mesh is totally eliminated.

10. In contrast to the CBM, the slCBM can simulate the variation of the crack band width due to loading and to crack growth. Importantly, the width of the crack band (or damage zone) remains nearly the same for different subdivisions. Mesh rotations by 5 deg, 15 deg, 30 deg make no discernible difference. Although the mesh orientation independence was demonstrated in the preceding study, it was at the cost of great programming difficulties as well as huge computational burden and enormous running time.

Partial financial support under NSF (Grant No. CMMI-202964) and ARO (Grant No. W911NF-19-1-003), both to Northwestern University, is gratefully acknowledged.

There are no conflicts of interest.

The program used to generate the computer results is posted on Z. P. Bažant’s website.