Phase-field modeling of fracture for quasi-brittle materials

TL;DR

This paper develops a phase-field model for fracture in quasi-brittle materials that incorporates a damage threshold, providing a more realistic representation of crack initiation and growth compared to traditional models.

Contribution

It introduces a new phase-field fracture model with a damage threshold, systematically formulated within a variational framework and validated through 2D and 3D benchmark simulations.

Findings

The model ensures an elastic stage before damage initiation.

Good agreement with existing literature demonstrates the model's validity.

Inclusion of a damage threshold improves fracture prediction accuracy.

Abstract

This paper addresses the modeling of fracture in quasi-brittle materials using a phase-field approach to the description of crack topology. Within the computational mechanics community, several studies have treated the issue of modeling fracture using phase fields. Most of these studies have used an approach that implies the lack of a damage threshold. We herein explore an alternative model that includes a damage threshold and study how it compares with the most popular approach. The formulation is systematically explained within a rigorous variational framework. Subsequently, we present the corresponding three-dimensional finite element discretization that leads to a straightforward numerical implementation. Benchmark simulations in two dimensions and three dimensions are then presented. The results show that while an elastic stage and a damage threshold are ensured by the present model, good agreement with the results reported in the literature can be obtained, where such features are generally absent.

Tables

Table 1. Table 1 :  Global primary variables

Primary Variable		Field	Type
$𝒖$	displacement field	vector	observable
$\alpha$	damage	scalar	internal

Table 2. Table 2 :  Constitutive state variables

State Variable		Field	Type
$ϵ$	total strain	second-order tensor	observable
$\alpha$	damage	scalar	internal
$\nabla \alpha$	damage gradient	vector	internal

Equations

$$\mathcal{W}(\mathbf{u},\alpha )=\mathcal{P}(\mathbf{u},\alpha )+\mathcal{D}(\alpha ).$$

$$\boldsymbol{\epsilon}=\boldsymbol{\epsilon}(\boldsymbol{u})=\frac{1}{2}\bigg{(}\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{\mathrm{T}}\bigg{)},$$

$$\alpha \in [0,1]\text{with}\dot{\alpha }\ge 0,$$

$$\mathcal{E}(\mathbf{ϵ},\alpha )={\int }_{\Omega }\Psi (\mathbf{ϵ},\mathbf{\alpha })\mathrm{d}\Omega .$$

$$\mathbf{\sigma }(\mathbf{ϵ},\alpha )=\frac{\partial \Psi (\mathbf{ϵ},\mathbf{\alpha })}{\partial \mathbf{ϵ}}.$$

$$\Phi =w^{\prime }(\alpha )\dot{\alpha }+{\partial }_t\left(\frac{1}{2}{\eta }^2\nabla \alpha \cdot \nabla \alpha \right).$$

$$\mathcal{D}(\alpha )={\int }_{\Omega }\left(w(\alpha )+\frac{1}{2}{\eta }^2\nabla \alpha \cdot \nabla \alpha \right)\mathrm{d}\Omega .$$

$$\mathcal{W}(\mathbf{u},\alpha )=\mathcal{E}(\mathbf{u},\alpha )+\mathcal{D}(\alpha )={\int }_{\Omega }\Psi (\mathbf{ϵ},\alpha )\mathrm{d}\Omega +{\int }_{\Omega }\left(w(\alpha )+\frac{1}{2}{\eta }^2\nabla \alpha \cdot \nabla \alpha \right)\mathrm{d}\Omega .$$

