User login


You are here

Journal Club For April 2021: Variational phase-field modeling of brittle and cohesive fracture

phunguyen's picture

Journal Club For April 2021: Variational phase-field modeling of brittle and cohesive fracture

Vinh Phu Nguyen (Monash University) and Jian-Ying Wu (South China University of Technology)

1. Griffith’s linear elastic fracture mechanics

Without being historically exact one can say that fracture mechanics has started with the work of Alan Arnold Griffith (1893–1963). As Griffith’s work is the basic of variational phase-field fracture model developed 80 years later, let’s examine what he did.

In 1921, Griffith conducted experiments on fracture of glass fibers [1]. He found two things: (1) the fracture strength of glass is significantly smaller than the theoretical value (coming from breaking the atomic bonds) and (2) small glass fibers are stronger than larger fibers. He concluded that small naturally occurring defects existing in the fibers make them weak. Precisely these defects amplify the stress field in front of their tips and thus rendering the fracture stress much smaller compared with the theoretical strength. Actually Griffith was aware of the work of Charles Inglis (1875–1952) conducted seven years ago about stress concentration due to an elliptical hole [2].

With defects now in his mind, he made specimens with artificial surface cracks (to overcome natural defects) of varying sizes and quantified the relationship between the remote tensile stress  and the crack size or length a. What he found is an inverse relation between the strength and the crack length. He found that the product of the applied stress with the square root of the crack length  is a constant. To find this constant, he carried out an energy-based analysis that basically led to the born of what is now known as fracture mechanics.

Griffith computed the energy of the system which consists of the stored elastic strain energy and the surface energy i.e., energy due to the creation of the new crack surfaces. He considered a unit thickness infinite plate with a surface crack of length a subjected to a remote tensile stress σ  normal to the crack. The energy of this system is given by

where he used Inglis’ solution to obtain the elastic strain energy released due to the crack’s presence; U0 is the elastic strain energy of the plate without crack. What is interesting here is the last term -- the surface energy associated with a crack of length a where gamma_s is the energy required to creare a unit surface area.

The first derivative of U with respect to the crack length a is

And the vanishing derivative condition gives us

which is the well-known Griffith’s equation relating the remote stress to the crack length.

The work of Griffith which is applicable only to brittle materials (e.g. glasses) was ignored for almost 20 years. It was not until the modifications made by Orowan and particularly George Rankin Irwin (1907–1998) that, a new field has emerged: Linear Elastic Fracture Mechanics (LEFM). Irwin introduced the concept of energy release rate G which is the negative of the derivative of the elastic strain energy with respect to the crack length and he and Orowan replaced the surface energy by Gf — the critical energy release rate or fracture energy — to take into account other dissipative processes such as plastic deformation, i.e.,

Another significant contribution to fracture mechanics was made by James Rice (1940) in 1968, the famous J -integral [3] which is equal to the fracture energy release rate G.

By recalling the irreversibility of the crack propagation that the time derivative of the crack length is non-negative, the Griffith crack propagation criterion is given then:

That is, for quasi-static loading case, the crack propagates if G > Gf and otherwise remains stationary for G < Gf.

Fracture mechanics has been a great success as it provides the engineers a continuum mechanics tool to quantita- tively predict the structural integrity of large structures using data such as fracture toughness (a concept introduced by Irwin which is related to the fracture energy) which can be experimentally measured using laboratory scale specimens. Furthermore it helps material scientists to improve existing materials and design new ones by looking at their fracture toughness.

2. Barenblatt’s cohesive zone model

In 1959, Grigory Isaakovich Barentblatt (1927–2018) [4] proposed the celebrating cohesive zone model (CZM). Barentblatt’s CZM solved two major issues of classical fracture mechanics: crack nucleation and stress singularity at the crack tip. Its idea is to lump what is going on in the fracture process zone (FPZ) ahead of a traction-free crack tip on a surface, and approximate this zone by a traction-separation law that relates the cohesive traction transfered across the crack surfaces σ and the crack opening w. In this case, the fracture energy release rate G(w) is no longer a constant but a non-convex function of the crack opening w as shown in Figure 1.

Thus, instead of instant dissipation upon the creation of a unit fracture surface in Griffith’s theory, in Barenblatt’s CZM, the energy is released gradually:

The traction corresponding to zero crack opening is the tensile strength ft of the material. It is then able to define a material characteristic length


which is called Irwin’s internal length [6], with E0 being Young’s modulus of the material. It characterizes (is proportional to) the size of the FPZ. Hence, it also measures the brittleness of the material: the smaller (compared to the structural size) it is, the more brittle the material behaves. Note that the CZM applies not only to cohesive cracks but also to brittle ones with a small Irwin’s length.

It was A. Hillerborg et al. [7] who assumed that the cohesive crack may develop anywhere, even if no a priori existing macrocrack is actually present and named this extension as the fictitious crack model. They also implemented the CZM in a finite element framework to model fracture of concrete beams [7]. Another notable work is of Xu and Needleman [8] for dynamic fracture. CZM is usually implemented using the so-called interface elements. Recent works of Paulino’s group and Papoulia’s demonstrate that CZM and advanced interface elements are indeed a powerful tool for fracture simulations [9, 10].

3. The variational approach to fracture

In 1998 i.e., 78 years after Griffith’s work Francfort and Marigo reformulated Griffith’s energetic theory in a variational framework and coined it the variational approach to brittle fracture [11]. This approach generalizes LEFM by allowing crack nucleation and arbitrary crack propagation within a single framework.

In this energetic approach, crack propagation results from the competition between the bulk energy away from the crack and the surface energy on the crack. From this point of view, the total energy functional E in a quasi-static loading regime reads

for the external potential energy P.

Motivated by Griffith theory, Francfort and Margio [11] recast brittle fracture as an energy minimization problem:

Note that for the case in which the crack path  may be not smooth “enough”, an infimum energy rather than a minimal one should be sought. In such a variational approach to brittle fracture, cracks should propagate along the path of least energy. In particular, for an a priori existing crack constrained to propagate along a pre-defined crack path, Griffith’s criterion is exactly retrieved [12]. However, the most significant merit is that with this variational approach to fracture it is possible to deal with crack nucleation in an initially perfectly sound solid and to determine intrinsically crack paths in a variationally consistent manner, bypassing the underlying assumptions of LEFM.

4. The variational phase-field model for fracture

Though there exist other schemes, e.g., the eigenerrosion method [13, 14], the adaptive finite elements, to approxi- mate the variational approach to brittle fracture, the variational phase-field model for fracture, set forth in [15], might be the most versatile. This numerically more amenable counterpart of the Francfort–Marigo’s variational approach to brittle fracture is motivated by the Ambrosio–Tortorelli elliptic regularization [16] of the Mumford–Shah functional [17] in image segmentation problems. The idea is to introduce a continuous scalar field — the crack phase-field or damage field d(x) — that takes either a value of zero for intact material, or a value of unity for completely broken material or a value between 0 and 1 for partially broken material. This field helps to approximate the Griffith surface energy–a surface integral–as a volume integral over the computational domain Ω. So, in this regularized model, for quasi-static fracture of solids under the infinitesimal strain regime, the displacement field u and damage field d are minimizers of the following total energy functional of the solid

