Phase-field modeling of fracture for quasi-brittle materials

02/27/2019 ∙ by Jacinto Ulloa, et al. ∙ Inria 0

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.



There are no comments yet.


page 18

page 19

page 20

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 (comi1999, ). Methods developed in the framework of classical fracture mechanics, typically rooted in the pioneering works of griff , irwin1957 and barenblatt1962 , 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 (FrancMar1998, ). Moreover, their application requires nontrivial finite element techniques such as local refinement near the crack tip (nguyen2013, ). 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 (DeBorstMull1992, ). Consequently, modeling strain localization has been a constant challenge in computational mechanics (OliHueSamCha2004, ). 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 SimOliArm1993 , OliHueSamCha2004 proposed a continuum SDA to describe failure in geomaterials. In turn, starting with the work of MoeDolBely1999 , 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 (MoeBely2002, )

. Moreover, a method based on local maximum-entropy (LME) interpolation was studied in conjunction with XFEM

(amiri2014, ). 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 (MieWelHof2010, ), 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 deborst2016 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 FrancMar1998 , 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 MumfShah1989 for image segmentation. Subsequently, to avoid the numerical difficulties imposed by the free discontinuity problem of FrancMar1998 , BourFrancMar2000 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 convergence (DalToa2002, ), and was inspired by the work of AmbTor1990 for the regularization of the Mumford and Shah potential. The reader is referred to BourFrancMar2008 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 MieWelHof2010 , and further developed in MieHofWel2010 , 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 Bord-Hugh2012 . Moreover, ductile behavior has been considered in several works. For instance, AmbGerDeL2015 combines local hardening plasticity with gradient damage and also considers anisotropic damage behavior. Furthermore, MieTeiAld2016 ; MieAldTei2016 extended these formulations to gradient plasticity combined with gradient damage. In ambati2016shells , fracture in shells was approached using a phase-field model with isogeometric analysis (hughesiso, ). Phase-field modeling combined with LME interpolation was proposed by amiri2014phase for the study of thin shells, and by amiri2016 using a fourth-order phase-field model. The above-mentioned contributions were generally built upon BourFrancMar2000 , 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 AmorMarMau2009 and pham2011 . The general and rigorous framework of this approach is detailed in the survey of MarMauPham2016 . An important feature of this framework is that the models considered were related to the work of (BourFrancMar2000, ). 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 Alessi2013 and AleMarVid2014 ; AleMarMauVid2017 , and incorporated in a ductile fracture model. The formulation proposed therein can capture a variety of macroscopic fracture behaviors. In UllRodSam2016 and further in rodriguez2018 , 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 Mielke2006 and Mielke2015 . As discussed in Alessi2016 , 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 (Mielke2006, ). 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 and the damage variable as primary fields.

Following the theory of generalized standard materials (halphen1975, ), an energy functional is defined as the sum of potential and dissipative energy terms:


The minimization of this functional with respect to and 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 (AleMarMauVid2017, ).

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 (deborst2016, ). 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

Figure 1:  Diffuse crack description.

Consider a solid body with Neumann boundary and Dirichlet boundary . Let be the displacement of a point at time

. 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


where we have dropped the explicit dependence of the involved fields on and for the sake of notational simplicity.

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


which indicates the damaged/undamaged points in the solid. A value of corresponds to an undamaged material state, while defines a completely broken material state. In the formulation presented herein, regularization is attained using the gradient of the damage variable , which allows for the formation of a diffuse crack whose width is finite and depends on the damage internal length scale . 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 and . Meanwhile, the constitutive state variables, which define the material behavior of each point within the solid, are , and (see Table 2).

Primary Variable Field Type
displacement field vector observable
damage scalar internal
Table 1:  Global primary variables
State Variable Field Type
total strain second-order tensor observable
damage scalar internal
damage gradient vector internal
Table 2:  Constitutive state variables

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 (halphen1975, ; germ, ), 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 as


The Cauchy stress tensor is given by the constitutive equation


The dissipative nature of the internal variables is characterized by the definition of the dissipation potential (AleMarMauVid2017, ). 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


