Variational damage model: A new paradigm for fractures

Abstract

The computational modeling of fracture in solids using damage mechanics faces challenges with complex crack topology. This can be addressed by using a variational framework to reformulate the damage mechanics. In this paper, we propose several mathematically elegant variational damage models (VDMs) for fracture mechanics without explicitly using damage variables. Based on the energy density ϕ, the fracture energy density is formulated as \tilde{\varphi } = ϕ/(1 + ℓϕ/G_c) and the damage variable is expressed as s = ϕ/(ϕ + G_c/ℓ), which satisfy \tilde{\varphi }|_ϕ→∞ = G_c/ℓ and s|_ϕ→∞ = 1 as ϕ approaches infinity. These limits demonstrate that the new energy density converges to the Griffith energy release rate at full-damage state. The VDM profoundly modified the energy functional, implicitly incorporating the damage field. As a generalization of previous model, we propose a family of VDMs of varying orders. Additionally, we develop a multi-damage model to account for different types of energy densities, such as elastic thin plate and gradient elasticity. Using this functional, it is straightforward to deduce the governing equation for automatically evolving fractures. These formulations can be employed in conventional finite element method or other numerical methods with minimal modifications. Compared to the phase field method with the same mesh density, a sharper crack interface can be achieved. We demonstrate the capabilities of the proposed variational damage formulations using representative numerical examples.

Graphical abstract

Keywords

fracture / damage mechanics / multi-damage / variational principle

1 Introduction

The phenomenon of fracture evolution can be attributed to the interplay between crack topology and the distribution of deformation energy. Fracture can result in severe and disastrous failures of structures and components and is thus of major concern in engineering and material science. The precursors of fracture are commonly cracks which nucleate, propagate and interact resulting in a notable substantial topological change, which is accompanied by the presence of stress singularities at the crack tip/front. In the past decades, numerous models and computational methods for fracture have been developed. Among those models and methods are the extended finite element method (XFEM) [1–3], meshless methods such as the reproducing kernel particle method [4] or the element-free Galerkin method [5], the cracking particle methods (CPM) [6,7], the cracking element method [8], remeshing techniques [9,10], cohesive elements [11,12], phase field methods [13–16], peridynamics (PD) [17,18], nonlocal theories [19,20] and various others.

These numerical methods can be categorized into two main groups: topology-modifying methods and field-modifying methods. Direct modification of topology, such as remeshing in the XFEM, is computationally challenging, especially in three dimensions (3D) and for complex crack patterns including crack branching and interactions. To our knowledge, such complex scenarios require ad hoc schemes to determine, for example, whether a crack will branch. Topology-modifying methods like CPM can naturally capture these phenomena, but this comes at the cost of the accuracy of the crack’s topology. The phase field method, as one of the field-modifying methods, employs an additional field to represent the spatial location of fractures. The proposed approach is astute, as it guarantees that the degradation of strength is mathematically continuous. Without artificial perturbation, the fracture fields exhibit enhanced stability and a smoother pattern.

The basis of fracture modeling relies on the description of fracture evolution [21,22]. The phase field technique utilizes a volume representation of the surface through the application of Γ-convergence [23], which is achieved by updating the phase field in a variational framework through the solution of an additional differential equation. The phase field approach has been applied to a range of complex engineering problems [24–26]. However, the phase field leads to a multi-physics problem, even for a purely mechanical problem resulting in increased computational expenses. Furthermore, in order to adequately represent the crack width, a highly dense mesh in the crack region is necessary.

Several authors have ‘linked’ phase field models to (nonlocal) damage models [27], which do not require the solution of an additional differential equation and the introduction of another primary variable. Damage mechanics has proven to be quite effective in many practical applications. The damage variable is commonly defined as the decrease in the cross-sectional area due to the evolution of fracture. Mathematically, a damage variable for an isotropic homogeneous solid, which we designate by s, is introduced as

(1)s_{\text{v}} =1-v^{\ast }/v, v\in \{ E,v,\mu ,K\} ,

where v denotes the original value in the absence of damage, and v^{\ast } represents the effective value after the occurrence of damage; E is the elasticity modulus, v is the Poisson ratio, µ is the shear modulus and K is the bulk modulus. The principles in damage mechanics show some similarities to plasticity theory which studies dissipative processes. The establishment of the variational nature of the dissipative process remains incomplete, necessitating additional consideration such as the enforcement of irreversibility conditions [28] and the principle of maximum dissipation [29].

Conventional continuum damage mechanics (CDM) follows the routine of plasticity mechanics, where the damage field is set locally and its evolution is determined by the Kuhn-Tucker loading−unloading conditions. On the other hand, the phase field approach is a theoretical framework that encompasses both the displacement field and phase field into a single energy functional. The determination of the phase field is finally achieved through the solution of a differential equation, which requires the gradient of the phase field. For this reason, the phase field is commonly related to gradient damage models. By borrowing certain concepts from the phase field model, this research aims to present a variational damage model (VDM) for fracture modeling in solid. The proposed model offers an effective variational framework for studying the failure of quasi-brittle solids. The energy density is significantly altered in this model through the concept of damage variables as a function of energy density.

The subsequent sections are structured as follows. Section 2 provides a brief review of damage models. Section 3 introduces the damage approximation of crack topology using volume representation, without considering the influence of surface density gradients. In Section 4, we construct the governing equations of fracture in solids, using elasticity theory as an illustrative example, and present the accompanying boundary conditions. Section 5 outlines the development of a variational framework for damage mechanics, establishing simple governing equations that focus solely on the mechanical field to model damage evolution. Section 6 generalizes the first-order VDM to arbitrary orders and thoroughly examines the features of the one-dimensional (1D) case. Section 7 presents a theoretical framework that extends existing damage models to incorporate multi-damage patterns across different physical domains. Section 8 discusses the implementation of this approach using the finite element method (FEM). In Section 9, we demonstrate the model’s capabilities through a variety of fracture problems in statics. Finally, the concluding section draws several conclusions and proposes prospective avenues for future research.

2 Review of damage mechanics and gradient damage model

2.1 Continuum damage mechanics

The CDM framework offers a sophisticated approach to simulating material degradation, particularly relevant in fracture modeling [30–34]. In scalar-valued CDM, the damage variable s serves to quantify the degree of material damage, ranging from 0 for intact material to 1 indicating complete failure.

Within this framework, the stress state is adjusted to account for material damage, as depicted by the following relationship:

(2)\text{σ} =(1- s){\mathbb{D}}^{\text{e}} :\text{ε}

where \text{σ} represents the stress tensor, s is the damage variable, {\mathbb{D}}^{\text{e}} denotes the fourth-order elastic stiffness tensor, and \text{ε} is the strain tensor. This expression captures the degradation of material stiffness due to evolving microstructural imperfections [35].

To model the damage evolution, CDM introduces a damage loading function,

(3)f=f( \tilde{\epsilon },\kappa )=\tilde{\epsilon }-\kappa ,

where \tilde{\epsilon } is a scalar function derived from the strain tensor, encapsulating the strain’s effect on damage, and κ is the historical variable of the material. In the above damage evolution function, the simplest form is assumed. The evolution of damage is governed by the Kuhn-Tucker loading−unloading conditions: f\le 0,\dot{\kappa }\ge 0, \dot{\kappa }f= 0, ensuring the irreversibility of damage. The historical variable, κ, initiates at a predefined threshold, κ_i, and evolves as dictated by f = 0 to signify ongoing damage progression.

Furthermore, the relationship s = s(κ) governs the evolution of the damage variable, analogous to the plastic flow rule in plasticity mechanics. It is noteworthy that the aforementioned damage models utilize linear degradation functions, as opposed to nonlinear degradation functions in stress relations. For a comprehensive overview of damage models, please consult Ref. [27].

2.2 Regularized damage model

A notable example of a damage model with regularization, as described in Ref. [36], is based on the stored energy functional

(4)\mathrm{\Phi }(\mathit{u},s)={\int }_{\mathrm{\Omega }}(1-s)\varphi (\text{ε})- \mathit{g}\cdot \mathit{u}+ \frac{1}{2}{G}_{\text{c}} {\ell }^2\mathrm{\nabla }s \cdot \mathrm{\nabla }s \text{d}V,

where \mathit{g} is the gravity density, \mathit{u} is the displacement field, G_{\text{c}} is the critical energy release rate and ℓ is the crack length scale.

The variational derivation of the energy functional outlined above results in

(5)\begin{array}{c}\mathrm{\nabla }\cdot \left((1-s){\displaystyle \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\text{ε}}}\right)+\mathit{g}=0,\\ \varphi +G_{\text{c}}{\ell }^2\mathrm{\nabla }\cdot \mathrm{\nabla }s=0.\end{array}

The second equation above constitutes a partial differential equation, depicting the evolution of the damage field akin to a diffusion process, where the potential energy density serves as the “heat” source term. Evidently, the damage field spans the entire domain, and achieving damage localization necessitates additional constraints. Therefore, during the numerical implementation process, it becomes imperative to incorporate the dissipation distance function