$${\alpha }_{t_{n+1}}(\mathbf{x})=\{\begin{array}{ll}{\alpha }_{n+1}(\mathbf{x}),& \text{if}{\alpha }_{n+1}(\mathbf{x})\ge {\alpha }_n(\mathbf{x})\\ {\alpha }_n(\mathbf{x}),& \text{if}{\alpha }_{n+1}(\mathbf{x})<{\alpha }_n(\mathbf{x}),\end{array}$$

$$\left.\frac{\mathrm{d}}{\mathrm{d}h}\mathcal{W}(\boldsymbol{u}+h\tilde{\boldsymbol{u}},\alpha+h\tilde{\alpha})\right|_{h=0}=\int_{\mathit{\Omega}}\bigg{[}\boldsymbol{\sigma}(\boldsymbol{\epsilon},\alpha):{\boldsymbol{\epsilon}}(\tilde{\boldsymbol{u}})+\bigg{(}\frac{\partial\Psi(\boldsymbol{\epsilon,\alpha})}{\partial{\alpha}}+w^{\prime}(\alpha)\bigg{)}\tilde{\alpha}+\eta^{2}\nabla\alpha\cdot\nabla\tilde{\alpha}\bigg{]}\mathrm{d}{\mathit{\Omega}}\geq 0.$$

$${\int }_{\Omega }\mathbf{\sigma }(\mathbf{ϵ},\alpha ):\mathbf{ϵ}(\tilde{\mathbf{u}})\mathrm{d}\Omega =0,$$

$$\int_{\mathit{\Omega}}\bigg{[}\left(\frac{\partial\Psi(\boldsymbol{\epsilon,\alpha})}{\partial{\alpha}}+w^{\prime}(\alpha)\right)\tilde{\alpha}+\eta^{2}\nabla\alpha\cdot\nabla\tilde{\alpha}\bigg{]}\mathrm{d}{\mathit{\Omega}}\geq 0,$$

$$f_d(\mathbf{u},\alpha )=$$

$$\mathbf{n}\cdot \nabla \alpha \ge 0\text{on}\partial \Omega$$

$$\displaystyle\int_{\mathit{\Omega}}\bigg{[}\boldsymbol{\sigma}(\boldsymbol{\epsilon},\alpha):{\boldsymbol{\epsilon}}(\dot{\boldsymbol{u}})+\bigg{(}\frac{\partial\Psi(\boldsymbol{\epsilon,\alpha})}{\partial{\alpha}}+w^{\prime}(\alpha)-\eta^{2}\nabla\cdot\nabla\alpha\bigg{)}\dot{\alpha}\bigg{]}\mathrm{d}{\mathit{\Omega}}=0.$$

$$f_d(\mathbf{u},\alpha )\dot{\alpha }=0\text{in}\Omega$$

$$\dot{\alpha }\mathbf{n}\cdot \nabla \alpha =0\text{on}\partial \Omega .$$

$$\mathbf{n}\cdot \nabla \alpha =0\text{on}\partial \Omega ,$$

$${\Psi }_0(\mathbf{ϵ})$$

$${\mathbf{C}}_0=\lambda 1\otimes 1+2\mu \mathbf{I},$$

$${\Psi }_0^{+}(\mathbf{ϵ})=\frac{1}{2}K⟨\mathrm{tr}(\mathbf{ϵ}){⟩}_{+}^2+\mu ({\mathbf{ϵ}}_{\mathrm{dev}}:{\mathbf{ϵ}}_{\mathrm{dev}})\text{and}{\Psi }_0^{-}(\mathbf{ϵ})=\frac{1}{2}K⟨\mathrm{tr}(\mathbf{ϵ}){⟩}_{-}^2,$$

$$\Psi (\mathbf{ϵ},\alpha )=f(\alpha ){\Psi }_0^{+}(\mathbf{ϵ})+{\Psi }_0^{-}(\mathbf{ϵ}).$$

$$\mathbf{\sigma }(\mathbf{ϵ},\alpha )=\frac{\partial \Psi }{\partial \mathbf{ϵ}}=$$

$$f^{\prime }(\alpha )<0\text{and}f(1)=0.$$

$$f(\alpha )=(1-\alpha {)}^2.$$

$$w(1)=w_0<\infty .$$

$$w(\alpha )=\{\begin{array}{ll}{w}_0\alpha & \text{model with an elastic stage},\\ w_0{\alpha }^2& \text{model without an elastic stage.}\end{array}$$

$$G_{\mathrm{c}}=2l{\int }_0^1\sqrt{2w_{0}w(\beta )}\mathrm{d}\beta =c_{\mathrm{w}}\frac{l{w}_0}{\sqrt{2}}\text{with}{c}_{\mathrm{w}}=4{\int }_0^1\sqrt{\frac{w_0}{2w(\beta )}}\mathrm{d}\beta ,$$

$$l=\frac{\eta }{\sqrt{w_0}}.$$

$$G_{\mathrm{c}}=\frac{4\sqrt{2}}{3}{w}_{0}l.$$

Full text

††journal: Underground Space

1 Introduction

Crack initiation and propagation in quasi-brittle materials is both of significant scientific interest and paramount importance in engineering applications. Examples of these materials include concrete and rocks that exhibit microcracks that localize in narrow bands [1]. Methods developed in the framework of classical fracture mechanics, typically rooted in the pioneering works of [2], [3] and [4], have been relatively successfull in several engineering applications. Nevertheless, these theories present major limitations, including the inability to naturally describe crack initiation and propagation without initial defects and a prescribed crack path [5]. Moreover, their application requires nontrivial finite element techniques such as local refinement near the crack tip [6]. In view of these drawbacks, a variety of methods have emerged in the past few decades with the objective of providing a more convenient description of fracture.

From a continuum mechanics perspective, the study of post-critical behavior in solids can be approached using softening constitutive models. This leads to strain localization, which signals the appearance of the so-called fracture process zone. Thus, some form of representation of the high gradients in the displacement field is required to describe the material behavior after the onset of strain localization. Classical continuum mechanics cannot handle the ill-posedness of the evolution problem, thus resulting in numerical simulations with a localization zone that tends to measure zero and vanishing dissipated energy as the mesh size decreases [7]. Consequently, modeling strain localization has been a constant challenge in computational mechanics [8]. To overcome these problems, the correct description of crack kinematics, topology, and energetics is required. A possible classification of these descriptions considers two subsets: i) kinematics with sharp discontinuities and ii) regularized kinematics.