The function 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 (MarMauPham2016, ). With Eq. (6), assuming a smooth evolution, the dissipated work follows:


In this work, no external forces are considered (although their inclusion would be straightforward), and displacements are imposed on . Thus, the total energy functional is defined using Eq. (7) and Eq. (4) as


2.2.2 Energetic formulation

The building blocks of the variational framework adopted in this work are the following principles (pham2011, ; Alessi2013, ):

  1. Stability condition,

  2. Energy balance,

  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:


where is the damage value for the current time step in .

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:


We refer the reader to pham_mara , pham2011 , and MarMauPham2016 for more details on the stability condition in the context of gradient-enhanced damage and to Mielke2006 and Mielke2015 for more general definitions.

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

  • For :


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

  • For :


    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


    Moreover, it can be easily shown that the boundary conditions


    also follow, where is the unit outward normal vector to .

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


For , and using Eq. (13), the damage consistency conditions


are obtained, along with the boundary conditions


Thus, for any ,


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


where is the undamaged Cauchy stress tensor. The symbols and are the Lamé parameters, and is the undamaged elastic fourth-order tensor, given for an isotropic elastic material by


where is the second-order identity tensor and 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 (AmorMarMau2009, ). 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 faria2004 to develop a local damage model for concrete structures under cyclic loading conditions. In the context of phase-field fracture, AmorMarMau2009 proposed an energy decompositon based on a volumetric-deviatoric split. Alternative approaches have been proposed in the literature, such as a spectral decomposition by MieWelHof2010 , a no-tension split by freddi2010 , and recently, a stress-based split by steinke2018 . 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 AmorMarMau2009 . Therein, the elastic energy density is expressed in terms of the volumetric and deviatoric components, and the positive and negative contributions are expressed as


where the ramp function is used, is the bulk modulus and denotes the deviatoric part of the strain tensor. Anisotropic material degradation can now be described by a stored energy density of the form


Using the Heaviside step function , the Cauchy stress tensor can be decomposed as follows:


Regarding the elastic energy degradation, we require that


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:


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

Bord-Hugh2016 ; 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 that represents the dissipated energy of a volume element throughout the damage process. As described by MarMauPham2016 , this dissipation is set to increase to a critical value


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


The critical damage dissipation represents the energy dissipated during a complete damage process for a volume element, and is related to the fracture toughness used in other formulations by (MarMauPham2016, ):


where the internal length is denoted by and can be expressed in terms of as


From the previous expressions and the energetic definitions (7) and (8), note that has the units of .

The model without an elastic stage has been consistently used in phase-field models, as shown in AmbGerDeL2015 ; AmbKruDel2016 and Bord-Hugh2012 ; Bord-Hugh2016 . For this choice, the absence of an elastic stage is due to . In this work, we adopt as the constitutive choice, for which . Therefore, this function allows for the existence of an initial elastic stage before the damage evolution. Moreover, the parameter 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), can be related to the fracture toughness by


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


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 and individually. Introducing the constitutive assumptions of the previous section into Eq. (8), the global energy functional reads


The alternate minimization follows.

  • Minimization with respect to the displacement field:


  • Minimization with respect to the damage field:


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


where the global vector of nodal displacements and the shape function matrix are obtained using an assembly operator :




The superscript denotes an element vector or matrix, and the index indicates the node number with representing the corresponding shape functions. The vector contains the nodal displacements in each spatial direction.

The strain vector inside an element is defined as


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




The damage and gradient damage fields are approximated by




where the assembly operator is used again:




Hereinafter, all element vectors and matrices are assumed to be assembled into their corresponding global forms through .

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


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


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


the local elastic matrix is defined in the decomposed form


where is the Heaviside step function. The global damaged stress vector can than be readily obtained from


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


from where the global displacement vector can be readily obtained after imposing boundary conditions.

Finally, the global damage field vector is obtained from the discrete version of Eq. (34), which is given by