(6)\mathcal{D}\left(s_1,s_2\right)=\{\begin{array}{l}{{\displaystyle \int }}_{\text{Ω}}{G}_{\text{c}}\left(s_2-s_1\right)\text{d}V,\text{if}{s}_2\ge s_1\text{in Ω},\\ +\mathrm{\infty },\text{otherwise}\text{.}\end{array}

In the study by Jirásek and Zeman [36], the analysis of the aforementioned damage fields was confined to 1D scenarios. The research employs intricate mathematical derivations to qualitatively analyze the evolution process of these damage fields. However, extending these analyses to 3D scenarios is challenging, with limited validation in such cases. The numerical solution of conventional variational damage methods necessitates specialized techniques, such as local iteration. Cracks exhibit diffusive behavior, and their distribution is extensive. Significant disparities exist between the physical reality of cracks and their representation in models, given that cracks lack thickness and only manifest in localized regions––a phenomenon commonly referred to as localization.

3 Damage approximation of crack topology

3.1 Motivation: one-dimensional bar with a crack

Let us consider an infinitely long bar with a cross-sectional area denoted by Γ. The bar is situated in a domain Ω = Γ × L with L = [−∞,+∞]. The location along the axis of the bar is denoted by x, which belongs to the interval L. It is assumed that a crack occurs at the axial location x = 0 in a bar. The sharp crack topology can be mathematically represented by introducing an auxiliary field variable, denoted as s(x), which takes values in the interval [0,1]. The function s(x) is defined as follows:

(7)s(x):=\{\begin{array}{ll}1,& x=0,\\ 0,& \text{otherwise}.\end{array}

The undamaged state is denoted by s = 0, while the fully-damaged state of the material is denoted by s = 1, as seen in Fig.1(a). Nevertheless, it is important to note that a sharp fracture can be characterized as a surface with a negligible thickness, hence exhibiting a distinct topological difference when compared to volumetric structures. Considering the aforementioned, we can approximate the non-continuous damage field represented by Eq. (7) as a piecewise function denoted by:

(8)s(x):=\{\begin{array}{ll}1,& x\in [-\ell /2,\ell /2],\\ 0,& \text{otherwise},\end{array}

as shown in Fig.1(b). The distribution of the sharp damage field is confined to a limited region, which is characterized by the length-scale parameter ℓ. As ℓ approaches zero, the function described by Eq. (7) is obtained. Hence, the utilization of a simple discontinuous piecewise function can be an effective means to represent crack.

The quadratic functional associated with the given piecewise function is assumed as

(9)I(s)={{\displaystyle \int }}_{\text{Ω}}{s}^2\text{d}V.

It is important to note that the power of surface density is not limited to 2, as 1^{\alpha }=1 for any α. This topic will be further explored in subsequent sections.

By employing the solution presented in Eq. (8) and considering the differential volume element dV as equal to Γdx, we may establish the following identification:

(10)I(s)=\ell \mathrm{\Gamma }.

Therefore, the relationship between the functional I and the crack surface Γ is established.

Consequently, the functional can be introduced

(11){\text{Γ}}_{\ell }(s):=\frac{1}{\ell }I( s)=\frac{1}{\ell }{{\displaystyle \int }}_{\text{Ω}}{s}^2\text{d}V.

Evidently, the functional yields the fracture topology as illustrated in Fig.1(b). The functional Γ_ℓ(s) can be seen as the fracture surface that is dispersed throughout a domain described by ℓ. By employing volume representation, it is possible to depict the sharp crack surface through the damage field s. The aforementioned attribute holds significant importance in the formulation of constitutive models for the analysis of fracture propagation through the utilization of damage fields.

In the multi-dimensional space, the crack functional denoted by damage field reads

(12){\text{Γ}}_{\ell }(s)={{\displaystyle \int }}_{\text{Ω}}\gamma ( s)\text{d}V,

where we introduce the crack surface density function per unit volume

(13)\gamma (s)=\frac{1}{\ell }{s}^2.

This function depends on the damage field only, and no spatial gradient is required. The crack surface density neglects the contribution of gradient, which contributes to a sharper description of interface.

3.2 Dissipation function for crack evolution

The evolution of cracks is driven by some energetic driving forces as outlined below. Because the cracking is fully dissipative in nature and the cracking process is irreversible, the growth of crack surface in time should satisfy

(14){\dot{\mathrm{\Gamma }}}_{\ell }:=\frac{\text{d}}{\text{d}t}{\mathrm{\Gamma }}_{\ell }(s(x,t) )⩾0.

More specifically,

(15){\dot{\mathrm{\Gamma }}}_{\ell }={{\displaystyle \int }}_{\text{Ω}} \dot{\gamma }\text{d}V={{\displaystyle \int }}_{\text{Ω}} \frac{\mathrm{\partial }\gamma }{\mathrm{\partial }s}\dot{s}\text{d}V\ge 0.

We observe that the global irreversibility constraint is satisfied locally by the positive evolution of damage field

(16)\frac{\mathrm{\partial }\gamma }{\mathrm{\partial }s}=\frac{2s}{\ell }\ge 0\text{and}\overline{s}\ge 0.

The later constraint is related to the non-reversible evolution of micro-cracks and micro-voids. The evolution of damage field can be achieved by a global constitutive dissipation functional as

(17)D(\dot{s};s)={{\displaystyle \int }}_{\text{Ω}}{G}_{\text{c}}\dot{\gamma }(\dot{s};s )\text{d}V,

where G_c is a constitutive threshold value related to the critical Griffith-type fracture energy [37].

4 Fracture in the elastic solid

4.1 Degradation of strain energy

For the small-strain deformation, the motion of fracturing solid can be described by displacement field \mathit{u}:= \mathit{u}(\mathit{x},t) with \mathit{x}\in \mathrm{\Omega } and t denoting time. The strain tensor for the geometrically linear theory is defined as

(18)\text{ε}(\mathit{u}):= \frac{1}{2}\left(\mathrm{\nabla }\mathit{u}+{\mathrm{\nabla }}^{\text{T}}\mathit{u}\right).

The global potential energy functional for elasticity reads

(19)E(\mathit{u},s)={{\displaystyle \int }}_{\text{Ω}}\varphi ( \text{ε}(\mathit{u}),s)\text{d}V,

where ϕ is the energy stored in the bulk of solid per unit volume. Similar to the phase field theory, the degradation of stored bulk energy by damage field can be isotropic degradation or anisotropic degradation. For the case of isotropic degradation due to damage field, the form is

(20)\varphi (\mathit{\epsilon },s) =g(s )\varphi (\mathit{\epsilon }),

where \varphi (\text{ε}) is a standard free energy function of an undamaged elastic solid

(21)\varphi (\text{ε})=\lambda \text{t}{\text{r}}^2(\text{ε})/2+\mu \text{tr}\left({\text{ε}}^2\right),

with elastic constants λ > 0 and µ > 0. The monotonically decreasing function g(s) reflects the degradation of strain energy because of evolving damage field by the properties

(22)g(0)=1, g( 1)=0 ,g^{\mathrm{\prime }}(1)=0 .

A simple function having these properties is

(23)g(s)=( 1-s{)}^2.

Then the evolution of stored energy becomes

(24)\dot{E}(\dot{\mathit{u}},\dot{s})={{\displaystyle \int }}_{\text{Ω}} \dot{\varphi }\text{d}V={{\displaystyle \int }}_{\text{Ω}} \left(\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }s}\dot{s}+\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\text{ε}}:\dot{\text{ε}}\right)\text{d}V,

with the constitutive expression

(25)\text{σ} :=\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\text{ε}}=(1-s{)}^2 {\text{σ}}_0 (\text{ε}),

in terms of stress of undamaged solid

(26){\sigma }_0(\epsilon ):=\lambda \text{tr}(\epsilon )\mathit{I}+ 2\mu \epsilon ,

where I is the identity matrix. On the other hand, the anisotropic degradation describes the energy by

(27)\varphi (\text{ε}, s)=g (s){\varphi }^{+}(\text{ε})+{\varphi }^{-}(\text{ε}),

based on the spectral decomposition of the stored energy of isotropic linear elasticity as

(28)\varphi (\text{ε})={\varphi }^{+}(\text{ε})+ {\varphi }^{-}(\text{ε}).

The spectral decomposition of strain tensor \text{ε} has the following properties

(29)\begin{array}{l}\varphi ={\displaystyle \frac{1}{2}}\text{σ}:\text{ε}={\displaystyle \frac{1}{2}}\left({\sigma }_i{\mathit{n}}_i\otimes {\mathit{n}}_i\right):\left({\epsilon }_i{\mathit{n}}_i\otimes {\mathit{n}}_i\right)={\displaystyle \frac{1}{2}}\sum _i{\sigma }_i{\epsilon }_i.\\ {\varphi }^{\pm }=\lambda {⟨{\epsilon }_1+{\epsilon }_2+{\epsilon }_3⟩}_{\pm }^2+\mu \left({⟨{\epsilon }_1⟩}_{\pm }^2+{⟨{\epsilon }_2⟩}_{\pm }^2+{⟨{\epsilon }_3⟩}_{\pm }^2\right).\end{array}

Based on above energy decomposition, the total stress tensor can be written as

(30)\text{σ} ={\text{σ}}^{+}+{\text{σ}}^{-},

and

(31)\begin{array}{rl}{\text{σ}}^{\pm }=& \frac{\mathrm{\partial }{\varphi }^{\pm }}{\mathrm{\partial }\text{ε}}=\lambda {⟨{\epsilon }_1+{\epsilon }_2+{\epsilon }_3⟩}_{\pm }\mathit{I}\\ & +2\mu \left({⟨{\epsilon }_1⟩}_{\pm }{\mathit{n}}_1\otimes n_1+{⟨{\epsilon }_2⟩}_{\pm }{\mathit{n}}_2\otimes {\mathit{n}}_2+{⟨{\epsilon }_3⟩}_{\pm }{\mathit{n}}_3\otimes {\mathit{n}}_3\right).\end{array}

where n_1,n_2,n_3 are the eigenvectors of the associated principal strains {\epsilon }_1, {\epsilon }_2,{\epsilon }_3 of \text{ε}\left(=\sum _{i=1}^3{\epsilon }_i{\mathit{n}}_i\otimes {\mathit{n}}_i\right).

4.2 Governing balance equations

For the displacement field, the exterior surface of the body Ω comprises of boundaries ∂Ω = ∂Ω_u∪∂Ω_t,∂Ω_u∩∂Ω_t = 0, and the corresponding Dirichlet condition

(32)\mathit{u}(\mathit{x},t)={\mathit{u}}_{\text{D}} (\mathit{x},t)\text{at} \mathit{x}\in \mathrm{\partial }{\mathrm{\Omega }}_{\text{u}},

and the Neumann-type boundary conditions with traction force t_{\mathrm{N}}(x ,t). The external mechanical loading is defined as

(33)P(\mathit{u}):={{\displaystyle \int }}_{\text{Ω}} \mathit{g}\cdot \mathit{u}\text{d}V +{{\displaystyle \int }}_{\mathrm{\partial }{\text{Ω}}_{\text{t}}}{\mathit{t}}_{\text{N}}(\mathit{x},t)\cdot \mathit{u}\text{d}A.

Based on the principle of virtual power, the rate of energy storage functional, the dissipation functional in Eq. (17) and the external power functional in Eq. (33) satisfy

(34)\dot{E}(\dot{\mathit{u}},\dot{s})+D (\dot{s})-P(\dot{\mathit{u}})=0.

For all admissible rates \dot{\mathit{u}} and \overline{s}, they meet the homogeneous form of Dirichlet boundary conditions,