Regarding the first subset, discontinuities in the displacement field are explicitly introduced to describe sharp crack kinematics. Well-known examples of this strategy are the strong discontinuity approach (SDA) and the extended finite element method (XFEM). Building upon [9], [8] proposed a continuum SDA to describe failure in geomaterials. In turn, starting with the work of [10], XFEM has been an attractive alternative for the modeling of arbitrary crack paths without remeshing. It has, for instance, been applied to describe crack growth with a cohesive law [11]. Moreover, a method based on local maximum-entropy (LME) interpolation was studied in conjunction with XFEM [12]. The results suggest that the method proposed therein is a competitive alternative in terms of computational cost when compared to the standard XFEM.

Although sharp crack modeling has yielded noteworthy results, major limitations can be highlighted, particularly, the inability to naturally describe complex crack topologies. In addition, the numerical treatment of crack propagation requires some form of crack-path tracking, which can be a very cumbersome task. To mitigate these difficulties, phase-field regularized models have been shown to be a competitive alternative [13], which can be implemented straightforwardly. In these formulations, a continuous variable, namely, the phase-field variable, is used to describe a smooth transition between the damaged/undamaged phases. Moreover, regularization based on phase fields can be viewed as a bridge between damage mechanics and a diffuse approximation of the sharp crack topology, which avoids the need to explicitly model discontinuities in the displacement field. In the context of fracture, the phase-field variable is represented by the scalar-valued damage quantity, whose gradient is introduced in the formulation. Thus, a clear link can be established between phase-field fracture and gradient damage. We refer to [14] for a comparison between gradient-enhanced damage models and the phase-field approach to fracture.

Several studies that apply gradient-based phase-field regularizations for rate-independent systems can be related to the work of [5], where a variational formulation was introduced to overcome the limitations of the Griffith model for brittle fracture, particularly, the need for a priori constraints on the crack topology. This formulation resulted in an energy functional reminiscent of the potential of [15] for image segmentation. Subsequently, to avoid the numerical difficulties imposed by the free discontinuity problem of [5], [16] proposed an energy functional based on phase-field regularization, where a damage gradient term was introduced (although it was not originally viewed as such). This regularization has been shown to converge to the Griffith fracture model through $\Gamma$ convergence [17], and was inspired by the work of [18] for the regularization of the Mumford and Shah potential. The reader is referred to [19] for an overview of the regularized formulation of brittle fracture.

In the computational mechanics community, several contributions to the phase-field modeling of fracture have emerged. In [13], and further developed in [20], an alternative phase-field formulation was proposed, also based on thermodynamic principles. An attractive feature was incorporated in this formulation: the definition of a realistic anisotropic stored energy, obtained by defining the bulk energy density as an additive decomposition of positive (due to tension) and negative (due to compression) contributions. In this setting, damage is allowed to act on the positive component only, disallowing fracture due to compression. This formulation was extended to the dynamic case by [21]. Moreover, ductile behavior has been considered in several works. For instance, [22] combines local $J_{2}$ hardening plasticity with gradient damage and also considers anisotropic damage behavior. Furthermore, [23, 24] extended these formulations to gradient plasticity combined with gradient damage. In [25], fracture in shells was approached using a phase-field model with isogeometric analysis [26]. Phase-field modeling combined with LME interpolation was proposed by [27] for the study of thin shells, and by [28] using a fourth-order phase-field model. The above-mentioned contributions were generally built upon [16], using a particular method of approximation of the Griffith model using elliptic functionals that implies, from a mechanical perspective, the lack of a damage threshold. Therefore, an elastic stage was not included in the resulting evolution.