where the element form of reads


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

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

  • Initialization: initializing convergence vectors and storage matrices.

  • One module for each primary variable: obtaining the global vectors using two separate finite-element problems.

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

  • Postprocessing: displaying the primary fields illustratively.

3.2.1 Modules for primary variables

The nodal vectors corresponding to the global primary fields, and , 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:


where is a vector containing the nodal values of a primary variable. Therefore, the solution consists of determining the local coefficient matrix and the local right hand side . Subsequently, the global coefficient matrix and the global right hand side 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 incremental displacements imposed on , are found iteratively using two independent tolerances and , for a maximum number of iterations .

1:  for  do
2:     Data input and initialization
3:     k=0
4:     while (   )     do
5:        k=k+1
6:        Obtain from module elast3D
7:        Obtain from module dam3D
8:     end while
9:  end for
Algorithm 1 General routine

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:

  • MPa;

  • ;

  • MPa mm and

  • (MPa mm) mm.

Figure 2:  General scheme for Experiment I and Experiment II.
Experiment I

The tension case is considered first, with and (Fig. 2). Vertical displacements are imposed on the top boundary from 0 to mm, with increments of 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 MieHofWel2010 , where a finer mesh of 20000 linear triangles was used.

(A)  Crack initiation close-up and crack propagation
(B) Deformed specimen
Figure 3:  Brittle crack evolution and deformed specimen for Experiment I: 2D simulation.
Figure 4:  Force-displacement curve for Experiment I: 2D simulation and comparison with the results of (MieWelHof2010, ).

For the same experiment, a three-dimensional (3D) simulation is shown in Figure 5. A uniform thickness of 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 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 liu2016 , where the implementation of the phase-field model without a damage threshold was performed in the commercial software Abaqus using 96896 hexahedral elements.

Figure 5:  Brittle crack for Experiment I: 3D simulation.
Figure 6:  Force-displacement curve for Experiment I: 3D simulation and comparison with the results of liu2016 .
Experiment II

We subject the same square specimen with a notch from Experiment I to shear loading, with and (Fig. 2). Incremental horizontal displacements are imposed on the top boundary, while the bottom boundary is fixed. Displacements are imposed from to mm, with increments of 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 MieHofWel2010 and Bord-Hugh2012 . The former used 30000 linear triangles, while the latter applied cubic T-splines with 5587 cubic basis functions.