(35)\begin{array}{rl}& \dot{\mathit{u}}\in {\mathcal{V}}_{\text{u}}:=\left\{\dot{\mathit{u}}\mid \dot{\mathit{u}}=\mathbf{0}\text{on}\mathrm{\partial }{\mathrm{\Omega }}_{\mathrm{u}}\right\},\\ & \text{and}\dot{s}\in {\mathcal{V}}_{\text{s}}:=\left\{\dot{s}\mid \dot{s}=0\text{on}\mathrm{\partial }{\mathrm{\Omega }}_{\text{s}}\right\}.\end{array}

Inserting this value and applying the Gauss theorem give

(36)\begin{array}{rl}& {{\displaystyle \int }}_{\text{Ω}}(-\left(\mathrm{\nabla }\cdot \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\text{ε}}+\text{g}\right)\cdot \dot{\mathit{u}}+\left(\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }s}+G_{\text{c}}\frac{\mathrm{\partial }\gamma }{\mathrm{\partial }s}\right)\cdot \dot{s}\text{d}V\\ +& {{\displaystyle \int }}_{\mathrm{\partial }{\text{Ω}}_{\text{t}}}\left(\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\text{ε}}\cdot \mathit{n}-{\mathit{t}}_{\text{N}}\right)\cdot \mathit{u}\text{d}A=0,\end{array}

where n is the outward unit normal direction of the Neumann boundary. Therefore, the damage model of fracture in elastic solids is

(37)\begin{array}{rl}& \mathrm{\nabla }\cdot \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\mathit{z}}+\mathit{g}=\mathbf{0},\\ & \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }s}+G_{\text{c}}\frac{\mathrm{\partial }\gamma }{\mathrm{\partial }s}=0,\end{array}

along with the Neumann-type boundary conditions

(38)\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\text{ε}}\cdot \mathit{n}={\mathit{t}}_{\text{N}}\text{on}\mathrm{\partial }{\text{Ω}}_{\text{t}}.

5 Variational framework of damage model

We use the surface density by Eq. (13) and formulate the total potential energy of the elastic body as

(39)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}g(s){\varphi }^{+}+{\varphi }^{-}+ G_{\text{c}}\frac{s^2}{\ell }-\mathit{g} \cdot \mathit{u}\text{d} V,

where {\varphi }^{\pm } =1/2{\text{σ}}^{\pm }:{\text{ε}}^{\pm }. It is worth noting that the representation of a sharp surface in this VDM is different from the phase field model. The phase field fracture model uses \mathrm{\Gamma }( d)=\left(d^2+{\ell }^2\mathrm{\nabla }d\cdot \mathrm{\nabla }d\right)/2\ell to describe the surface, where the gradient term of phase field d is compulsory. Current model uses the damage variable only without considering its gradient effect.

The variation of Φ reads

(40)\begin{array}{rl}\delta \text{Φ}=& {{\displaystyle \int }}_{\text{Ω}}\left(g(s)\frac{\mathrm{\partial }{\varphi }^{+}}{\mathrm{\partial }\text{ε}}+\frac{\mathrm{\partial }{\varphi }^{-}}{\mathrm{\partial }\text{ε}}\right):\delta \text{ε}+{\varphi }^{+}{g}^{\mathrm{\prime }}(s)\delta s\\ & +2G_{\text{c}}\frac{s}{\ell }\delta s-\mathit{g}\cdot \delta \mathit{u}\text{d}V\\ =& {{\displaystyle \int }}_{\text{Ω}}\left(-\mathit{g}-\mathrm{\nabla }\cdot \left(g(s)\frac{\mathrm{\partial }{\varphi }^{+}}{\mathrm{\partial }\text{ε}}+\frac{\mathrm{\partial }{\varphi }^{-}}{\mathrm{\partial }\text{ε}}\right)\right)\delta \mathit{u}\\ & +\left(g^{\mathrm{\prime }}(s){\varphi }^{+}+2G_{\text{c}}\frac{s}{\ell }\right)\delta s\text{d}V.\end{array}

For any \delta \mathit{u} and δs, δΦ = 0 yields

(41)\mathrm{\nabla }\cdot \left(g(s)\frac{\mathrm{\partial }{\varphi }^{+}}{\mathrm{\partial }\text{ε}}+\frac{\mathrm{\partial }{\varphi }^{-}}{\mathrm{\partial }\text{ε}}\right)+\mathit{g}= \mathbf{0},

(42)g^{\prime }(s){\varphi }^{+}+2G_{\text{c}} \frac{s}{\ell }=0.

In the case of g(s)=( 1-s{)}^2, Eq. (42) becomes

(43)\begin{array}{rl}& 2(s-1){\varphi }^{+}+2s{G}_{\text{c}}/\ell =0\to \\ & s=\frac{{\varphi }^{+}}{{\varphi }^{+}+G_{\text{c}}/\ell }=\frac{1}{1+G_{\text{c}}/\left(\ell {\varphi }^{+}\right)}.\end{array}

This gives an explicit formula for the evolution of the damage field. It is obvious that 0 ≤ s < 1 and s satisfies:

(44)\underset{{\varphi }^{+}\to 0}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{1}{1+G_{\text{c}}/\left(\ell {\varphi }^{+}\right)}=0, \underset{{\varphi }^{+}\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{1}{1+G_{\text{c}}/\left(\ell {\varphi }^{+}\right)}=1.

Inserting Eq. (43) back into function g(s), we arrive at

(45)g(s)=\frac{1}{{\left(1+\ell {\varphi }^{+}/G_{\text{c}}\right)}^2}.

Hence, the current model incorporates the dependence of both g(s) and s on the strain energy density that is classified as “positive”. The determination of these quantities can be conducted at a local level by the analysis of strain energy density. To ensure the irreversibility, the positive strain energy density is substituted with the historical maximum positive energy density

(46)H=\underset{\hat{t}\le t}{\mathrm{m}\mathrm{a}\mathrm{x}}{\varphi }^{+}(\text{ε}(\hat{t})).

The governing equations for solid domain become

(47)\mathrm{\nabla }\cdot \left(\frac{1}{{\left(1+\ell H/G_{\text{c}}\right)}^2}\frac{\mathrm{\partial }{\varphi }^{+}}{\mathrm{\partial }\text{ε}}+\frac{\mathrm{\partial }{\varphi }^{-}}{\mathrm{\partial }\text{ε}}\right)+\mathit{g}= \mathbf{0},

or

(48)\mathrm{\nabla }\cdot \left(\frac{{\sigma }^{+}}{{\left(1+\ell H/G_{\text{c}}\right)}^2}+{\sigma }^{-}\right)+\mathit{g}= \mathbf{0}.

For other form of g(s), for example g(s)= (1-s{)}^2(1-(m-2 )s) [38], the derivative can be written as g^{\mathrm{\prime }}(s)=(1-s )((3m- 6)s- m). In this case, the damage field equation becomes (1-s)( (3m- 6)s- m){\varphi }^{+}+2s{G}_{\text{c}} /\ell = 0, which is a quadratic equation. The solution is

(49)s=\frac{G_{\text{c}}/\left(\ell {\varphi }^{+}\right)+2m-3\pm \left|G_{\text{c}}/\left(\ell {\varphi }^{+}\right)+2m-3\right|}{3(m-2)}.

When m=2, it recovers s={\varphi }^{+}/\left({\varphi }^{+}+G_{\text{c}}/\ell \right) . It can be observed that the damage model for fracturing is contingent solely upon the mechanics field, whereas the damage field can be locally established by the assessment of energy density. The utilization of this approach presents significant benefits in terms of numerical implementation. The conventional finite element approach, with minor adjustments, can be utilized to numerically apply the current VDM.

By utilizing the equations s^2=1/{\left(1+G_{\text{c}}/\ell {\varphi }^{+}\right)}^2 and g(s)=1/{\left(1+\ell {\varphi }^{+}/G_{\text{c}}\right)}^2, the functional expression can be reduced to

(50)\begin{array}{rl}& \text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{1}{{\left(1+\ell {\varphi }^{+}/G_{\text{c}}\right)}^2}{\varphi }^{+}+{\varphi }^{-}+\frac{1}{{\left(1+G_{\text{c}}/\left(\ell {\varphi }^{+}\right)\right)}^2}\frac{G_{\text{c}}}{\ell }-\mathit{g}\cdot \mathit{u}\text{d}V\\ & \to \text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{{\varphi }^{+}}{1+\ell {\varphi }^{+}/G_{\text{c}}}+{\varphi }^{-}-\mathit{g}\cdot \mathit{u}\text{d}V.\end{array}

By doing so, a super concise functional of fracture mechanics is achieved for the first time. Interesting, the functional density is coincidentally related to the damage field in Eq. (43) by

(51)\frac{{\varphi }^{+}}{1+\ell {\varphi }^{+}/G_{\text{c}}}=\frac{G_{\text{c}}}{\ell }s =\frac{G_{\text{c}}}{\ell }\left(\frac{{\varphi }^{+}}{{\varphi }^{+}+G_{\text{c}}/\ell }\right).

The strain energy density is positive and the limit of the effective strain energy density of a fully crack solid is

(52)\underset{{\varphi }^{+}\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{{\varphi }^{+}}{1+\ell {\varphi }^{+}/G_{\text{c}}}=\frac{G_{\text{c}}}{\ell }.

Equation (52) shows that the functional is consistent with the Griffith principle for fracture for the full-damaged region.

The variation of Φ is

(53)\begin{array}{rl}\delta \text{Φ}=& {{\displaystyle \int }}_{\text{Ω}}\frac{\delta {\varphi }^{+}}{{\left(1+\ell {\varphi }^{+}/G_c\right)}^2}+\delta {\varphi }^{-}-\mathit{g}\cdot \delta \mathit{u}\text{d}V\\ =& {{\displaystyle \int }}_{\text{Ω}}\frac{1}{{\left(1+\ell {\varphi }^{+}/G_{\text{c}}\right)}^2}{\text{σ}}^{+}:\delta \text{ε}+{\text{σ}}^{-}:\delta \text{ε}-\mathit{g}\cdot \delta \mathit{u}\text{d}V\\ =& {{\displaystyle \int }}_{\text{Ω}}-\mathrm{\nabla }\cdot \left(\frac{1}{{\left(1+\ell {\varphi }^{+}/G_{\text{c}}\right)}^2}{\text{σ}}^{+}+{\text{σ}}^{-}\right)\cdot \delta \mathit{u}-\mathit{g}\cdot \delta \mathit{u}\text{d}V\\ & +{{\displaystyle \int }}_{\mathrm{\partial }{\text{Ω}}_{\text{t}}}\underset{\text{traction load:}{t}_{\text{N}}}{\underbrace{\left(\frac{1}{{\left(1+\ell {\varphi }^{+}/G_{\text{c}}\right)}^2}{\text{σ}}^{+}+{\text{σ}}^{-}\right)\cdot \mathit{n}}}\cdot \delta \mathit{u}\text{d}S.\end{array}

For any \delta \mathit{u}, \delta \mathrm{\Phi }= 0 yields the governing equations

(54)\mathrm{\nabla }\cdot \left(\frac{1}{{\left(1+\ell {\varphi }^{+}/G_{\text{c}}\right)}^2}{\text{σ}}^{+}+{\text{σ}}^{-}\right)+\mathit{g}= \mathbf{0},

which is almost the same as Eq. (48).

When considering the irreversibility of fracture evolution and using s^2= 1/ {\left(1+G_{\text{c}}/(\ell H)\right)}^2 and g(s)=1/(1+\ell H/G_{\text{c}}{)}^2, with some manipulation, we have

(55)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{{\varphi }^{+}}{{\left(1+\ell H/G_{\text{c}}\right)}^2}+{\varphi }^{-}-\mathit{g} \cdot \mathit{u}+ \frac{\ell H^2/G_{\text{c}}}{{\underset{\text{a constant}}{\underbrace{\left(1+\ell H/G_{\text{c}}\right)}}}^2} \text{d}V.

Since the last term is viewed as a constant, the energy functional can be simplified as

(56)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{{\varphi }^{+}}{{\left(1+\ell H/G_{\text{c}}\right)}^2}+{\varphi }^{-}-\mathit{g} \cdot \mathit{u}\text{d} V.

The limit of the constant is

(57)\underset{H\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\ell H^2/G_{\text{c}}}{{\left(1+\ell H/G_{\text{c}}\right)}^2}=\frac{G_{\text{c}}}{\ell }.

Meanwhile, the strain energy density with crack in Eq. (56) has the limit

(58)\underset{H\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\underset{{\varphi }^{+}\to H^{-}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{{\varphi }^{+}}{{\left(1+\ell H/G_{\text{c}}\right)}^2}=0.

Here the calculation of limit is from right to left and {\varphi }^{+}\to H^{-} denotes that ϕ⁺ approaches H from the left side. We can see that Eq. (58) is consistent with Eq. (52), where the effective energy reduces to zero and crack surface energy converges to a constant. By introducing the historical maximal energy density, the functional is a standard quadratic form and the variational derivation is the same as conventional elasticity.

6 Higher order variational damage model

In Section 2, it is seen that the numerical value “2” in Eq. (13) does not affect the integration process for the distribution of the damage field depicted in Fig.1(b). In this section, we substitute the variable “2” with an arbitrary positive numerical value in Eq. (13), and subsequently introduce a higher order VDM.

6.1 Higher order surface density functions

If we choose the function g(s)=( 1-s{)}^{m+1} and \mathrm{\Gamma } (s)=s^{m+1} with real number m>0, the functional becomes

(59)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}{(1-s)}^{m+1}\varphi +s^{m+1}\frac{G_{\text{c}}}{\ell }\text{d}V.

Here we do not distinguish the “negative”/“positive” parts of the strain energy density. For different \text{m}, (1-s{)}^{m+1} and s^{m+1} are plotted in Fig.2. One can see that these curves are symmetric with respect to the line s = 0.5. The larger m is, the steeper the curve is. The first variation of Φ for any δs yields

(60)s=\frac{{\left(\ell \varphi /G_{\text{c}}\right)}^{1/m}}{1+{\left(\ell \varphi /G_{\text{c}}\right)}^{1/m}}.

The crack energy density becomes

(61)\frac{G_{\text{c}}}{\ell }\mathrm{\Gamma }(s)=\frac{G_{\text{c}}}{\ell }{\left(\frac{{\left(\ell \varphi /G_{\text{c}}\right)}^{1/m}}{1+{\left(\ell \varphi /G_{\text{c}}\right)}^{1/m}}\right)}^{m+1}.

Inserting the value of s back into Φ, Φ can be simplified as

(62)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{\varphi }{{\left(1+{\left(\ell \varphi /G_{\text{c}}\right)}^{1/m}\right)}^m}\text{d}V.

Hence, a nonlinear energy functional with varying orders for fractures is achieved, derived from the concepts of energy density and surface density orders. These surface density orders can be tailored to model materials with unconventional fracture behaviors. The formulation of this function employs a nonlinear polynomial, as opposed to a regular polynomial. The functional incorporates the crucial parameters \left(G_{\text{c}},\ell \right) in a creative manner.

If we set G_{\text{c}} /\ell =1, the function \varphi /{\left(1+{\varphi }^{1/m}\right)}^m for different m is plotted in Fig.3. One can see that the smaller m, the faster conversion of strain energy density into fracture surface energy density. For larger m, the fracture requires more energy to develop. However, for the case of m = 1, the function is most concise and the conversion velocity is moderate.

The series expansions of 1/ (1+x ) and 1/(1+x{)}^m are

(63)\begin{array}{rl}\frac{1}{1+x}=& \sum _{i=0}^{\mathrm{\infty }}(-x{)}^t=1-x+x^2-x^3+x^4-\cdots ,\\ \frac{1}{(1+x{)}^m}=& \sum _{i=0}^{+\mathrm{\infty }}{C}_{m+i-1}^t(-x{)}^i=1-mx+\frac{m(m+1)}{2!}{x}^2\\ & -\frac{m(m+1)(m+2)}{3!}{x}^3+\cdots ,\end{array}

where C_{m+i-1}^i= {\displaystyle \frac{1}{i!}} {\mathrm{\Pi }}_{j=0}^{t-1}(m+j), C_{m-1}^0 =1. These two series are valid in .

Based on the series expansion above, the Taylor series of {\displaystyle \frac{\varphi }{{\left(1+{\left(\ell \varphi /G_{\text{c}}\right)}^{1/m}\right)}^m}} is

(64)\begin{array}{rl}& \frac{\varphi }{1+\ell \varphi /G_{\text{c}}}=\varphi \sum _{i=0}^{+\mathrm{\infty }}{\left(-\frac{\ell \varphi }{G_{\text{c}}}\right)}^i=\varphi \left(1-\frac{\ell \varphi }{G_{\text{c}}}+{\left(\frac{\ell \varphi }{G_{\text{c}}}\right)}^2-\cdots \right),m=1,\\ & \frac{\varphi }{{\left(1+{\left(\ell \varphi /G_{\text{c}}\right)}^{1/m}\right)}^m}=\varphi \sum _{i=0}^{+\mathrm{\infty }}{(-1)}^i{C}_{m+i-1}^i{\left(\frac{\ell \varphi }{G_{\text{c}}}\right)}^{\frac{i}{m}},m>0.\end{array}

It is evident that the aforementioned energy functional exhibits nonlinearity and encompasses an endless number of higher order polynomial terms pertaining to the energy density. The series expression converges when the condition is satisfied.

The first variation of Φ becomes

(65)\delta \text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{\delta \varphi }{{\left(1+{\left(\ell \varphi /G_{\text{c}}\right)}^{1/m}\right)}^{m+1}}\text{d} V.

When the “positive”/“negative” energy and the maximal history “positive” energy density are used, the fracture energy functional is

(66)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{{\varphi }^{+}}{{\left(1+{\left(\ell H/G_{\text{c}}\right)}^{1/m}\right)}^{m+1}}+{\varphi }^{-}\text{d}V

The governing equation of solid becomes

(67)\mathrm{\nabla }\cdot \left(\frac{{\text{σ}}^{+}}{{\left(1+{\left(\ell H/G_{\text{c}}\right)}^{1/m}\right)}^{m+1}}+{\text{σ}}^{-}\right)+\mathit{g}= \mathbf{0}

where {\sigma }^{\pm }={\displaystyle \frac{\mathrm{\partial }{\varphi }^{\pm }}{\mathrm{\partial }\epsilon }}.

In the present model, the fracture length scale ℓ is a material parameter that is primarily associated with the theoretical framework rather than the practical numerical implementation.

6.2 Variational damage model in one dimension

For the purpose of illustration, we consider the dimensionless energy density of a bar in 1D as

(68)\tilde{\varphi }=\frac{\varphi }{1+\varphi },\text{with}\varphi =\frac{1}{2}{\epsilon }^2,\epsilon =\frac{\text{d}u}{\text{d}x},

where ε represents the strain, u is the displacement, fracture parameter \ell /G_{\text{c}} =1 and the stiffness coefficient k = 1 are assumed. The evolution of \tilde{\varphi } and its gradient with respect to strain energy density and strain are depicted in Fig.4. From Fig.4(a), it is evident that there is a positive relationship between the variables \tilde{\varphi } and ϕ. However, the rate of increase diminishes rapidly. When the strain ε is selected as the independent variable, Fig.4(b) illustrates that the effective stress exhibits an increasing trend for small values of ε, but falls after ε beyond a specific threshold value. There is a direct relationship between the reduction of effective stress and the expansion of the damage field.

In the 1D case, the effective stress is

(69)\sigma =\frac{\mathrm{\partial }\tilde{\varphi }}{\mathrm{\partial }\epsilon }=\frac{E\epsilon }{{\left(1+{\left({\displaystyle \frac{1}{2}}E{\epsilon }^2\ell /G_{\text{c}}\right)}^{1/m}\right)}^{m+1}} .

With the aid of Mathematica [39], one can verify that for any order m, there is

(70){\int }_0^{+\mathrm{\infty }}\sigma \text{d}\epsilon =\frac{G_{\text{c}}}{\ell }.

In the case of E=1,\ell /G_{\text{c}} =1, we plot σ for different values of m in Fig.5, which shows that the larger m, the flatter the σ curve is. Fig.5(a) presents the curves for m < 1. It is interesting to note that as m approaches zero, the curve increases linearly before abruptly dropping to zero when the peak value reaches \sqrt{2}. This behavior indicates a full-brittle damage at that point. In this extreme scenario, the curve comprises three line segments, resembling to the bond-cutting process in bond-based PD, where the internal force is immediately removed upon exceeding the critical strain [40,41]. In addition, the case of m = 0 corresponds exactly to the degradation function used in the damage model in Eq. (4). This might be the reason why conventional damage model is tedious in implementation, while m = 0 is usually “neglected” in the current VDM so as to achieve the most efficient implementation of calculating the damage field. Fig.5(b) illustrates the effective stress curves for m ∈ {1,2,3,4,5,6}. It is evident that the peak stress value is reduced approximately by half for each unit increase in m.

In Fig.4(b), we can see that the effective stress increases monotonically and then decreases monotonically after reaching a maximum of s=\sqrt{2/3}. For the general case in 1D, the maximal stress {\sigma }_m corresponds to the critical strain {\epsilon }_{\text{c}}satisfying {\displaystyle \frac{\text{d}\sigma }{\text{d}\epsilon }}=0 which gives

(71){\epsilon }_{\text{c}} =\sqrt{\frac{2G_{\text{c}}}{E\ell }{\left(\frac{m}{m+2}\right)}^m}.

The critical stress, critical strain energy density and the damage can accordingly be formulated as

(72)\begin{array}{rl}& {\sigma }_{\text{max}}=\sqrt{\frac{2E{G}_{\text{c}}}{\ell }}\frac{{\left({\displaystyle \frac{m}{m+2}}\right)}^{m/2}}{{\left(1+{\displaystyle \frac{m}{m+2}}\right)}^{m+1}},\\ & {\varphi }_{\text{c}}=\frac{1}{2}E{\epsilon }_{\text{c}}^2=\frac{G_{\text{c}}}{\ell }{\left(\frac{m}{m+2}\right)}^m,\\ & s_{\text{c}}=\frac{m}{2(m+1)}.\end{array}

The maximum stress is equivalent to the critical load of the 1D structure. The plot of the maximum stress as a function of the order of damage model is depicted in Fig.6. The data clearly indicates a decrease in the maximal effective stress as the parameter m increases.

7 Further generalization

7.1 Generalization to other physical fields

The general framework is based on the fundamental modification of the energy functional. The degradation of the energy density is inherently or implicitly described by the energy density. For the simplest case, the modified fracture energy functional are written as

(73)\begin{array}{rl}& \text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{{\varphi }^{+}}{1+\ell {\varphi }^{+}/G_{\text{c}}}+{\varphi }^{-}-\mathit{g}\cdot \mathit{u}\text{d}V,\\ & \text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{{\varphi }^{+}}{{\left(1+\ell H/G_{\text{c}}\right)}^2}+{\varphi }^{-}-\mathit{g}\cdot \mathit{u}\text{d}V.\end{array}

The damage model described above is founded upon the concept of energy density. The utilization of the equation \varphi = 1/ 2\sigma :\epsilon allows for the natural modeling of fractures in solid materials. The energy density is not solely confined to the elastic mechanics of solids. There exist many functionals that bear resemblance to elasticity mechanics, for instance, the bending energy density in plate/shell theories, the energy density in gradient elasticity. Hence, the incorporation of damage evolution in other theories such as gradient elasticity and plate/shell can be seamlessly integrated within the existing framework.

We take the thin plate theory for an example

(74)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{\varphi (\text{κ})}{1+\ell \varphi (\text{κ})/G_{\text{c}}}+pw\text{d} V,

where bending energy density \varphi (\text{κ})=1/2\mathit{m}:\kappa, curvature tensor \text{κ}=\mathrm{\nabla }\mathrm{\nabla } w, w is the deflection and p is the transverse load on the mid-surface of the plate. ℓ and G_{\mathrm{c}} are material parameters related to crack in plate. The bending moment is expressed as \mathit{m}=\mathbb{D}:\kappa ={\displaystyle \frac{\mathrm{\partial }\varphi (\text{κ})}{\mathrm{\partial }\text{κ}}}, \mathbb{D} is the material tensor of the thin plate.

The variation of Φ yields

(75)\begin{array}{rl}\delta \text{Φ}=& {{\displaystyle \int }}_{\text{Ω}}\frac{\delta \varphi (\text{κ})}{{\left(1+\ell \varphi (\text{κ})/G_{\text{c}}\right)}^2}+p\delta w\text{d}V\\ =& {{\displaystyle \int }}_{\text{Ω}}\frac{\mathit{m}:\delta \kappa }{{\left(1+\ell \varphi (\text{κ})/G_{\text{c}}\right)}^2}+p\delta w\text{d}V\\ =& {{\displaystyle \int }}_{\text{Ω}}\mathrm{\nabla }\cdot \mathrm{\nabla }\left(\frac{\mathit{m}}{{\left(1+\ell \varphi (\text{κ})/G_{\text{c}}\right)}^2}\right)\delta w\\ & +p\delta w\text{d}V+{{\displaystyle \int }}_{S}BC\text{terms}.\end{array}

For any δw, δΦ = 0 gives the governing equations of thin plate embedded with VDM

(76)\mathrm{\nabla }\cdot \mathrm{\nabla }\left(\frac{m}{{\left(1+\ell \varphi (\text{κ})/G_{\text{c}}\right)}^2}\right)+p=0.

The damage variable of the plate at any point is given by

(77)s=\frac{\varphi (\text{κ})}{\varphi (\text{κ})+G_{\text{c}}/\ell }.

For the case of gradient elasticity, the gradient energy density is \varphi (\text{ε},\mathit{\chi })=1 /2 (\text{σ}:\text{ε}+\mathrm{\Sigma }:\mathit{\chi } ), where \mathit{\chi }=\mathrm{\nabla } \text{ε} and gradient stress tensor \mathbf{\Sigma }=\mathbb{E}:\mathit{\chi } with \mathbb{E} being the 6-th order material tensor. Using the total energy for the evolution of crack, the energy functional with variational damage can be written as

(78)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{\varphi (\text{ε},\mathit{\chi })}{1+\ell \varphi (\text{ε},\mathit{\chi })/G_{\text{c}}}+\mathit{g}\cdot \mathit{u}\text{d} V.

The first variation of Φ reads

(79)\begin{array}{rl}\delta \text{Φ}=& {{\displaystyle \int }}_{\text{Ω}}\frac{\text{σ}:\delta \text{ε}+\mathrm{\Sigma }:\delta \mathit{\chi }}{{\left(1+\ell \phi (\text{ε},\mathit{\chi })/G_{\text{c}}\right)}^2}-\mathit{g}\cdot \delta \mathit{u}\text{d}V\\ =& {{\displaystyle \int }}_{\text{Ω}}\mathrm{\nabla }\cdot \mathrm{\nabla }\cdot \left(\frac{\mathrm{\Sigma }}{{(1+\ell (ϵ,\mathit{\chi })/G_{\text{c}})}^2}\right)\cdot \delta \mathit{u}-\mathit{g}\cdot \delta \mathit{u}\text{d}V\\ & -\mathrm{\nabla }\cdot \left(\frac{\text{σ}}{{(1+\ell (\text{ε},\mathit{\chi })/G_{\text{c}})}^2}\right)\cdot \delta \mathit{u}\text{d}V+{{\displaystyle \int }}_S\text{BC terms}.\end{array}

For any \delta \mathit{u}, δΦ = 0 leads to

(80)\mathrm{\nabla }\cdot \left(\frac{\text{σ}}{{\left(1+\ell \varphi (\text{ε},\mathit{\chi })/G_{\text{c}}\right)}^2}\right)-\mathrm{\nabla }\cdot \mathrm{\nabla }\cdot \left(\frac{\mathbf{\Sigma }}{{\left(1+\ell \varphi (\text{ε},\mathit{\chi })/G_{\text{c}}\right)}^2}\right)+\mathit{g}= \mathbf{0}.

From above examples, we can see that the VDMs are very versatile in accounting for different fracture problems.

7.2 Variational multi-damage model

There are many forms of energies for multi-physics problems. For example, the gradient elasticity theory, the functional is

(81)\mathrm{\Phi }={\int }_{\mathrm{\Omega }}\underset{{\varphi }_1(\text{ε})}{\underbrace{\frac{1}{2}\text{σ}:\text{ε}}}+\underset{{\varphi }_2(\mathit{\chi })}{\underbrace{\frac{1}{2}\mathbf{\Sigma }:\mathit{\chi }}}\text{d}V .

The degradation of the energy of gradient strain and the strain energy could be different, which implies that two different damage states might exist. For this end, we devise a new energy functional based on the VDM as

(82)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{{\varphi }_1}{1+{\ell }_1{\varphi }_1/G_1}+\frac{{\varphi }_2}{1+{\ell }_2{\varphi }_2/G_2}\text{d}V,

where G₁ and G₂ are the energy release rate for solid structure and gradient solid structure, respectively; ℓ₁ and ℓ₂ are the crack length scale for solid structure and gradient solid structure, respectively. Then the governing equations become

(83)\mathrm{\nabla }\cdot \left(\frac{\text{σ}}{{\left(1+{\ell }_1{\varphi }_1(\text{ε})/G_1\right)}^2}\right)-\mathrm{\nabla } \cdot \mathrm{\nabla }\cdot \left(\frac{\mathbf{\Sigma }}{{\left(1+{\ell }_2{\varphi }_2(\mathit{\chi })/G_2\right)}^2}\right) +\mathit{g}=\mathbf{0}.

Equation (80) is similar to Eq. (83) except that the latter uses multiple damage models.

More generally, for different kinds of energies, the new functional can be written as

(84)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\sum _{i=1}^N\frac{{\varphi }_i}{1+{\ell }_i{\varphi }_i/G_{\text{c}i}}\text{d}V,

which can be easily derived from functional Φ = {{\displaystyle \int }}_{\text{Ω}} {\displaystyle \sum _{i=1}^N{\left(1-s_i\right)}^2{\varphi }_i+{\displaystyle \frac{G_{\text{c}i}}{{\ell }_i}}{s}_i^2\text{d}V} by calculating the first variation, where G_{\text{c}i} and {\ell }_1 are the fracture parameters related to energy density {\varphi }_1.

We take the elastic solid for an example. In the elastic solid, the strain energy density can be decomposed into deviatoric part and volumetric part, which indicates that degradation of the deviatoric deformation and the degradation of the volumetric part might be enforced separately as:

(85)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{{\varphi }_{\text{d}}}{1+{\ell }_{\text{d}}{\varphi }_{\text{d}}/G_{\text{c}}^{\text{d}}}+\frac{{\varphi }_{\text{v}}}{1+{\ell }_{\text{v}}{\varphi }_{\text{v}}/G_{\text{c}}^{\text{v}}}\text{d}V,

where \left({\ell }_{\text{d}},{\ell }_{\text{v}},G_{\text{c}}^{\text{d}},G_{\text{c}}^{\text{v}}\right) are the crack parameters associated to the energy densities.

The generalization to space with different dimensions is straightforward. We consider the combination of different energies in different dimensions. The fracture energy functional of an elastic solid reinforced with N bars can be written as

(86)\text{Φ}={{\displaystyle \int }}_{\text{Ω}}\frac{{\varphi }_1}{1+\ell {\varphi }_1/G_{\text{c}}}\text{d}V+\sum _{i=1}^N{{\displaystyle \int }}_{L_i}\frac{{\psi }_i}{1+{\ell }_i{\psi }_i/G_i}\text{d}L,

where {\varphi }_1 is the energy density of elastic solid and {\psi }_t is the energy density of the bar per unit length, L_i is the length of ith bar. The reinforced effect of the bar on the solid domain can be considered by attaching the stiffness of the bar to the respected degree of freedom in the solid domain. The damages in elastic solid and bars are considered directly.

The single-field framework of VDM can be applied to multi-physical problems such as modeling batteries and solar cells [42], which involve displacement, electric, chemical, and thermal fields. These fields, as fundamental aspects of the problem, can induce various types of damage or material degradation, not just fractures in solids. Therefore, it is possible to use VDM to model discontinuities induced by these fields without introducing an additional field for the discontinuity, simplifying the complexity of multi-physical problems.

8 Finite element implementation

The VDMs proposed in this study are classified as theoretical models for describing fracture phenomena. They are independent of the numerical methods. Various numerical methods such as FEM [43–45], isogeometric analysis [46,47], meshless method [7,48,49], PD [17,50], virtual element method [51], nonlocal operator methods [19], physics-informed neural networks [52–54] can be readily utilized for numerical implementation. To streamline the analysis, we employ conventional finite element techniques to investigate the behavior of cracks within the field of elasticity mechanics.

8.1 Shape function of finite element method

Given that the energy density fully determines the damage field at a local level, it becomes necessary to discretize solely the displacement field. The displacement within a quadrilateral element consisting of four nodes is interpolated as:

(87){\mathit{u}}^{\text{h}}(\mathit{x})=\mathit{N}(\mathit{x}){\left({\mathit{u}}_1^{\text{T}},{\mathit{u}}_2^{\text{T}},{\mathit{u}}_3^{\text{T}},{\mathit{u}}_4^{\text{T}}\right)}^{\text{T}},

where \mathit{N}(\mathit{x}) is the vector of shape functions, {\left({\mathit{u}}_1^{\text{T}},{\mathit{u}}_2^{\text{T}},{\mathit{u}}_3^{\text{T}},{\mathit{u}}_4^{\text{T}}\right)}^{\text{T}} is the nodal displacement vector. The two dimensions (2D) strain tensor in vector-form can be formulated as

(88)\left(\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\epsilon }_{xy}\end{array}\right)=\underset{{\mathit{B}}_{\text{e}}}{\underbrace{\left(\begin{array}{llll}{\hat{\mathit{N}}}_1& {\hat{\mathit{N}}}_2& {\hat{\mathit{N}}}_3& {\hat{\mathit{N}}}_4\end{array}\right)}}{\left({\mathit{u}}_1^{\text{T}},{\mathit{u}}_2^{\text{T}},{\mathit{u}}_3^{\text{T}},{\mathit{u}}_4^{\text{T}}\right)}^{\text{T}},

where the term {\hat{\mathit{N}}}_i in strain matrix {\mathit{B}}_{\text{e}} and nodal unknowns are expressed as

(89){\hat{\mathit{N}}}_i=\left(\begin{array}{cc}{\mathrm{\partial }}_x{N}_i& 0\\ 0& {\mathrm{\partial }}_y{N}_i\\ {\mathrm{\partial }}_y{N}_i& {\mathrm{\partial }}_x{N}_i\end{array}\right),{\mathit{u}}_i=\left(\begin{array}{c}{u}_i\\ v_i\end{array}\right).

For simplicity, we overlook the detailed FEM formulation in this paper. The interested reader can resort to standard textbook of FEM. The elements used in the numerical examples include triangular element (denoted by T1), quadrilateral element (denoted by Q1), tetrahedral element (denoted by O1), hexahedral element (denoted by H1).

8.2 Spectral decomposition of strain tensor

We employ the spectral decomposition of strain tensor to calculate the “positive” and “negative” parts of the strain energy, where the “positive” energy drives the crack propagation, while stresses from the “negative” energy density remain intact. For a strain tensor with eigenvalue decomposition (e.g. \text{ε}={\displaystyle \sum _{j=1}^3{\text{ε}}_j{\mathit{n}}_j\otimes {\mathit{n}}_j}), the stress tensor in Eq. (31) can be calculated with ease

(90)\begin{array}{rl}{\text{σ}}^{\pm }=& \lambda {⟨{\epsilon }_1+{\epsilon }_2+{\epsilon }_3⟩}_{\pm }\mathit{I}\\ & +2\mu \left({⟨{\epsilon }_1⟩}_{\pm }{\mathit{n}}_1\otimes {\mathit{n}}_1+{⟨{\epsilon }_2⟩}_{\pm }{\mathit{n}}_2\otimes {\mathit{n}}_2+{⟨{\epsilon }_3⟩}_{\pm }{\mathit{n}}_3\otimes {\mathit{n}}_3\right).\end{array}

We assume that these three eigenvalues are distinct. For repeated eigenvalues, a small perturbation can be applied to make them distinct. Since an implicit implementation is used, the derivative of the stress tensor with respect to the strain tensor must be calculated. The variables in {\text{σ}}^{\pm } include {\epsilon }_f and {\mathit{n}}_i for i\in \{1, 2,3\}. Therefore, the derivatives of eigenvalues and eigenvectors with respect to the strain tensor ε should be computed, which can be expressed as

(91)\mathrm{\partial }{\epsilon }_i= {\mathit{n}}_i^{\text{T}}\mathrm{\partial }(\text{ε}){\mathit{n}}_i,\mathrm{\partial }{\mathit{n}}_i={\left({\epsilon }_i\mathbf{I}-\text{ε}\right)}^{-1}\mathrm{\partial }( \text{ε}){\mathit{n}}_i.

The generalized inverse of \left({\epsilon }_i\mathbf{I}-\text{ε}\right) has the explicit form: {\left({\epsilon }_i\mathbf{I}-\text{ε}\right)}^{-1}={\displaystyle \sum _{j\ne i}^3{\displaystyle \frac{1}{{\epsilon }_i-{\epsilon }_j}}{\mathit{n}}_j\otimes {\mathit{n}}_j}.

Then the matrix derivative can be written as well as

(92)\frac{\mathrm{\partial }{\epsilon }_i}{\mathrm{\partial }\text{ε}}= {\mathit{n}}_i\otimes {\mathit{n}}_i, \frac{\mathrm{\partial }{\mathit{n}}_i}{\mathrm{\partial }\text{ε}}= \sum _{j\ne i}^3\frac{1}{{\epsilon }_i-{\epsilon }_j}{\mathit{n}}_j\otimes {\mathit{n}}_j\otimes {\mathit{n}}_i.

Based on Eq. (92), an explicit expression of the fourth-order material tensor {\mathbb{D}}^{\pm }={\displaystyle \frac{\mathrm{\partial }{\text{σ}}^{\pm }}{\mathrm{\partial }\text{ε}}} can be derived by hand or symbolic mathematical software such as Mathematica [39]. The Voigt notation can be employed to convert {\mathbb{D}}^{+} into matrix form of material tensor, denoted by {\mathit{D}}^{\pm }. Consequently, the tangent stiffness matrix of one element becomes

(93){\mathit{K}}_{\text{e}}={\int }_{{\mathrm{\Omega }}_{\text{e}}} {\mathit{B}}_{\text{e}}^{\text{T}}\left(\frac{{\mathit{D}}^{+}}{{\left(1+\ell H/G_{\text{c}}\right)}^2}+{\mathit{D}}^{-}\right){\mathit{B}}_{\text{e}} \text{d}V,

where {\mathit{K}}_{\text{e}} denotes the elemental tangent stiffness matrix and {\mathrm{\Omega }}_{\text{e}} is the region of one element. For a more detailed implementation of spectral decomposition, the reader is referred to Ref. [55].

8.3 Enforcement of initial crack/discontinuity

In VDM, the maximal historical “positive” energy density determines the distribution of the crack field. There are two approaches to express the initial crack or discontinuity in solid. The first one is to model the crack directly in the meshing process. Another approach is to modify the historical “positive” energy density without meshing the crack in the geometrical model. One can assign H_i= {10}^3{G}_i/4, where the corresponding damage value at Gauss point i is approximately s_i=H_i/(H_i+G_i/\ell )={10}^5/\left({10}^5+1\right)\approx 1.

8.4 Adaptive time increment

During the progression of the damage field, the distribution of strain energy is modified at every increment. The rate of damage progression accelerates significantly when the structural integrity of a system is compromised at a subsequent phase. To effectively represent the damage, we propose using of an adaptive load increment

(94)\mathrm{\Delta }t= \mathrm{\Delta }{t}_0\times max\left(\frac{1}{{\left(1+\ell H_{max}/G_{\text{c}}\right)}^2},{10}^{-3}\right),

where \mathrm{\Delta }{t}_0 represents the initial load increment and H_{max} is defined as the maximum value of historical energy field H in the computational domain. In the initial stages, where the magnitude of strain energy is relatively low, it is possible to utilize a significantly larger load increment, hence resulting in computational time savings. During the later stage, the load increment is constrained to tiny magnitudes in order to effectively monitor the rapid progression of the damage field.

9 Numerical examples

In this section, we present a series of numerical tests designed to evaluate the capabilities of the VDM, as implemented in the previous section. These tests include: a single-edge-notched plate subjected to tensile/shear displacement boundary conditions in both 2D and 3D, an L-shaped panel test, a symmetric three-point bending test and a non-planar crack growth. Each example aims to demonstrate the effectiveness and robustness of the VDM under various loading conditions and geometrical configurations. All numerical examples are based on Eq. (48).

9.1 Single-edge-notched tension test in two dimensions

In this subsection, we model the single-edge-notched tension test, which is a squared plate with initial notched crack as shown in Fig.7. The material parameters are set as λ = 121.1538 kN/mm² and µ = 80.7692 kN/mm² for elastic constants, G_c = 2.7 × 10⁻³ kN/mm for the critical energy release rate. These parameters are identical to that used in the small strain brittle fracture phase field in Ref. [13]. Two displacement conditions are tested: Case a) for tensile boundary condition and Case b) for shear boundary conditions.