An interesting alternative conceived from the standpoint of continuum gradient-enhanced damage mechanics was proposed by [29] and [30]. The general and rigorous framework of this approach is detailed in the survey of [31]. An important feature of this framework is that the models considered were related to the work of [16]. Nonetheless, the possibility of having a damage threshold was considered, which is generally not present in most studies using phase fields to model fracture. This option was adopted by [32] and [33, 34], and incorporated in a ductile fracture model. The formulation proposed therein can capture a variety of macroscopic fracture behaviors. In [35] and further in [36], within the same framework, hardening effects and gradient plasticity with variable internal length were considered. The variational approach used in these studies follows a rigorous energetic formulation, formalized by [37] and [38]. As discussed in [39], the energetic formulation presents several advantages with respect to classical formulations. For instance, the governing equations are naturally derived using the calculus of variations and three physical principles: the stability condition, energy balance and the fulfilment of the second law of thermodynamics. Moreover, the definition of a global energy functional leads to a robust numerical implementation that can be solved using a simple staggered scheme.

In this study, we perform an analysis of this alternative method to describe fracture using phase fields for quasi-brittle materials. In Section 2, the governing equations of the problem are systematically derived following the variational framework. Then, Section 3 presents the details of the corresponding finite element implementation. The numerical simulations are presented in Section 4 and are compared with the results reported in the literature, where the most popular model for quasi-brittle materials was used.

2 Formulation

We adopt an energetic framework for the description of the behavior of deformable solids in the rate-independent case [37]. In this work, we assume evolutions under small strains, with the exception of certain localized regions. Our goal is to describe quasi-brittle fracture with an elastic stage, which results in a two-field formulation, considering displacements $\boldsymbol{u}$ and the damage variable $\alpha$ as primary fields.

Following the theory of generalized standard materials [40], an energy functional is defined as the sum of potential and dissipative energy terms:

[TABLE]

The minimization of this functional with respect to $\boldsymbol{u}$ and $\alpha$ separately entails the fulfillment of the momentum balance and damage criterion, respectively. The weak form of each of these equations is naturally obtained and can be discretized using the finite element method. It is noteworthy that the definition of a total energy quantity is not always straightforward: it strongly depends on the form of the dissipated work $\mathcal{D}(\alpha)$ [34].

The fact that gradient-enhanced damage models entail a regularization of the softening problem that leads to mesh-independent solutions with nonvanishing dissipation is well known. Our approach follows the process of phase-field models, which are strongly linked to gradient-enhanced damage. In fact, the primary difference between both formulations is how they are conceived: the primary idea of phase-field models is the description of the discontinuity of the crack using a continuous field, while gradient damage models are approached from a mechanical perspective, where gradients are included to regularize the ill-posed boundary value problem. Nevertheless, an important difference is that models developed in a purely gradient-damage framework typically result in broadening of the damage zone, whereas phase-field models capture a sharp transition zone in more naturally owing to the use of a degradation function [14]. Thus, in this study, the phase-field approach is employed, but the constitutive functions are defined from the perspective of damage mechanics to preserve the mechanical interpretation.

The ingredients to establish the energetic formulation are described in the following sections.

2.1 Primary and state variables

Consider a solid body ${\mathit{\Omega}}$ with Neumann boundary $\partial{\mathit{\Omega}}_{\sigma}$ and Dirichlet boundary $\partial{\mathit{\Omega}}_{u}$. Let $\boldsymbol{u}(\boldsymbol{x},t)$ be the displacement of a point $\boldsymbol{x}\in{\mathit{\Omega}}$ at time $t$. Because the small strain hypothesis is adopted in this work, no distinction is made between the original and the deformed configuration of the body. Consequently, the second-order total strain tensor is defined as

[TABLE]

where we have dropped the explicit dependence of the involved fields on $\boldsymbol{x}$ and $t$ for the sake of notational simplicity.

For the phase-field description of the crack topology, the internal scalar-valued damage variable $\alpha$ is characterized by

[TABLE]

which indicates the damaged/undamaged points in the solid. A value of $\alpha=0$ corresponds to an undamaged material state, while $\alpha=1$ defines a completely broken material state. In the formulation presented herein, regularization is attained using the gradient of the damage variable $\nabla\alpha$, which allows for the formation of a diffuse crack whose width is finite and depends on the damage internal length scale $\eta$. Figure 1 shows this type of crack description within the general problem setting.

The above-mentioned variables are fundamental in the multifield model. As summarized in Table 1, the primary global variables for the model are $\boldsymbol{u}$ and $\alpha$. Meanwhile, the constitutive state variables, which define the material behavior of each point within the solid, are $\boldsymbol{\epsilon}$, $\alpha$ and $\nabla\alpha$ (see Table 2).

2.2 Variational approach

2.2.1 Energy functional

This section describes the energetic formulation in the context of gradient damage. Adopting the theory of generalized standard materials [40, 41], the variational approach to fracture begins with the description of the basic energetic quantities. The global stored energy is defined using the elastic energy density $\Psi$ as