where the first integral is the stored strain energy influenced now by the crack phase-field, the second one denotes the fracture energy a` la Griffith. The phase-field length scale b is a regularization parameter that controls the crack band i.e., where d(x) is non-zero; ω(d) is the energetic degradation function and α(d) represents the crack geometric function. This positive/negative decomposition plays an important role in capturing the tension/compression asymmetry of fracture and in removing spurious compressive fracture – fracture does not occur in domains under compression; see [18] for the discussion of several positive/negative decomposition schemes of the elastic strain energy.

Where does the name phase field come from? In the similar phase-field models for phase transformation based upon the Ginzburg–Landau equation [19] in physics and materials science, for multiphase fluid based upon the Cahn–Hillard [20] or Allan–Cahn [21] equation in fluid mechanics, the so-called order parameter is also introduced to discriminate the distinct states or phases. Here, the crack phase-field d(x) bears the same role in discriminating two phases of the solid: the intact phase and the completely broken phase, Kuhn and R. Mu ̈ller called this model a phase-field fracture model (PFM) [22] in 2010, ten years after the work of Bourdin, Francfort and Marigo! To highlight its variational argument, we believe a better term is variational phase-field fracture model. In the same year, the later mechanician, Christian Miehe, published two papers [23, 24] presenting a more intuitive approach alternative to the formal and mathematically demanding formulation of Francfort–Marigo–Bourdin paper. Indeed Miehe was successful in making PFMs much more attractive in the engineering community: we have seen a surge in publications on PFM since 2010. Another key player is probably Thomas Hughes who proposed the fourth-order PFM [25] and presented many keynote lectures on applications of PFM to dynamic fracture and ductile fracture at various Complas (Computational Plasticity) conferences.

Looking at the expression of the free energy functional one can see that there are many possibilities in choosing the degradation function  and the crack geometric function. We refer to [18] for an extensive discussion. The three most common (second order) PFMs are listed in Table 1. The AT2 model, developed in [15, 23] is probably the most widely used PFM in engineering even though it lacks an elastic domain i.e., damage becomes non-zero immediately when the load (no matter how small it is) is applied. The AT1 model of [26], possessing an elastic domain, is getting more attention. While both the AT1 and AT2 models apply only to brittle fracture, the PF-CZM (phase field cohesive zone model) of [27] is the first PFM that applies to both brittle and cohesive fracture.

Herein we need to stress the most important difference between AT1/2 and PF-CZM. That is, the latter converges to the CZM while all the others to the Griffith LEFM, in the Gamma-convergence [28], to be discussed shorly. 

5. Role of the length scale and crack nucleation

Using the tools from free discontinuity problems, it can be shown that when the phase-field length scale b approaches zero in the vanishing limit, the PFM solution converges to the solution of the original problem. This is known as Gamma-convergence [28]. Therefore, in the early days of PFMs the length scale was considered as a merely numerical parameter. However during verification attempts against experiment data, it was soon realized that for the AT1/2 models the length scale must be considered as a material parameter relating to Irwin’s internal length lch so that the failure strengthft can be matched in a uni axial traction state [26,30]. Due to this new interpretation of the length scale, Marigo re-named it as a gradient damage model [26] since it belong to the gradient damage model developed by Fre ́mond and Nedjar in 1996 [31]. Still, as recently pointed out in [32] such models cannot capture crack nucleation in intact solids with no pre-existing defects. To demonstrate this issue, we model the uniaxial tensile test of a biological tissue [33], and we used different values for the length scale. To match the experiment a large length scale have to be used, but this value results in an erroneously wide crack band (Figure 2). 

For the PF-CZM, on the other hand, the phase-field length scale is merely a numerical parameter (and it can also be treated as a physical one related to possibly the microscopic inhomogeneity of the material) so long it is small enough compared to the structural size. As seen from Figure 2c, all length scales yield the same load-displacement responses and the peak is similar to the experiment. Therefore, we can use a small length scale to guarantee the Gamma-convergence. See [33] for details regarding the formulation for large deformation anisotropic fracture of hyperelastic solids.

We present another example reported by Bourdin himself [34] so that our discussion would be more convincing. This example is the popular L-shape panel test of concrete fracture. The specimen size is about 500 mm. Using a small length scale (b=3.125 mm), Bourdin got a crack path quite matching the experiment, but the peak load is highly overestimated. To match the failure strength of concrete, a much larger phase-field length scale b=126.32 mm has to be used. With this length scale, a better load–displacement curve is obtained, but the crack path is completely polluted by the too diffuse regularization (Figure 3a). Comparatively, this is not an issue for the PF-CZM at all – both the crack path and the load – displacement curve can be well captured, both insensitive to the phase-field length scale as shown in Figure 3b. This is because the PF-CZM incorporates the stress-based crack nucleation criterion, the energy-based crack propagation criterion and the variational principle based crack path chooser into one standalone framework.

6. Numerical implementation of variational phase-field fracture models

Compared with discrete fracture numerical methods such as XFEM [35], the implementation of a PFM is quite simple. We do not need to introduce the cumbersome crack tracking algorithm nor any ad hoc criteria for determining the crack orientation/branching. Only standard continuum finite elements with both displacement and damage degrees of freedom (dofs) are needed in the numerical implementation of a PFM. Compared to those classical continuum damage models in which particular strategies (e.g., the mixed stabilized finite elements recently advocated by Cervera and coworkers [36]) have to be employed, the solution of PFMs do not depend on the underlying spatial discretization. That is, just as elastostatic problems, almost all the existing numerical methods, e.g., irreducible and mixed FEs, meshless or meshfree methods, particle methods (MPM, RKPM) and even the more recent peridynamics, can be employed. There is no stability issue nor mesh bias dependence [37]. Moreover, for the PF-CZM the solution is also insensitive to the phase-field length scale. In other words, users need only to focus on the physical problems of interest rather than worry about those numerical issues/parameters.

As other coupled multi-physics problems, the solver used to solve these equations is more involved as the monolithic solver using the standard Newton-Raphson method usually does not work, since the energy functional E is non-convex with respect to (u,d). Therefore, the most common solver is the alternate minimization or staggered (AM/staggered) solver in which one first fixes the damage dof and solve for the displacement dof, and then solves for the damage dof using the updated displacement dof. Though the one-pass AM/staggered solver has been widely adopted, it would result in delay of damage evolution, inaccurate post-peak regimes and spurious damage widening when the crack arrives at the external boundary. These issues are largely alleviated in the iterative multi-pass AM/staggered solver.

Staggered solvers are easy to implement and stable but their convergence is extremely slow. Among many mono- lithic solvers proposed, the BFGS solver — a popular solver in nonlinear optimization problems — has proven to be promising: the BFGS-based monolithic solver with the most crude initial guess of the system matrix (BFGS-SM0) is about 4 to 8 times faster than the AM/staggered solver [44], and with the optimized system modification (BFGS-SM1) the computational efficiency can be further boosted [45]; see Figure 4 for the comparison. We have also demonstrated that this gain in computational efficiency also applies to thermo-mechanical fracture problems [46].

We note again that PFMs can be implemented in any numerical discretization method and the PF-CZM is insensitive to those numerical issues/parameters. The only requirement is that the phase-field length is sufficiently resolved by the spatial discretization.

7. Applications

Variational phase-field fracture models have been applied to many problems. And in [18] a quite comprehensive review has been provided, so we herein present some of our recent works covering static and dynamic fracture of both brittle and quasi-brittle materials.

7.1. Quasi-static brittle fracture

Figure 5 presents a PF-CZM simulation of the peeling test of a biological tissue. It is taken from our work in [33]. See [47, 48] for other works on finite deformation fracture of rubber-like materials using phase-field models. 

7.2. Dynamic brittle fracture

PFMs have been applied to dynamic fracture as reported in [49, 50, 51, 52, 53, 54, 55]. The idea is simple: the rigorously derived and well studied variational approach to quasi-static fracture is maintained and inertial effects are incorporated. Borden et al. used the Hamilton’s principle of least action to derive the governing equations [52].


From the results reported in the literature, and particularly from our work in [57], we can say that:

* Phase-field fracture simulation results are very encouraging. Indeed, PFMs can capture several dynamic fracture

phenomena: crack branching, crack arrest, fragmentation and multiple branching (Figure 6);

 * Phase-field simulations results are similar to predictions of peridynamics and discontinuous Galerkin extrinsic cohesive elements [57];

* The AT1 model and PF-CZM are spatial and temporal convergent for dynamic brittle fracture problems;

* Only a few quantitative assessment of dynamic brittle fracture has been done. More work are needed.

7.3. Quasi-static cohesive fracture

Only when it comes to 3D fracture problems with non-planar crack surfaces that the strength of PFMs is obvious. One can capture complex 3D crack paths involving merging, branching, twisting with a relatively simple implementa- tion (quite the same implementation as a 2D implementation, precisely) and the entire fracture process is done over a fixed finite element mesh. To illustrate this, we show in Figure 8  one simulation of mixed-mode I/III fracture taken from our recent work [58].

8. Conclusions

It has been a century since Griffith’s seminal work on brittle fracture. Let’s us summarize this adventure. Griffith set up an energy principle for brittle crack propagation. Irwin developed an equivalent approach with the concept of stress intensity factor and established the framework of linear elastic fracture mechanics. Barenblatt proposed the well-known cohesive zone model and greatly broadened the scope of fracture mechanics to nonlinear territories. Rice proposed the J -integral and proved its relation to the fracture energy release rate.

About eighty years after Griffith, with the help of modern variational calculus, Francfort and Marigo were able to cast Griffith’s energetic principle into a variational framework. This variational approach to fracture generalizes Griffith’s energetic principle by allowing crack nucleation and non-predefined crack paths without resorting to external criteria. Then comes Bourdin’s phase-field approximation of the variational approach. For simple problem, Gamma-comvergence says that when length scale is going to zero, Bourdin’s approximate model converges to the original problem with sharp cracks. Our recent work on the phase-field cohesive zone model then sheds light on the promising prospect in unifying fracture mechanics and continuum damage mechanics into a single framework.

Variational phase-field fracture models have been applied to more and more problems. We can cite multi-physics fracture (hydraulic fracture, hydrogen-assisted cracking, thermo-elastic fracture etc.), ductile fracture, fatigue, anisotropic fracture e.g. fracture of fiber reinforced composites, biological tissues. However, to be a practical tool one needs to lower down the computational cost of these models.

For a better reading experience, we also attach the pdf of this post (


[1] A. A. Griffith. The phenomena of rupture and flow in solids. Philosophical Transactions of the Royal Society of Londres, 221:163–198, 1920.

[2] C. E. Inglis. Stresses in plates due to the presence of cracks and sharp corners. Transactions of the Institute of Naval Architects, 55:219–241, 1913.

[3] J. R. Rice. A path independent integral and the approximate analysis of strain cncentrations by notches and cracks. J. Appl. Mech.-T. ASME, 35:379–386, 1968.

[4] G.I. Barenblatt. The formation of equilibrium cracks during brittle fracture. general ideas and hypotheses. axially-symmetric cracks. Journal of Applied Mathematics and Mechanics, 23:622–636, 1959.

[5] A. A. Griffith. The phenomena of rupture and flow in solids. Philosophical Transactions of the Royal Society of Londres, 221:163–198, 1920.

[6] G. R. Irwin. Fracture. In W. Fl”ugge, editor, Handbuch der Physik, Vol. VI, pages 551–590, Springer-Verlag, Berlin, 1958.

[7] A. Hillerborg, M. Mod´eer, and P.E. Petersson. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement Concr. Res., 6:773–781, 1976.

[8] X.P. Xu and A. Needleman. Numerical simulations of fast crack growth in brittle solids. Journal of the Mechanics and Physics of Solids, 42(9), 1994.

[9] D. W. Spring and G. H. Paulino. Achieving pervasive fracture and fragmentation in three-dimensions: an unstructuring-based approach. International Journal of Fracture, 210(1):113–136, Mar 2018.

[10] M Reza Hirmand and Katerina D Papoulia. Block coordinate descent energy minimization for dynamic cohesive fracture. Computer Methods in Applied Mechanics and Engineering, 354:663–688, 2019.

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

[12] B. Bourdin, G. Francfort, and J.-J. Marigo. The variational approach to fracture. Springer, Berlin, 2008.

[13] B. Schmidt, F. Fraternali, and M. Ortiz. Eigenfracture: An eigendeformation approach to variational fracture. Multiscale Modeling & Simulation, 7(3):1237–1266, 2009.

[14] A. Pandolfi and M. Ortiz. An eigenerosion approach to brittle fracture. International Journal for Numerical Methods in Engineering, 92(8):694–714, 2012.

[15] B. Bourdin, G.A. Francfort, and J.-J. Marigo. Numerical experiments in revisited brittle fracture. J. Mech. Phys. Solids, 48(4):797–826, 2000.

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

[17] D. Mumford and J. Shah. Optimal approximations by piecewise smooth functions and associated variational problems. Communications on Pure and Applied Mathematics, 42(5):577–685, 1989.

.[18] JY Wu, VP Nguyen, C. Nguyen-Thanh, D. Sutula, S. Sinaie, and S Bordas. Phase-field modeling of fracture. In Mahmoud I. Hussein, editor, Fracture Mechanics: Recent Developments and Trends, volume 53 of Advances in Applied Mechanics, pages 1 –x. Elsevier, 2019.

[19] L.D. Landau and E.M. Lifshitz. Statistical physics. Pergamon Press, Oxford, 1980.

[20] J. W. Cahn and J. E. Hilliard. Free energy of a nonuniform system. i. interfacial free energy. J. Chem. Phys., 28:258, 1958.

[21] S.M. Allen and J.W. Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metall., 27:1085–1095, 1979.

[22] C. Kuhn and R. M¨uller. A continuum phase field model for fracture. Engineering Fracture Mechanics, 77(18):3625 – 3634, 2010.

[23] C. Miehe, M. Hofacker, and 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-48):2765 – 2778, 2010.

[24] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multifield fe implementations. Int. J. Numer. Meth. Engng., 83:1273–1311, 2010.

[25] M. J. Borden, T. J.R. Hughes, C. M. Landis, A. Anvari, and 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:130 – 166, 2016.

[26] K. Pham, H. Amor, J.-J. Marigo, and C. Maurini. Gradient damage models and their use to approximate brittle fracture. International Journal of Damage Mechanics, 20:618–652, 2011.

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

[28] A. Braides. Approximation of Free-Discontinuity Problems. Springer science & Business Media, Berlin, 1998.

[29] Kyoungsoo Park and Glaucio H Paulino. Cohesive zone models: a critical review of traction-separation relationships across fracture surfaces. Applied Mechanics Reviews, 64(6):060802, 2011.

[30] T. T. Nguyen, J. Yvonnet, M. Bornert, C. Chateau, K. Sab, R. Romani, and R. Le Roy. On the choice of parameters in the phase field method for simulating crack initiation with experimental validation. International Journal of Fracture, 197(2):213–226, 2016.

[31] M. Fr´emond and B. Nedjar. Damage, gradient of damage and principle of virtual power. International Journal of Solids and Structures, 33(8):1083 – 1103, 1996.

[32] Aditya Kumar, Blaise Bourdin, Gilles A Francfort, and Oscar Lopez-Pamies. Revisiting nucleation in the phase-field approach to brittle fracture. Journal of the Mechanics and Physics in Solids, 2020.

[33] T. K. Mandal, V. P. Nguyen, and J.-Y. Wu. A length scale insensitive anisotropic phase field fracture model for hyperelastic composites. International Journal of Mechanical Sciences, (105941), 2020.

[34] A. Mesgarnejad, B. Bourdin, and M.M. Khonsari. Validation simulations for the variational approach to fracture. Computer Methods in Applied Mechanics and Engineering, 290:420 – 437, 2015.

[35] N. Mo¨es, J. Dolbow, and T. Belytschko. A finite element method for crack growth without remeshing. Int. J. Numer. Meth. Engng., 46:131–150, 1999.

[36] M. Cervera, M. Chiumenti, and R. Codina. Mesh objective modeling of cracks using continuous linear strain and displacement interpolations. Int. J. Numer. Meth. Engng., 87(10):962–987, 2011.

[37] T K Mandal, V P Nguyen, and J-Y Wu. Length scale and mesh bias sensitivity of phase-field models for brittle and cohesive fracture. Engineering Fracture Mechanics, 217(106532), 2019.

[38] S. J. Benson and T. S. Munson. Flexible complementarity solvers for large-scale applications. Optimization Methods and Software, 21:155–168, 2006.

[39] H. Amor, J.J. Marigo, and C. Maurini. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments.J. Mech. Phys. Solids, 57:1209–1229, 2009.

[40] P.E. Farrell and C. Maurini. Linear and nonlinear solvers for variational phase-field models of brittle fracture. Int. J. Numer. Meth. Engng., 109(5):648–667, 2017.

[41] J. Y. Wu and V. P. Nguyen. A variationally consistent anisotropic phase-field damage model for fracture. Journal of the Mechanics and Physics of Solids, page submitted, 2018.

[42] Mary F. Wheeler, Thomas Wick, and W. Wollner. An augmented-lagrangian method for the phase-field approach for pressurized fractures.Comput. Methods Appl. Mech. Engrg., 271:69–85, 2014.

[43] T. Gerasimov and L. De Lorenzis. On penalization in variational phase-field models of brittle fracture. Computer Methods in Applied Mechanics and Engineering, 354:990–1026, 2019.

[44] J. Y. Wu, Y. Huang, and V. P. Nguyen. On the bfgs monolithic algorithm for the unified phase field damage theory. Comput. Methods Appl. Mech. Engrg., page 112704, 2020.

[45] J. Y. Wu, Y. Huang, H. Zhou, and V. P. Nguyen. Three-dimensional phase-field modeling of mode i + ii/iii failure in solids. Computer Methods in Applied Mechanics and Engineering, 373.

[46] T. K. Mandal, V. P. Nguyen, J.Y. Wu, C. Nguyen-Thanh, and A. de Vaucorbeil. Fracture of thermo-elastic solids: Phasefield modeling and new results with an efficient monolithic solver. Computer Methods in Applied Mechanics and Engineering, 376:113648, 2021.

[47] C. Miehe and L.-M. Sch¨anzel. Phase field modeling of fracture in rubbery polymers. part i: Finite elasticity coupled with brittle failure. Journal of the Mechanics and Physics of Solids, 65:93 – 113, 2014.

[48] O. G¨ultekin, H. Dal, and G. A. Holzapfel. A phase-field approach to model fracture of arterial walls: Theory and finite element analysis. Computer Methods in Applied Mechanics and Engineering, 312:542 – 566, 2016.

[49] C. J. Larsen, C. Ortner, and E. Sali. Existence of solutions to a regularized model of dynamic fracture. Mathematical Models and Methods in Applied Sciences, 20(07):1021–1048, 2010.

[50] B. Bourdin, C. J. Larsen, and C. L. Richardson. A time-discrete model for dynamic fracture based on crack regularization. International Journal of Fracture, 168(2):133–143, 2011.

[51] M. Hofacker and C. Miehe. Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation. International Journal of Fracture, 178(1):113–129, 2012.

[52] M. J. Borden, C. V. Verhoosel, M. A. Scott, T.J.R. Hughes, and C. M. Landis. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering, 217-220:77 – 95, 2012.

[53] A. Schl¨uter, A. Willenb¨ucher, C. Kuhn, and R. M¨uller. Phase field approximation of dynamic brittle fracture. Computational Mechanics, 54(5):1141–1161, 2014.

[54] Tianyi Li, Jean-Jacques Marigo, Daniel Guilbaud, and Serguei Potapov. Gradient damage modeling of brittle fracture in an explicit dynamics context. Int. J. Numer. Meth. Engng., 108(11):1381–1405, 2016.

[55] J. Bleyer, C. Roux-Langlois, and J.-F. Molinari. Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms. International Journal of Fracture, 204:79–100, 2017.

[56] Kazuo Arakawa and Kiyoshi Takahashi. Branching of a fast crack in polymers. International Journal of fracture, 48(4):245–254, 1991.

[57] T. K. Mandal, V. P. Nguyen, and J-Y Wu. Evaluation of variational phase-field models for dynamic brittle fracture. Engineering Fracture Mechanics, 235(107169), 2020.

[58] Jian-Ying Wu, Yuli Huang, Hao Zhou, and Vinh Phu Nguyen. Three-dimensional phase-field modeling of mode I + II/III failure in solids. Computer Methods in Applied Mechanics and Engineering, 373:113537, 2021.



Konstantin Volokh's picture

 Dear Dr. Nguyen,

This topic is interesting. Fig. 2 is not clear. You say you consider homogeneous tension. Then, how localization of damage can occur? You should consider truly homogeneous case in which failure is also homogeneous. This will give you strength.

Best regards,


phunguyen's picture

Dear Dr Volokh,

Thanks for your comment. According to Ref. [26], K. Pham, H. Amor, J.-J. Marigo, and C. Maurini. Gradient damage models and their use to approximate brittle fracture. International Journal of Damage Mechanics, 20:618–652, 2011, who used a second-order stability condition to show that the softening stage of the homogeneous solution is stable only if the bar is short enough. 

So it is possible to have localized damage even for a homogenous state.

Best regards,


Konstantin Volokh's picture

Dear Phu,

You consider uniform tension. The deformation should be uniform until the limit point on the stress-strain curve. The limit point is called strength. Strength in your theory (or any phase field formulation) will depend on the characteristic length (b).


Jian-Ying Wu's picture

Dear Kosta,

Thank you for your interest to our work. PF-CZM is insensititve to the length scale parameter (b) with the failure strength fixed as a material constant. This is different from other phase-field models. Please check our previous papers [27, 41] for more details on how this property is achieved.

Under uniform tension, though the position might be random, localization is indeed possible and there is no need to introduce any predefined defect. 


Konstantin Volokh's picture

Dear Jiang-Ying,

Consider your model before damage localization. For this purpose, you should assume that the gradient of the phase field is zero. Under this assumption, the differential equation for the phase field reduces to algebraic equation. You can solve it for the phase field. The phase field will depend on the characteristic length. Then, substitute the found phase field into the constituive law for stresses and you will get the dependence of the constitutive law on the characteristic length. The latter also means that both stiffness and strength depend on the characteristic length. Am I missing something?



phunguyen's picture

Hi Kosta,

Your analysis is correct for the original PFM of Bourdin et all, dubbed AT2 in our post. But in AT1 and PF-CZM, the phase field is identically zero before damage localisation. 



Konstantin Volokh's picture

Dear Phu,

My analysis is correct for any model :))

The phase field must be a solution of the phase field equation. If the solution is identically zero, then you cannot initiate the localization at all.


phunguyen's picture

Hi Kosta,

See this for details.

In short, PF-CZM is similar to CZM: initially the bar is linear elastic, the stress is homogeneous and its magnitude is smaller than the material strength f_t. This is the case until the stress equals f_t, a cohesive crack is initiated, and the stress is decreasing. 

PF-CZM is just a geometrically regularised CZM, in which the phase-field is used to represent the sharp cohesive cracks. 

Hope it helped.



Konstantin Volokh's picture

Sorry, I do not understand transition from (9.2) to (9.3). You should solve (9.2)_2 for d: Q(d)=0. You cannot assume d=0 in advance!

Jian-Ying Wu's picture

yes, d= 0 holds for AT1 and PF-CZM, but not for AT2. For the AT2 the damage is actiavated immediately once the load is applied no matter how small it is.

The damage is not activated until sigma = ft for AT1 and PF-CZM, but in the AT1 sigma is dependent on the length scale while the PF-CZM. I suggest you to read our AAM paper for the differences between AT1/AT2/PF-CZM.

Konstantin Volokh's picture

Sorry again. d=0 cannot "hold". It must be solution of Q(d)=0 in any theory. If that is not the case then you do not have theory - you have a tricky algorithm instead.

Jian-Ying Wu's picture

Note that d only needs to satisfy Q <=0, not necessarily Q = 0. So d = 0 is the solution. No trick at all.

Konstantin Volokh's picture

Then, there is no instability and strength of material.

Jian-Ying Wu's picture

yes, there is when sigma = ft for AT1/PF-CZM.

Konstantin Volokh's picture

You impose this strength as an extra constraint, and it does not come out of the constitutive law as in damage mechanics. No material instability.

Konstantin Volokh's picture

Well, if it is CZM then you cannot create stress-strain curve in uniform tension without imperfection. Something should break symmetry...

Blaise Bourdin's picture



I did not read the thread in much detail, but I think that the confusion comes from the fact that the variational phase field theory is variational by nature (unless its variational nature has been destroyed by some ill-advised modification). Because the total energy is non-convex, the first order optimality conditions: stability with respect to displacement, i.e. static equilibrium, and stability with respect to the phase field variable, i.e. a "smoothed" version of $G=G_c$) are onle _necessary_ conditions. The long series of papers initiated by the work of Pham, Marigo and Maurini in 2010 studied higher order optimality conditions to highlight the fact that above a given loading, the solution with no damage d=0 or uniform damage (d = cst) is no longer a stable critical point of the phase field energy, and that _full_ localisation (which can be interpreted as crack nucleation) of the damage variable takes place.


Konstantin Volokh's picture

Thanks. Do you mean that material becomes unstable even without damage (d=0)?

Blaise Bourdin's picture

I mean that for the AT1 model, the one dimensional response consists of an elastic phase (d=0) until a critical stress \sigma_c. Past sigma_c, a solution with d=cst (corresponding to negative rademing) satisfying the _first_ order optimality conditiuon can be constructed but can be shown to be "unstable" in teh sense that it does not satisfy second order optimality conditions. Further analysis sows that necessarily, in this condition, the maximum value of d must be 1, and a _stable_ solution (corresponding to the classical AMbrosio Tortorelli optimal profile for Gamma convergence) can be constructed.

For the AT2 model, the situation is slightly different as the elastic limit is 0. The soltion with homogeneous damage d=cst is stable in the positive hardening phase, but as before unstable in the negative hardening phase. Again, biffurcation towards fully localized state takes place.

I think that the exposition in Marigo, J.-J., Maurini, C., and Pham, K. (2016). An overview of the modelling of fracture by gradient damage models. Meccanica, 51(12):3107–3128. is quite clear (although technical).


Konstantin Volokh's picture

Let me ask it simpler. Can you create the stress-strain curve with a limit point (indicating strength) without introducing the strength explicitly as an additional constraint?

Blaise Bourdin's picture


Yes, provided that the regularization parameter \ell be chosen as sigma_c  = \sqrt{3G_c E/8\ell} for the AT1 model.

The essence of [Tanné et al 2018] is to illustrate the dual nature of the phase-field approach: an approximation of Griffith criterion as \ell approaches 0, which inherits from botyh Griffith's strenght and weakness, and a gradient damage model properly accounting for strength and toughness (and the transition from one to the other) when \ell is fixed, given by the relation above.  Note that the standard expression for the cohesive length is \ell_c = K_Ic^2/sigma_c^2, so up to a rescaling, the regularization parameter \ell is nothing but the standard cohesive legth, and the Griffith approximation makes sense when the size of teh structire is very large compared to this lebgth.


Konstantin Volokh's picture

Well, then the characterisitc length must be fixed and it is not a varying parameter anymore. No physically meaningful "convergence" under the decrease of the characteristic length will take place. Various lengths will trigger different physics.

Blaise Bourdin's picture

Well, it is a bit more complicated than this... Although Gamma convergence to teh Francfort-Marigo generalized Griffith energy requires that \ell converges to 0, Sicsic, P. and Marigo, J.-J. (2013). From gradient damage laws to Griffith’s theory of crack propagation. J. Elasticity, 113(1):55–74. showed that a propagating crack satisfies G=G_c for any \ell  "small" compared to the domain size (back to the idea that Griffith's approximation makes sense at a scale larger than the cohesize length). This can actually be seen in Figure 1 of Tanné, E., Li, T., Bourdin, B., Marigo, J.-J., and Maurini, C. (2018). Crack nucleation in variational phase-field models of brittle fracture. J. Mech. Phys. Solids, 110:80–99.

Konstantin Volokh's picture

Griffith should not bother you. His theory was pioneering and I appreciate that, yet it contradicts experiments. Cracks are not "mathematical" - they have finite thickness and the characteriistic length should not go to zero.

Jian-Ying Wu's picture

In the pf-czm, the length scale can also be fixed to match the damage bandwidth (provided it is available), without sacrifying the critical strength.

Konstantin Volokh's picture

I have to think more about your model. It is not completely clear to me. CZM, by its definition, is not mesh sensitive as damage mechanics. Camacho and Ortiz (1996) proposed CZM similar to yours yet not diffused. What are the benefits of diffused CZM? The fact that you need extra strength constraint is a disadvantage typical of CZM. The whole idea of continuum damage mechanics is the possibility to completely describe failure and fracture by using constitutive equations only. That makes CDM advantageous as compared to cohesive surfaces, which require additional conditions of their insertion. Your diffused CZM requires fine meshes while the classical CZM does not... 

Jian-Ying Wu's picture

I agree with Kosta -- if the phase-field length scale in AT1 is fixed as a material property related to lch in order to match the strength, then "convergence" is lost.

Konstantin Volokh's picture

By the way, you do not lose convergence in your theory because it is not realy phase-field in the sense of damage mechanics. Its a strange animal :))) If it is a sort of CSM then you seemengly need extra criteria for crack nucleation, branching, arrest... Your tension strength is one of such criteria. Do not be obsessed with it! :))) For example, a concrete cube in uniaxial compression fractures without having any tension stresses!

Jian-Ying Wu's picture

Sure. But I do not think any phase-field model with a a single damage variable is able to model various failure mechanisms simultaneously, like the compressive one you mentioned. If needed, it can be refined to include additional damage variables, but this is another story.

Konstantin Volokh's picture

Fracture is unsolved problem for a century. I doubt we will solve it soon :))) We are not alone in suffering. Fluid mechanicians are stuck with turbulence...

Jian-Ying Wu's picture

yes, I totally agree. This is what Prof. Bazant said in his speech for Timoshenko's medal when I was a visiting scholar in NU with him.

Jian-Ying Wu's picture

yes, you are right this time. Theoretically, things are more complicated and localization can occur anywhere. But numerically there always exist numerical errors and it is not necesary to introdue explicitly the imperfection as in other methods -- numerical errors trigs the localization interior somewhere if both ends are constrined to be free of damage (otherwise, localization occurs at either end).


Jian-Ying Wu's picture

Dear Kosta,

I asked Phu to upload a picture (I don't know how to do it) and you will understand why the PF-CZM is able to model crack nucleation.



Konstantin Volokh's picture

Yes, grad(d)=0 before the crack nucleation but d is not zero! It deviates from zero and this deviation leads to material instability (with the subsequent localization of damage into crack).

Emilio Martínez Pañeda's picture

Thank you Phu for this nice overview. It is probably worth emphasising that quasi-Newton methods such as the BFGS algorithm enable robust and efficient implementations of phase field fracture also for AT2 and AT1 models:

Jian-Ying Wu's picture

Thank you for the note, Emilio. You are right that BFGS monolithic algorithm is promising in solving the governing equations of PFMs.

Blaise Bourdin's picture

I strongly disagree with the statement that Figure 3 from ref 34 "demonstrates" that standard phase-field models cannot capture nucleation properly. As a matter of fact, [Tanné et al 2018], conveniently ignored in this post, makes a compelling case that this is not the case. If anythiong (and as noted in [34]), Figure 3 highlight the fact that the L-shaped experiment in concrete from Winkler thesis is a poor choice for the validation of a brittle fracture model, as representing concrete as a brittle material in a sample of 50cm x 50cm is inconsistent with the idea of a material length scale of the order of K_{Ic}/ \sigma_c, which ranges between 1 and 10 cm in typical concrete if memory serves me right.

Note that [Kumar et al., 2020] proposes an alternative approach to nucleation that does not require the introcudtion of a cohesize zone model, although at teh cost of teh variational nature of the problem.



[Kumar et al., 2020] Kumar, A., Bourdin, B., Francfort, G. A., and Lopez-Pamies, O. (2020). Revisiting nucleation in the phase-field approach to brittle fracture. J. Mech. Phys. Solids, 142:104027.

[Tanné et al 2018] Tanné, E., Li, T., Bourdin, B., Marigo, J.-J., and Maurini, C. (2018). Crack nucleation in variational phase-field models of brittle fracture. J. Mech. Phys. Solids, 110:80–99.

Jian-Ying Wu's picture

Dear Prof. Bourdin,

Thank you for your comments. The AT1 model is indeed able to deal with crack nucleation for most brittle fracture, provided the phase-field length scale is related to Irwin's characteristic length lch. In this case, the phase-field length scale is generally small enough not polluting the crack path. However, this is NOT always the case even for brittle fracture -- sometimes it is not possible to reconcile both needs to match the critical strength and to be small enough; see the example of indentation fracture [Strobl and Seelig, 2020]. 

About concrete, unless the structure size is exetrmely large, e.g., dams, concrete cannot be brittle and some quasi-brittle behavior is inevitable. For such cohesive fracture, I think it is not easy to fit both crack paths and peak loads for the standard AT1/AT2 models. If I am wrong, please let me know.

Anyway, we benifit a lot from variational phase-field model for fracture you, Francfort and Marigo pioneered over 20 years ago. And we are basically engineering oriented and might be not so rigorous sometimes. Thank you for your great contributions to this community.

[Strobl and Seelig, 2020] Strobl, M., Seelig, T., 2020. Phase field modeling of hertzian indentation fracture. Journal of the Mechanics and Physics of Solids 143, 104026.

Blaise Bourdin's picture

I don't disagree with your first statement, although I think that it is relatively easy to build generalized Ambrosio Tortorelly models such that the internal length to strength relation leads to a much smaller length. 


I think that a classical confusion is in the applicability of brittle models in nominally brittle situations. Very few materials are perfectly brittle. For nominally brittle materials, the brittle approximation only makes sense at a large enough scale (typically quantified by the size of the structure compared to that of the cohesize or internal length. So even though concrete at the scale of a dam can probably be seen as a brittle material, at the scale of Winkler's experiments, I have strong doubts. So using this experiment to "demonstrate" that phase-field models, seen as approximation of brittle fracture, cannot properly account for nucleation is a dubious claim in my opinion.

In order to study nucleation in L-shaped sample, I would actually recomment the problem from Figure 2(c) in Gómez, F. J., Elices, M., Berto, F., and Lazzarin, P. (2009). Fracture of V-notched specimens under mixed mode (I+II) loading in brittle materials. Int. J. Fracture, 159(2):121–135 is a much better choice that Winkler's experiment.

Jian-Ying Wu's picture

I agree with you on the second point. Application of the AT1/2 models to Concrete or more generally, quasi-brittle solids (even the so-called brittle solids exhibit quasi-brittle behavior at smaller scale) might be inappropriate.

For the first point, just as the cohesive zone model applies to both brittle fracture with a "small" internal length scale and quasi-brittle one with a "larger one", it might be more helpful if such a counterpart exists in the phase-field framework.

Dear Prof. Nguyen and Prof. Wu


Thank you for your sharing. I am reading your work recently and have several questions about the governing equations. According to your paper [1], the governing equations for phase-field are as follows,

However, these governing equations only apply in the localization band domain. The governing eq for the domain outside the band is not given. Besides, in the initial elastic stage, there is no localization band. So what is the governing equation then? If I understand the problem correctly, the governing equations for the domain outside the band and for the elastic stage are [2]

Due to α'(0)>0, we have d=0 outside the localization crack band all the time. This is also the reason why α(d)=2d-d^2 leads to an elastic stage while α(d)=d^2 doesn't.  If I made mistakes, please point them out. 

Here comes the first question, you said in your paper [3] that "However, this strategy is not necessarily mandatory. When the damage sub-domain Bcannot be easily selected, the whole computational domain can be used." Do you mean that the governing eq (1-2) can be applied to the whole computational domain? If so, can you explain the reason? From my point of view, eq(1-2) is not equivalent to eq(3-4) outside the crack band. 

The second question is for eq (1-2), there is a trivial solution d=1 in B. So how to avoid this solution in computation? 

I am new to this area and if you could answer my questions I will appreciate it very much.


Best Regards,




[1] Jian-Ying Wu, Yuli Huang, Hao Zhou, Vinh Phu Nguyen, Three-dimensional phase-field modeling of mode I + II/III failure in solids, Computer Methods in Applied Mechanics and Engineering, Volume 373, 2021, 113537, ISSN 0045-7825,

[2] Jian-Ying Wu, A unified phase-field theory for the mechanics of damage and quasi-brittle failure, Journal of the Mechanics and Physics of Solids, Volume 103, 2017, Pages 72-99, ISSN 0022-5096,

[3] Jian-Ying Wu, Jie-Feng Qiu, Vinh Phu Nguyen, Tushar Kanti Mandal, Luo-Jia Zhuang, Computational modeling of localized failure in solids: XFEM vs PF-CZM, Computer Methods in Applied Mechanics and Engineering, Volume 345, 2019, Pages 618-643, ISSN 0045-7825,

phunguyen's picture

Dear Tianyu,

Equations (1,2) hold for the entire domain, just that outside the damged region B, d=0 identically. The reason that Jian-Ying always write (1,2) because he wanted to emphasize that the damage equation is confined only to B, a much smaller sub-domain. Plus, in our implementation, we exploit this fact to have two sub-domains: B where elements with damage dofs and the remaining with elements without damage dofs.

We can have d=1, I do not see any problem with that.



Dear Prof. Nguyen and Prof. Wu


Thanks for your reply! I think you have made it clear. I can see now why your algorithm set a lower bound for the history variable H. So you are trying to use that lower bound for H to set d=0 outside the crack band so that you can create the band automatically. It looks very cool! But I still have one question, is this technique used just for a numerical purpose? Does it have any physical meaning?  


Best Regards,


phunguyen's picture

I think that, the H history variable is an easy way to impose that below certain level of deformation there is no fracture/damage. Is it ok for you?

Yes. Thank you again for your sharing! :)



L. Roy Xu's picture

Dear Profs. Nguyen and Wu,

Interesting simulation approach. As an experimentalist on fracture and impact mechanics, I always try to learn some new simulation tools,  especially XFE, CZ and PD cannot simulate my previous dynamic fracture experiments well.

For example, I have some unique crack nucleation (from zero length to a finite length), crack branching (branching angle is larger than the max. theoretical angle), sharply curved cracks which exceed the simulation capability of XFEM.  

I found you’re simulating old crack branching experiments in 1991. I suggest that we may have collaboration in the future. I attached three papers on dynamic fracture experiments co-authored with Prof. Rosakis at the California Institute of Technology.

1 L. R. Xu and A. J. Rosakis,  “An Experimental Study on Dynamic Failure Events in Homogeneous Layered Materials Using Dynamic Photoelasticity and High-speed Photography,” Optics and Laser in Engineering, Vol. 40, pp. 263-288,2003. ( PDF file)

2. L. R. Xu, Y. Y. Huang and A. J. Rosakis, “Dynamic Crack Deflection and Penetration at Interfaces in Homogeneous Materials: Experimental Studies and Model Predictions,” Journal of the Mechanics and Physics of Solids, Vol. 51, pp.425-460, 2003. ( PDF file)

3. L. R. Xu and A. J. Rosakis, “Real-time Experimental Investigation of Dynamic Crack Branching Using High-speed Optical Diagnostics,” SEM Experimental Techniques, Vol. 27,pp.23-26, 2003. (PDF file)

My email is Look forward to meeting you in the future.    

Roy Xu

phunguyen's picture

Dear Prof. Xu,

I have sent you an email regarding your suggestion. 

Best regards,


Subscribe to Comments for "Journal Club For April 2021:  Variational phase-field modeling of brittle and cohesive fracture"

Recent comments

More comments


Subscribe to Syndicate