9.1.1 Variational damage model vs. phase field method

For the same mesh, as shown in Fig.8, the crack field in the phase field is diffusive while that in the VDM is much more concentrated. Besides solving a sparse array of phase fields, the width of the crack field occupies the length of 4–6 elements. VDM can obtain a sharper fracture interface with the same discretization or equivalent sharp fracture interface with coarse mesh. In the phase field method, the phase field can be only obtained after assembling the matrix and solving the linear algebra equations. This additional work is not required in VDMs. We compared the computational times of the phase field method and VDM. In a laptop with i5-1135 G7 CPU, the computational time is 0.318 s/step, 0.225 s/step for phase field model and VDM, respectively, where the latter uses 70.8% of the time of the former. Hence, in a simpler theoretical framework, VDM is more efficient in computation as well as in numerical implementation.

9.1.2 Influence of length scale ℓ

In the VDM, ℓ is viewed as a parameter related to the length scale. Large \ell indicates the smaller energy density that induces damage, and smaller \ell means larger energy density or “stronger” material. With the same discretization of 60 × 60 elements, we alter the value of \ell and calculate the crack pattern and critical load. In the following cases, \ell = \{0.4,0.6,0.7, 0.8,1,2 ,3\} \mathrm{\Delta }x, where \mathrm{\Delta } x is the size of the elements. The load curves for different \ell are depicted in Fig.9. It can be seen that larger \ell corresponds to smaller critical load and smaller \ell conforms to larger critical load. We extracted the maximal values of these load curves and plotted them on Fig.10. Interesting, we found that the critical load is proportional to the inverse square root of \ell. The relation can be described by the fitting curve F_{\text{max}}( \ell )=0.068+0.527/\sqrt{\ell /\mathrm{\Delta }x}. It is worth noting that the varying of \elldid not affect the width of the crack interface, as shown in Fig.11. The crack patterns for \ell = 3\mathrm{\Delta } x and \ell = 0.4\mathrm{\Delta }x are almost the same, though the maximal loads are quite different. Based on the above analysis of \ell, we select the value of \ell as \ell = 0.7\mathrm{\Delta } x in the following calculations. We also test the influence of mesh. We employ three different mesh densities: 30 × 30, 60 × 60, and 120 × 120 elements. The crack interfaces are depicted in Fig.12. It can be seen that the finer mesh yields a sharper crack interface. In addition, we plot the load curves in Fig.13, where the mesh density has only a slight influence on load curves.