[TABLE]

The Cauchy stress tensor $\boldsymbol{\sigma}$ is given by the constitutive equation

[TABLE]

The dissipative nature of the internal variables is characterized by the definition of the dissipation potential [34]. Considering rate-independence, the dissipation potential is a first-order homogeneous convex function of the rate of the damage variable, and can be expressed as

[TABLE]

The function $w(\alpha)$ represents the dissipated energy in a homogeneous damage process. In addition, the term corresponding to the damage gradient with its corresponding internal length scale is included. This parameter is directly related to the size of the localization zone, this implying several advantages, including the ability to describe structural stability and size effects [31]. With Eq. (6), assuming a smooth evolution, the dissipated work $\mathcal{D}$ follows:

[TABLE]

In this work, no external forces are considered (although their inclusion would be straightforward), and displacements are imposed on $\partial{\mathit{\Omega}}_{u}$. Thus, the total energy functional is defined using Eq. (7) and Eq. (4) as

[TABLE]

2.2.2 Energetic formulation

The building blocks of the variational framework adopted in this work are the following principles [30, 32]:

Stability condition, 2. 2.

Energy balance, 3. 3.

Irreversibility condition.

The irreversibility condition is imposed on the damage variable to disallow material healing. It is applied numerically by considering the damage value corresponding to the previous load step as the minimum admissible level of damage for a given position in the body as follows:

[TABLE]

where $\alpha_{{n+1}}(\boldsymbol{x})$ is the damage value for the current time step ${n+1}$ in $\boldsymbol{x}\in{\mathit{\Omega}}$.

Stability condition

A directional stability condition may be defined by the Taylor expansion of the (sufficiently regular) energy functional. The first-order variations of Eq. (8) yield the first-order stability condition:

[TABLE]

We refer the reader to [42], [30], and [31] for more details on the stability condition in the context of gradient-enhanced damage and to [37] and [38] for more general definitions.

From Eq. (10), the following results are obtained.

For $\tilde{\alpha}=0$:

[TABLE]

representing the weak form of the equilibrium equation in the absence of external loads.

• 2.

For $\tilde{\boldsymbol{u}}=\boldsymbol{0}$:

[TABLE]

in which the weak form of the damage criterion is obtained. The last term is integrated by parts and the gradient-dependent damage yield criterion is recovered in the local form

[TABLE]

Moreover, it can be easily shown that the boundary conditions

[TABLE]

also follow, where $\boldsymbol{n}$ is the unit outward normal vector to $\partial{\mathit{\Omega}}$.

Energy balance

The energy balance states the need for the total energy to remain constant as the state variables evolve. Thus, it is essentially a restatement of the first law of thermodynamics. Following a procedure analogous to the treatment of the stability condition, the energy balance leads to

[TABLE]

For $\dot{\boldsymbol{u}}=\boldsymbol{0}$, and using Eq. (13), the damage consistency conditions

[TABLE]

are obtained, along with the boundary conditions

[TABLE]

Thus, for any $\dot{\alpha}>0$,

[TABLE]

representing the additional boundary conditions required for the gradient-enhanced evolution equations.

Eqs. (13), (16) and the irreversibility condition represent gradient-enhanced versions of the evolution equations of classical local models in the form of Karush-Kuhn-Tucker conditions.

2.3 Constitutive assumptions

This section defines the constitutive assumptions for the damaging process. The following definitions represent the functions and parameters that are intended to capture the behavior of quasi-brittle materials. The resulting model should describe the softening behavior and strain localization as the phase-field variable evolves.

The initial elastic energy density, corresponding to a sound material state, is expressed as

[TABLE]

where $\boldsymbol{\sigma}_{0}$ is the undamaged Cauchy stress tensor. The symbols $\lambda$ and $\mu$ are the Lamé parameters, and $\mathbf{C}_{0}$ is the undamaged elastic fourth-order tensor, given for an isotropic elastic material by

[TABLE]

where $\boldsymbol{1}$ is the second-order identity tensor and $\boldsymbol{I}$ is the fourth-order symmetric identity tensor.