(A) Crack initiation and crack propagation
(B) Deformed specimen
Figure 7:  Brittle crack evolution and deformed specimen for Experiment II.
Figure 8:  Force-displacement curves for Experiment II: comparison with the results of Bord-Hugh2012 .

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.


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.


  • (1) C. Comi, Computational modelling of gradient-enhanced damage in quasi-brittle materials, Mechanics of Cohesive-frictional Materials 4 (1) (1999) 17–36 (1999).
  • (2) A. A. Griffith, The phenomena of rupture and flow in solids, Philosophical transactions of the royal society of london. Series A, containing papers of a mathematical or physical character (1921) 163–198 (1921).
  • (3) G. Irwin, Analysis of stresses and strains near the end of a crack traversing a plate, Journal of Applied Mechanics 24 (1957) 361–364 (1957).
  • (4) G. Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture, in: Advances in applied mechanics, Vol. 7, Elsevier, 1962, pp. 55–129 (1962).
  • (5) G. A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, Journal of the Mechanics and Physics of Solids 46 (8) (1998) 1319–1342 (1998).
  • (6) H. Nguyen-Xuan, G. Liu, S. Bordas, S. Natarajan, T. Rabczuk, An adaptive singular es-fem for mechanics problems with singular field of arbitrary order, Computer Methods in Applied Mechanics and Engineering 253 (2013) 252–273 (2013).
  • (7) R. De Borst, H.-B. Mühlhaus, Gradient-dependent plasticity: formulation and algorithmic aspects, International Journal for Numerical Methods in Engineering 35, 521-539.(1992) (1992).
  • (8) J. Oliver, A. Huespe, E. Samaniego, E. Chaves, Continuum approach to the numerical simulation of material failure in concrete, International Journal for Numerical and Analytical Methods in Geomechanics 28 (7-8) (2004) 609–632 (2004).
  • (9) J. C. Simo, J. Oliver, F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids, Computational mechanics 12 (5) (1993) 277–296 (1993).
  • (10) N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International journal for numerical methods in engineering 46 (1) (1999) 131–150 (1999).
  • (11) N. Moës, T. Belytschko, Extended finite element method for cohesive crack growth, Engineering fracture mechanics 69 (7) (2002) 813–833 (2002).
  • (12) F. Amiri, C. Anitescu, M. Arroyo, S. P. A. Bordas, T. Rabczuk, Xlme interpolants, a seamless bridge between xfem and enriched meshless methods, Computational Mechanics 53 (1) (2014) 45–57 (2014).
  • (13) C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations, International Journal for Numerical Methods in Engineering 83 (10) (2010) 1273–1311 (2010).
  • (14) R. De Borst, C. V. Verhoosel, Gradient damage vs phase-field approaches for fracture: Similarities and differences, Computer Methods in Applied Mechanics and Engineering 312 (2016) 78–94 (2016).
  • (15) D. Mumford, J. Shah, Optimal approximations by piecewise smooth functions and associated variational problems, Communications on pure and applied mathematics 42 (5) (1989) 577–685 (1989).
  • (16) B. Bourdin, G. A. Francfort, J.-J. Marigo, Numerical experiments in revisited brittle fracture, Journal of the Mechanics and Physics of Solids 48 (4) (2000) 797–826 (2000).
  • (17) G. Dal Maso, R. Toader, A model for the quasi-static growth of brittle fractures: Existence and approximation results, Archive for Rational Mechanics and Analysis 162 (2) (2002) 101–135 (2002).
  • (18) L. Ambrosio, V. M. Tortorelli, Approximation of functional depending on jumps by elliptic functional via -convergence, Communications on Pure and Applied Mathematics 43 (8) (1990) 999–1036 (1990).
  • (19) B. Bourdin, G. A. Francfort, J.-J. Marigo, The variational approach to fracture, Journal of Elasticity 91 (1-3) (2008) 5–148 (2008).
  • (20) C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Computer Methods in Applied Mechanics and Engineering 199 (45) (2010) 2765–2778 (2010).
  • (21) M. J. Borden, C. V. Verhoosel, M. A. Scott, T. J. Hughes, C. M. Landis, A phase-field description of dynamic brittle fracture, Computer Methods in Applied Mechanics and Engineering 217 (2012) 77–95 (2012).
  • (22) M. Ambati, T. Gerasimov, L. De Lorenzis, Phase-field modeling of ductile fracture, Computational Mechanics 55 (5) (2015) 1017–1040 (2015).
  • (23) C. Miehe, S. Teichtmeister, F. Aldakheel, Phase-field modelling of ductile fracture: a variational gradient-extended plasticity-damage theory and its micromorphic regularization, Phil. Trans. R. Soc. A 374 (2066) (2016) 20150170 (2016).
  • (24) C. Miehe, F. Aldakheel, S. Teichtmeister, Phase-field modeling of ductile fracture at finite strains: A robust variational-based numerical implementation of a gradient-extended theory by micromorphic regularization, International Journal for Numerical Methods in Engineering (2017).
  • (25) M. Ambati, L. De Lorenzis, Phase-field modeling of brittle and ductile fracture in shells with isogeometric nurbs-based solid-shell elements, Computer Methods in Applied Mechanics and Engineering 312 (2016) 351–373 (2016).
  • (26) T. J. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement, Computer methods in applied mechanics and engineering 194 (39-41) (2005) 4135–4195 (2005).
  • (27) F. Amiri, D. Millán, Y. Shen, T. Rabczuk, M. Arroyo, Phase-field modeling of fracture in linear thin shells, Theoretical and Applied Fracture Mechanics 69 (2014) 102–109 (2014).
  • (28) F. Amiri, D. Millán, M. Arroyo, M. Silani, T. Rabczuk, Fourth order phase-field model for local max-ent approximants applied to crack propagation, Computer Methods in Applied Mechanics and Engineering 312 (2016) 254–275 (2016).
  • (29) H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments, Journal of the Mechanics and Physics of Solids 57 (8) (2009) 1209–1229 (2009).
  • (30) K. Pham, H. Amor, J.-J. Marigo, C. Maurini, Gradient damage models and their use to approximate brittle fracture, International Journal of Damage Mechanics 20 (4) (2011) 618–652 (2011).
  • (31) J.-J. Marigo, C. Maurini, K. Pham, An overview of the modelling of fracture by gradient damage models, Meccanica 51 (12) (2016) 3107–3128 (2016).
  • (32) R. Alessi, Variational approach to fracture mechanics with plasticity, Ph.D. thesis, Université Pierre et Marie Curie-Paris 6 (2013).
  • (33) R. Alessi, J.-J. Marigo, S. Vidoli, Gradient damage models coupled with plasticity and nucleation of cohesive cracks, Archive for Rational Mechanics and Analysis 214 (2) (2014) 575–615 (2014).
  • (34) R. Alessi, J.-J. Marigo, C. Maurini, S. Vidoli, Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: one-dimensional examples, International Journal of Mechanical Sciences (2017).
  • (35) J. Ulloa, P. Rodríguez, E. Samaniego, On the modeling of dissipative mechanisms in a ductile softening bar, Journal of Mechanics of Materials and Structures 11 (4) (2016) 463–490 (2016).
  • (36) P. Rodriguez, J. Ulloa, C. Samaniego, E. Samaniego, A variational approach to the phase field modeling of brittle and ductile fracture, International Journal of Mechanical Sciences 144 (2018) 517–502 (2018).
  • (37) A. Mielke, A mathematical framework for generalized standard materials in the rate-independent case, Multifield Problems in Solid and Fluid Mechanics 28 (2006) 399 (2006).
  • (38) A. Mielke, T. Roubícek, Rate-Independent systems. Theory and application, Springer, 2015 (2015).
  • (39) R. Alessi, Energetic formulation for rate-independent processes: remarks on discontinuous evolutions with a simple example, Acta Mechanica 227 (10) (2016) 2805–2829 (2016).
  • (40) B. Halphen, Q. Nguyen, Generalized standard materials, Journal de mécanique 14 (1) (1975) 39–63 (1975).
  • (41) P. Germain, Q. S. Nguyen, P. Suquet, Continuum thermodynamics, Journal of Applied Mechanics 50 (4b) (1983) 1010–1020 (1983).
  • (42) K. Pham, J.-J. Marigo, Approche variationnelle de l’endommagement: I. les concepts fondamentaux, Comptes Rendus Mécanique 338 (4) (2010) 191–198 (2010).
  • (43) R. Faria, J. Oliver, M. Cervera, Modeling material failure in concrete structures under cyclic actions, Journal of Structural Engineering 130 (12) (2004) 1997–2005 (2004).
  • (44) F. Freddi, G. Royer-Carfagni, Regularized variational theories of fracture: a unified approach, Journal of the Mechanics and Physics of Solids 58 (8) (2010) 1154–1174 (2010).
  • (45) C. Steinke, M. Kaliske, A phase-field crack model based on directional stress decomposition, Computational Mechanics (2018) 1–28 (2018).
  • (46) M. J. Borden, T. J. Hughes, C. M. Landis, A. Anvari, I. J. Lee, 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 312 (2016) 130–166 (2016).
  • (47) M. Ambati, R. Kruse, L. De Lorenzis, A phase-field model for ductile fracture at finite strains and its experimental verification, Computational Mechanics 57 (1) (2016) 149–167 (2016).
  • (48) J. Alberty, C. Carstensen, S. A. Funken, Remarks around 50 lines of matlab: short finite element implementation, Numerical Algorithms 20 (2-3) (1999) 117–137 (1999).
  • (49) G. Liu, Q. Li, M. A. Msekh, Z. Zuo, Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field model, Computational Materials Science 121 (2016) 35–47 (2016).