9.1.3 Shear boundary conditions

The quadrilateral mesh (Q1) and triangle mesh (T1) are employed for discretization of the notched plate under shear-displacement boundary conditions. The study examines four different mesh settings. Both Case 1 and Case 2 utilized a uniform Q1 mesh. In Case 3 and Case 4, the T1 mesh was employed and subsequently refined along the anticipated crack path. The first case has 10⁴ Q1 elements, the second case has 4 × 10⁴ Q1 elements, the third case has 10⁴ T1 elements, and the fourth case has 4 × 10⁴ T1 elements. The effective minimal element sizes in all four cases are given by \mathrm{\Delta } x=\{ 1.0\times {10}^{-2},5\times {10}^{-3},6.9\times {10}^{-4},3.5\times {10}^{-4}\} mm. For all situations, the crack length scale parameter is chosen as ℓ = 0.7∆x. Fig.14 displays the triangle meshes corresponding to Case 4.

The fracture patterns are depicted in Fig.15. In the depicted illustration, it is evident that all observed crack patterns have a consistent path. When the mesh size is reduced, the propagation of cracks is confined to smaller regions. In the context of T1 mesh, the distribution of the mesh has a minor influence on the fracture path, however the overall trajectory of the crack remains accurate. The load curves are plotted in Fig.16(a) and Fig.16(b).