A differentiated damaging behavior is considered by decomposing the elastic energy density into positive (due to tension) and negative (due to compression) energies. Damage is allowed to act on the positive part of the elastic energy only, thus disallowing failure due to compression. Consequently, fracture only takes place in regions under tension, while crack interpenetration is disallowed in the compressed regions [29]. This energy split allows to describe quasi-brittle fracture in materials that exhibit different strengths in tension and compression. Examples of these are concrete and other geomaterials, where cracking is associated with tensile stresses. A decomposition of stress into tension and compression was carried out by [43] to develop a local damage model for concrete structures under cyclic loading conditions. In the context of phase-field fracture, [29] proposed an energy decompositon based on a volumetric-deviatoric split. Alternative approaches have been proposed in the literature, such as a spectral decomposition by [13], a no-tension split by [44], and recently, a stress-based split by [45]. Although each approach owns advantages in certain cases, the definition of an optimal split remains an open issue. In this work, we adopt the volumetric-deviatoric model of [29]. Therein, the elastic energy density is expressed in terms of the volumetric and deviatoric components, and the positive and negative contributions are expressed as

[TABLE]

where the ramp function \langle\operatorname{tr}(\boldsymbol{\epsilon})\rangle_{\pm}=\frac{1}{2}\big{(}\operatorname{tr}(\boldsymbol{\epsilon})\pm\lvert\operatorname{tr}(\boldsymbol{\epsilon})\rvert\big{)} is used, $K=\lambda+\frac{2}{3}\mu$ is the bulk modulus and $\boldsymbol{\epsilon}_{\mathrm{dev}}$ denotes the deviatoric part of the strain tensor. Anisotropic material degradation can now be described by a stored energy density of the form

[TABLE]

Using the Heaviside step function $\operatorname{\mathcal{H}}$, the Cauchy stress tensor can be decomposed as follows:

[TABLE]

Regarding the elastic energy degradation, we require that

[TABLE]

The first condition ensures material degradation, while the second defines a fractured material state. For this purpose, we adopt the quadratic function of damage that has been widely used in phase-field models in the literature:

[TABLE]

The main advantage of using a quadratic degradation function is that the damage criterion is a (constrained) linear partial differential equation, as will be clear in the following section. Advantages of using a cubic degradation function were reported by [46]; however, we adopt the quadratic function to maintain our claim of a simple numerical implementation.

The dissipation due to local damage evolution is defined by the positive-valued function $w(\alpha)$ that represents the dissipated energy of a volume element throughout the damage process. As described by [31], this dissipation is set to increase to a critical value

[TABLE]

Among the existing phase-field models, two alternatives can be found for the local damage dissipation:

[TABLE]

The critical damage dissipation $w_{0}>0$ represents the energy dissipated during a complete damage process for a volume element, and is related to the fracture toughness $G_{\mathrm{c}}$ used in other formulations by [31]:

[TABLE]

where the internal length is denoted by $l$ and can be expressed in terms of $\eta$ as

[TABLE]

From the previous expressions and the energetic definitions (7) and (8), note that $\eta$ has the units of $\text{energy}^{1/2}\text{length}$.

The model without an elastic stage has been consistently used in phase-field models, as shown in [22, 47] and [21, 46]. For this choice, the absence of an elastic stage is due to $w^{\prime}(0)=0$. In this work, we adopt $w(\alpha)=w_{0}\alpha$ as the constitutive choice, for which $w^{\prime}(0)=w_{0}>0$. Therefore, this function allows for the existence of an initial elastic stage before the damage evolution. Moreover, the parameter $w_{0}$ can be considered as a damage threshold from the viewpoint of classical continuum damage mechanics. We consider this an attractive feature because it forbids the onset of damage before a critical stress state is reached. From this choice and using Eq. (28), $w_{0}$ can be related to the fracture toughness by

[TABLE]

With these definitions, the damage criterion from Eq. (13) becomes

[TABLE]

2.4 Alternate minimization

With the definition of the main ingredients of the variational approach, a numerical solution can be readily obtained. For this purpose, an alternate minimization algorithm is applied, which naturally emerges from the energetic principles. This procedure takes advantage of the fact that although the global energy is nonconvex, it is convex with respect to $\boldsymbol{u}$ and $\alpha$ individually. Introducing the constitutive assumptions of the previous section into Eq. (8), the global energy functional reads

[TABLE]

The alternate minimization follows.

Minimization with respect to the displacement field:

[TABLE]

• 2.

Minimization with respect to the damage field:

[TABLE]

The last two expressions can be easily viewed as weak forms of the Euler equations of the underlying global minimization problem. In this work, the staggered solution is applied incrementally, which results in a straightforward implementation, given the choice of constitutive functions. The solution of the problem is subsequently completed with the imposition of the irreversibility condition in the algorithmic form (9). The numerical solution is described in the following section.

3 Numerical solution and implementation

This section presents a straightforward implementation of the numerical solution of Eqs. (33) and (34). Linear tetrahedral finite elements are used throughout this section, although the use of higher-order elements would be equally straightforward.

The overall procedure can be summarized as follows: Eqs. (33) and (34) are first written in discrete form using the Voigt notation by projection over finite elements. The elastic equilibrium equation is subsequently used to obtain the displacement field, and the damage field is obtained using the updated displacements. This process is repeated iteratively to approximate the two primary fields at a specific pseudo-time. The discrete versions are then solved at the following pseudo-time, entering a temporal incremental scheme.

3.1 Numerical solution

3.1.1 Finite element approximation

The displacement field approximation is defined as

[TABLE]

where the global vector of nodal displacements $\mathbf{u}$ and the shape function matrix $\mathbf{N}_{\mathrm{v}}$ are obtained using an assembly operator $\operatorname*{\scalerel*{\mathsf{A}}{\sum}}$:

[TABLE]

and

[TABLE]

The superscript ${el}$ denotes an element vector or matrix, and the index $i\in[1,2,3,4]$ indicates the node number with $N_{i}$ representing the corresponding shape functions. The vector $\mathbf{u}_{i}^{el}$ contains the nodal displacements in each spatial direction.

The strain vector inside an element is defined as

[TABLE]

where $\gamma$ represents the shear strain components. The global strain vector is obtained from the nodal displacements as

[TABLE]

where

[TABLE]

The damage and gradient damage fields are approximated by

[TABLE]

and

[TABLE]

where the assembly operator is used again:

[TABLE]

and

[TABLE]

Hereinafter, all element vectors and matrices are assumed to be assembled into their corresponding global forms through $\operatorname*{\scalerel*{\mathsf{A}}{\sum}}$.

3.1.2 Discrete forms

The discrete versions of the evolution equations consist of systems of linear equations, from where the solution to each primary variable can be easily obtained. Before presenting the discrete forms, we introduce the definitions that follow. The volumetric strain is defined as

[TABLE]

and the deviatoric part of the strain vector inside an element reads

[TABLE]

For the anisotropic model, the constitutive matrix is expressed in terms of the bulk and shear moduli. Introducing the volumetric and deviatoric operators

[TABLE]

the local elastic matrix is defined in the decomposed form

[TABLE]

where $\operatorname{\mathcal{H}}$ is the Heaviside step function. The global damaged stress vector can than be readily obtained from

[TABLE]

Using the previous definitions, the discrete form of Eq. (33) can be expressed globally as

[TABLE]

from where the global displacement vector $\mathbf{u}$ can be readily obtained after imposing boundary conditions.

Finally, the global damage field vector $\boldsymbol{\alpha}$ is obtained from the discrete version of Eq. (34), which is given by

[TABLE]

where the element form of $\boldsymbol{\upsigma}_{0}^{+}$ reads

[TABLE]

3.2 Implementation details

The implementation of the proposed model is summarized in this section. We emphasize that our goal in this study is not improving computational efficiency, but to present a clear and straightforward implementation (in the spirit of [48]) of a rather complex problem.

The overall architecture of the proposed implementation consists of the following modules:

Preprocessing: generating the mesh. Details on mesh generation are outside the scope of this paper.

• 2.

Data input: reading data files and defining material constants, convergence factors and imposed displacements.

• 3.

Initialization: initializing convergence vectors and storage matrices.

• 4.

One module for each primary variable: obtaining the global vectors $[\mathbf{u},\boldsymbol{\alpha}]$ using two separate finite-element problems.

• 5.

Main module: it sequentially calls the modules for each primary variable into an iterative procedure for each load step.

• 6.

Postprocessing: displaying the primary fields illustratively.

3.2.1 Modules for primary variables

The nodal vectors corresponding to the global primary fields, $\mathbf{u}$ and $\boldsymbol{\alpha}$, are obtained in separate modules: elast3D and dam3D, respectively. All of them share the same general structure. Given that Eqs. (50) and (51) are linear, the subroutines consist of standard finite element procedures. Equations (50) and (51) exhibit the general form:

[TABLE]

where $\mathbf{x}$ is a vector containing the nodal values of a primary variable. Therefore, the solution consists of determining the local coefficient matrix $\mathbf{K}^{el}$ and the local right hand side $\mathbf{b}^{el}$. Subsequently, the global coefficient matrix $\mathbf{K}$ and the global right hand side $\mathbf{b}$ are assembled, from where the solution is obtained after imposing boundary conditions.

3.2.2 Main module

The numerical setting results in an incremental staggered algorithm, as described in Algorithm 1. For $n_{\mathrm{desp}}$ incremental displacements imposed on $\partial{\mathit{\Omega}}_{u}$, $[\mathbf{u}_{{n+1}},\boldsymbol{\alpha}_{{n+1}}]$ are found iteratively using two independent tolerances $t_{\mathrm{u}}$ and $t_{\mathrm{d}}$, for a maximum number of iterations $k_{\mathrm{max}}$.