9.1.4 Tensile boundary conditions with refinement

To enhance computational efficiency, the mesh has been refined in regions anticipated to undergo crack propagation. The mesh is depicted in Fig.17, utilizing a three-level refinement technique. The discrete model comprises a total of 10345 nodes and 10010 Q1 elements. The maximum size of the element is denoted as ∆x_max and is equal to 0.05 mm. Meanwhile, the minimum size of the element is ∆x_min = ∆x_max/27. The initial fracture is formed through the process of eliminating specific elements in the region of the initial crack. Without using a mesh refinement, the total number of elements exceeds 290000. The length scale is selected as ℓ = 0.7∆x_min. The material parameters utilized in this example are the same with those employed in prior cases. The depiction of the crack path with a high level of sharpness can be observed in Fig.18. The region near the crack front exhibits a higher deformation gradient in fracture problems, whereas the deformation gradient in other regions is low. Therefore, the implementation of mesh refinement in the crack path can significantly enhance computational efficiency while preserving the high resolution of the crack interface.

9.2 Three dimensional notched plate in tensile boundary

To reduce computational cost, we employ the refinement technique of H1 mesh alongside the crack path. Three-level refinement is used and the model contains 2.44 × 10⁵ nodes and 2.26 × 10⁵ H1 elements, which is shown in Fig.19. Without refinement, the number of H1 mesh would be around 1.96 million. The crack paths at different displacement loads are depicted in Fig.20. The crack initiated at u_y = 5.82 × 10⁻³ mm and propagated faster as the displacement increased. The crack pattern in 3D is very similar to the 2D cases. We can see that the crack paths at different stages fall inside the refine region. Very sharp crack paths are achieved in this model. We also plot the load displacement curve in Fig.21. The maximal load max(F_N) = 0.0365 kN occurs when u_y ≈ 5.9 × 10⁻³ mm, which is approximately 1/20 of the 0.74 kN of 2D case in Fig.13. The reason is that the thickness of the 2D model is 1 mm, while the thickness of the 3D model is selected as 0.05 mm, which is 1/20 of the thickness of the 2D model.

9.3 Three dimensional notched plate subjected to shear boundary

The present study investigates the response of a notched plate subjected to shear boundary conditions in a 3D setting. The material parameters utilized in this study are consistent with those employed in prior cases. Two different discretizations are being considered. In the first case, the plate is discretized into a total of 410000 H1 elements, having a uniform mesh size of ∆x = 1/160 mm. In the second case, the plate is meshed into 3.77 × 10⁵ O1 elements. A dense mesh is employed along the possible fracture path, with a minimum mesh size of ∆x = 1/160 mm. The distribution of mesh for the example is illustrated in Fig.22. The length scale is selected as ℓ = 1.4∆x. The 3D crack patterns of both cases are illustrated in Fig.23. The correspondence between the fracture paths in 3D and 2D scenarios is evident. In Case 2, the distribution of the O1 mesh has a small influence on the fracture path. Nevertheless, the overall trajectory remains consistent.

9.4 Symmetric three point bending test in three dimensions

The example handles a three-point bending test on a simply supported notched beam in a 3D setting. The beam has dimensions of 8 mm in length, 2 mm in width, and 0.5 mm in height. Fig.24 depicts both the geometric arrangement and the mesh distribution. The mesh is refined by three times in the anticipated crack propagation region, resulting in a discretization including 46848 elements and a minimum element size of approximately ∆x = 1.85 × 10⁻³ mm. The bulk modulus and shear modulus are selected as k = 12.00 kN/mm² and µ = 8.0 kN/mm², while the critical energy release rate is denoted as G_{\mathrm{c}} = 5.0 × 10⁻⁴ kN/mm. The crack length scale is selected as ℓ = 0.6∆x. The calculation is conducted by monotonic displacement and utilizes adaptive displacement increments.

Fig.25 illustrates the crack surface in 3D. Fig.26 displays the force−displacement curves. The figure shows that a very sharp crack pattern can be obtained by current VDM.

9.5 L-shaped panel

The unique geometry and potential for stress concentrations at the “L”’s corner make simulating the spread of cracks in L-shaped panels challenging. Numerous computational methods for fractures have been used to analyze the rupture in the L-shaped panel. This subsection examines this problem in two dimensions.

Fig.27(a) illustrates the configuration of the L-shaped panel test, including the loading conditions and the material parameters. Plane stress condition is assumed. The associated experimental experiments are documented in Ref. [56].

Two mesh settings are considered. The first case used 29907 Q1 elements, and case 2 has 101545 Q1 elements. Along the crack paths, the mesh is refined by three-levels. The mesh of case 1 is shown in Fig.27(b). The length scales for two cases are selected as ℓ = 9.24 mm and ℓ = 4.96 mm, respectively.

The crack pattern of case 2 is shown in Fig.28, which agrees well with the phase field method and the experiments. The load–displacement curves are depicted in Fig.29. The peak loads of both cases agree well with the experiments.

9.6 Non-planar crack growth

Non-planar crack development is a more sophisticated case. Gravouil et al. [57] investigated this example using the XFEM. As shown in Fig.30, a beam is notched at a 45° angle with the x-axis and loaded in bending by a moment at the ends. The dimensions of the beam are h = 0.02 m, l = 0.1 m, and d = 0.01 m, and the notch length is a = 0.01 m. Gravouil et al. [57] discovered that when the crack extended downwards, it tended to become normal to the x-axis.

In this example, the material parameters are elastic modulus E = 20 GPa, Poisson ratio ν = 0.2, and G_c = 113 N/m. The damage length scale is selected as ℓ = 1.4∆x, where ∆x = 0.185 mm. The 3D beam is discretized into 390096 H1 elements and 408184 nodes.

In the middle region of the beam, dense mesh is used. With a uniform mesh of mesh size ∆x, the number of elements will reach 3.15 million H1 elements.

The crack pattern is illustrated in Fig.31, and the zoomed crack pattern is shown in Fig.32. One can see the curved crack surface clearly, which matches well with those by the XFEM [57] or by CPM [7].

10 Conclusions

In this paper, we have presented a VDM for the purpose of fracture modeling in structures and materials. The model utilizes the energy density exclusively and represents the fracture energy density as a function of the energy density, while expressing the damage variable as a function of the energy density. The study started from the fundamental modification of energy density functional and devised new fracture energy functional, which forms the basis of VDM of different orders. Meanwhile, we have put forth a multi-damage model that aims to incorporate several types of energy densities, including those associated with thin elastic plates and gradient elasticity.