4 Numerical simulations

A square specimen with a horizontal notch in the middle is adopted for experiments involving both tension and shear loading. Monotonic displacements are imposed in the numerical experiments, as schematized in Fig. 2. For both tension and shear loading, the following parameters were adopted:

$K=121030$ MPa;

• 2.

$\nu=0.227$;

• 3.

$w_{0}=75.94$ MPa mm and

• 4.

$\eta=0.052$ (MPa mm)1/2 mm.

Experiment I

The tension case is considered first, with $\bar{u}_{x}=0$ and $\bar{u}_{y}\geq 0$ (Fig. 2). Vertical displacements are imposed on the top boundary from 0 to $6\text{\times}{10}^{-3}$ mm, with increments of $1\text{\times}{10}^{-4}$ mm, while the bottom boundary is fixed in both directions.

The crack propagation and the deformed specimen are shown in Figure 3. A single crack branch is initiated and propagates horizontally. As expected in brittle fracture, the specimen experiences a catastrophic-like failure mode, where the crack is initiated at the tip of the notch and propagates horizontally after a few load steps. This is reflected in the force-displacement curve shown in Figure 4, where an abrupt drop in the load-carrying capacity is observed. The simulation was performed using 3401 bilinear quadrilateral elements. As shown in Figure 4, the results closely resemble the experiment of [20], where a finer mesh of 20000 linear triangles was used.

For the same experiment, a three-dimensional (3D) simulation is shown in Figure 5. A uniform thickness of $0.1$ mm was used for the model, which was discretized with a mesh of 16500 hexahedral elements. The same displacements as in the two-dimensional simulation were imposed on the top boundary, while the bottom and lateral bottom boundary was fixed in all three directions and the lateral boundaries were fixed in the $z$ direction. In Fig. 5, the abrupt crack evolution can be observed better. Figure 6 shows the resulting force-displacement curve, which is compared with the results obtained in [49], where the implementation of the phase-field model without a damage threshold was performed in the commercial software Abaqus using 96896 hexahedral elements.

Experiment II

We subject the same square specimen with a notch from Experiment I to shear loading, with $\bar{u}_{x}\geq 0$ and $\bar{u}_{y}=0$ (Fig. 2). Incremental horizontal displacements are imposed on the top boundary, while the bottom boundary is fixed. Displacements are imposed from [math] to $12.7\text{\times}{10}^{-3}$ mm, with increments of $1\text{\times}{10}^{-5}$ mm. The lateral boundaries are fixed in the vertical direction. For the simulation, a mesh of 3532 bilinear quadrilateral elements was used.

The crack propagation and the deformed specimen shown in Fig. 7 are a direct result of the decomposition into positive (due to tension) and negative (due to compression) energies in Eq. (21). A single crack branch is initiated and propagates through regions of intense positive stress. Along with the force-displacement curve shown in Fig. 8, these results resemble the results of [20] and [21]. The former used 30000 linear triangles, while the latter applied cubic T-splines with 5587 cubic basis functions.

5 Conclusions

We herein presented a multidimensional model for the description of brittle fracture using the phase-field approach, which was written based on a consistent energetic formulation. As opposed to the most popular approaches of phase-field fracture, we adopted constitutive assumptions in the spirit of gradient-enhanced damage models, which allow for a damage threshold. We consider this to be an attractive feature because an elastic stage is included in the deformation process. The simulations presented indicate good agreement with the results obtained in the literature using the model without an elastic stage.

One of the main features of the formulation explored in this work is its variational character. This allows for an elegant approach, resulting in a powerful tool to describe the behavior of solids undergoing brittle fracture. Moreover, the resulting numerical problem can be solved in a relatively straightforward manner. We have explicitly presented a simple 3D finite element formulation that could be easily implemented in standard finite element routines. Although the staggered approach is known to result in slow convergence when compared to monolithic solution schemes, no major convergence issues were observed in the performed simulations other than the expected increase in iterations at the onset of the softening regime.

The obvious extension of this study is the introduction of ductile behavior, which is the subject of another article. In addition, since high run times were experienced, particularly for 3D simulations, efficient numerical implementations are crucial for further developments, which is a topic of further study.

Acknowledgements

We gratefully acknowledge the financial support of CYTED (Ibero- American Program to Promote Science and Technology) through the CADING network (Ibero-American Network for High Performance Computing in Engineering). The support of the EU H2020-MSCA-RISE-2016 project “BESTOFRAC: Environmentally best practices and optimisation in hydraulic fracturing for shale gas/oil development” is also acknowledged.