The model is formulated within a single physical field rather than multi-physical fields, resulting in significant computational cost savings and reduced complexity in numerical implementation. These formulations are applicable in the standard finite element approach as well as in other numerical methods with minimal modification required. The general framework can be readily applied to several physical fields. The single-field formulation and energy density focus in the VDM streamline computational processes, reduce numerical implementation complexity, facilitate efficient calculations, unify governing equations, minimize data handling, and enhance compatibility with various numerical methods. This contributes to cost-effective and straightforward fracture modeling in structures and materials.

In terms of efficiency of fracture modeling, compared with the phase field method of the same mesh density, VDM attains a fracture interface of a higher degree of sharpness. The sharper crack interface offers several advantages, such as improved resolution, accuracy, identification of critical regions, fine-grained fracture analysis, and computational performance, enhancing its utility in simulating complex fracture phenomena and advancing the understanding of fracture behavior in solid materials. We demonstrated the capabilities of the variational models by the utilization of sample 2D and 3D numerical cases.

This damage model is formulated solely based on the energy density, without taking into account the surface gradient. By doing so, one can determine the damage models at a local level. The absence of a surface gradient in the model presents certain advantages. In many physical models, defining an energy density is more straightforward compared to defining the surface gradient. In certain instances, such as bond-based PD, potential function in Molecule Dynamics [42], it is impossible to define the surface gradients, which renders the phase field model infeasible for these applications.

References

[1] Belytschko T, Möes N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering, 2001, 50(4): 993–1013

[2] Belytschko T, Chen H, Xu J, Zi G. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. International Journal for Numerical Methods in Engineering, 2003, 58(12): 1873–1905

[3] Shen Y X, Lew A. An optimally convergent discontinuous Galerkin-based extended finite element method for fracture mechanics. International Journal for Numerical Methods in Engineering, 2010, 82(6): 716–755

[4] Liu W K, Jun S, Zhang Y F. Reproducing kernel particle methods. International Journal for Numerical Methods in Fluids, 1995, 20(8–9): 1081–1106

[5] Belytschko T, Lu Y Y, Gu L. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering, 1994, 37(2): 229–256

[6] Rabczuk T, Areias P M A, Belytschko T. A simplified mesh-free method for shear bands with cohesive surfaces. International Journal for Numerical Methods in Engineering, 2007, 69(5): 993–1021

[7] Rabczuk T, Belytschko T. A three dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[8] Zhang Y M, Zhuang X Y. Cracking elements method for dynamic brittle fracture. Theoretical and Applied Fracture Mechanics, 2019, 102: 1–9

[9] Areias P, Msekh M A, Rabczuk T. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Engineering Fracture Mechanics, 2016, 158: 116–143

[10] Zhang Y M, Lackner R, Zeiml M, Mang H A. Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations. Computer Methods in Applied Mechanics and Engineering, 2015, 287: 335–366

[11] Ruiz G, Pandolfi A, Ortiz M. Three-dimensional cohesive modeling of dynamic mixed-mode fracture. International Journal for Numerical Methods in Engineering, 2001, 52(1–2): 97–120

[12] Radovitzky R, Seagraves A, Tupek M, Noels L. A scalable 3D fracture and fragmentation algorithm based on a hybrid, discontinuous galerkin, cohesive element method. Computer Methods in Applied Mechanics and Engineering, 2011, 200(1–4): 326–344

[13] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering, 2010, 83(10): 1273–1311

[14] Wu J Y. A unified phase-field theory for the mechanics of damage and quasi-brittle failure. Journal of the Mechanics and Physics of Solids, 2017, 103: 72–99

[15] Zhou S W, Zhuang X Y, Zhu H H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theoretical and Applied Fracture Mechanics, 2018, 96: 174–192

[16] Bie Y H, Ren H L, Yan H H, Chen J Y. The unified nonlocal peridynamics-based phase-field damage theory. Theoretical and Applied Fracture Mechanics, 2023, 126: 103980

[17] Silling S A, Epton M, Weckner O, Xu J, Askari E. Peridynamic states and constitutive modeling. Journal of Elasticity, 2007, 88(2): 151–184

[18] Ren H L, Zhuang X Y, Cai Y C, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[19] Ren H L, Zhuang X Y, Rabczuk T. A nonlocal operator method for solving partial differential equations. Computer Methods in Applied Mechanics and Engineering, 2020, 358: 112621

[20] Ren H L, Zhuang X Y, Rabczuk T. A higher order nonlocal operator method for solving partial differential equations. Computer Methods in Applied Mechanics and Engineering, 2020, 367: 113132

[21] Francfort G A, Marigo J J. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids, 1998, 46(8): 1319–1342

[22] Karma A, Kessler D A, Levine H. Phase-field model of mode III dynamic fracture. Physical Review Letters, 2001, 87(4): 045501

[23] Ambrosio L, Tortorelli V M. Approximation of functional depending on jumps by elliptic functional via T-convergence. Communications on Pure and Applied Mathematics, 1990, 43(8): 999–1036

[24] Ulmer H, Hofacker M, Miehe C. Phase field modeling of brittle and ductile fracture. Proceedings in Applied Mathematics and Mechanics, 2013, 13(1): 533–536

[25] Wheeler M F, Wick T, Wollner W. An augmented-Lagrangian method for the phase-field approach for pressurized fractures. Computer Methods in Applied Mechanics and Engineering, 2014, 271: 69–85

[26] Aldakheel F, Hudobivnik B, Hussein A, Wriggers P. Phase-field modeling of brittle fracture using an efficient virtual element scheme. Computer Methods in Applied Mechanics and Engineering, 2018, 341: 443–466

[27] de Borst R, Verhoosel C V. Gradient damage vs phase-field approaches for fracture: Similarities and differences. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 78–94

[28] SimoJ CHughes T J R. Computational Inelasticity. New York: Springer New York, 2006

[29] Hackl K, Fischer F D. On the relation between the principle of maximum dissipation and inelastic evolution given by dissipation potentials. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2008, 464(2089): 117–132

[30] KachanovL. Introduction to Continuum Damage Mechanics. Dordrecht: Springer Dordrecht, 1986

[31] OchsnerA. Continuum Damage Mechanics. Berlin: Springer, 2016

[32] ZhangW HCai Y Q. Continuum Damage Mechanics and Numerical Applications. Berlin: Springer, 2010

[33] Feng X Q, Yu S W. Damage micromechanics for constitutive relations and failure of microcracked quasi-brittle materials. International Journal of Damage Mechanics, 2010, 19(8): 911–948

[34] Ren X D, Chen J S, Li J, Slawson T R, Roth M J. Micro-cracks informed damage models for brittle solids. International Journal of Solids and Structures, 2011, 48(10): 1560–1571

[35] MurakamiS. Continuum Damage Mechanics: A Continuum Mechanics Approach to the Analysis of Damage and Fracture. Dordrecht: Springer Dordrecht, 2012

[36] Jirásek M, Zeman J. Localization study of a regularized variational damage model. International Journal of Solids and Structures, 2015, 69: 131–151

[37] Griffith A A. The phenomena of rupture and flow in solids. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1921, 221: 163–198

[38] Borden M J, Hughes T J R, Landis C M, Anvari A, Lee I J. A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 130–166

[39] Wolfram S. The mathematica book. Assembly Automation, 1999, 19(1): 77

[40] Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209

[41] Ren H L, Zhuang X Y, Fu X L, Li Z Y, Rabczuk T. Bond-based nonlocal models by nonlocal operator method in symmetric support domain. Computer Methods in Applied Mechanics and Engineering, 2024, 418: 116230

[42] Varma S S, Mangipudi K R, Budarapu P R. A coupled quantum-molecular mechanics approach for performance analysis of defective silicon based photovoltaic solar cells. Physica Scripta, 2023, 98(3): 035007

[43] BatheK J. Finite element method. In: Wah B W, ed. Wiley Encyclopedia of Computer Science and Engineering. Hoboken, NJ: John Wiley & Sons, Inc., 2007, 1–12

[44] ReddyJ N. Introduction to the Finite Element Method. New York: McGraw-Hill Education, 2019

[45] Wu Y C, Xiao J Z. Implementation of the multiscale stochastic finite element method on elliptic PDE problems. International Journal of Computational Methods, 2017, 14(1): 1750003

[46] Hughes T J R, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39–41): 4135–4195

[47] Nguyen V P, Anitescu C, Bordas S P A, Rabczuk T. Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation, 2015, 117: 89–116

[48] LiuG RGu Y T. An Introduction to Meshfree Methods and Their Programming. Dordrecht: Springer Dordrecht, 2005

[49] Ren H L, Zhuang X Y, Rabczuk T, Zhu H H. Dual-support smoothed particle hydrodynamics in solid: Variational principle and implicit formulation. Engineering Analysis with Boundary Elements, 2019, 108: 15–29

[50] Bie Y H, Ding K J, Zhao Z F, Wei Y G. The adaptive coupling of dual-horizon peridynamic element and finite element for the progressive failure of materials. International Journal of Fracture, 2024, 245: 89–114

[51] Hudobivnik B, Aldakheel F, Wriggers P. A low order 3D virtual element formulation for finite elasto–plastic deformations. Computational Mechanics, 2019, 63(2): 253–269

[52] Raissi M, Perdikaris P, Karniadakis G E. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 2019, 378: 686–707

[53] Samaniego E, Anitescu C, Goswami S, Nguyen-Thanh V M, Guo H W, Hamdia K, Zhuang X Y, Rabczuk T. An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications. Computer Methods in Applied Mechanics and Engineering, 2020, 362: 112790

[54] Siruvuri S V, Budarapu P R, Paggi M. Influence of cracks on fracture strength and electric power losses in silicon solar cells at high temperatures: Deep machine learning and molecular dynamics approach. Applied Physics A: Materials Science & Processing, 2023, 129(6): 408

[55] Miehe C, Lambrecht M. Algorithms for computation of stresses and elasticity moduli in terms of Seth–Hill’s family of generalized strain tensors. Communications in Numerical Methods in Engineering, 2001, 17(5): 337–353

[56] WinklerB J. Traglastuntersuchungen von unbewehrten und bewehrten Betonstrukturen auf der Grundlage eines objektiven Werkstoffgesetzes für Beton. Dissertation for the Doctoral Degree. Innsbruck: University of Innsbruck, 2001

[57] Gravouil A, Möes N, Belytschko T. Non-planar 3D crack growth by the extended finite element and level sets––Part II: Level set update. International Journal for Numerical Methods in Engineering, 2002, 53(11): 2569–2586

RIGHTS & PERMISSIONS

The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn