Implementation aspects of a phase-field approach for brittle fracture

Abstract

This paper provides a comprehensive overview of a phase-field model of fracture in solid mechanics setting. We start reviewing the potential energy governing the whole process of fracture including crack initiation, branching or merging. Then, a discretization of system of equation is derived, in which the key aspect is that for the correctness of fracture phenomena, a split into tensile and compressive terms of the strain energy is performed, which allows crack to occur in tension, not in compression. For numerical analysis, standard finite element shape functions are used for both primary fields including displacements and phase field. A staggered scheme which solves the two fields of the problem separately is utilized for solution step and illustrated with a segment of Python code.

Keywords

phase-field modeling / FEM / staggered scheme / fracture

Introduction

The prediction of failure mechanism due to fracture is of great importance in engineering designs. From the physical point of view, the fracture failure is a multiscale phenomenon [¹,²], where the fracture starts at the microscopic level, i.e., crack nucleation, and propagates to the macroscopic level in the form of observable fracture. The fracturing process is characterized by the reduction of load-bearing capacity of the structure, leading to partially or totally collapse at the macroscopic level. Because of the complex character of the fracturing process, a wide variety of fracture models have been suggested. Griffith [³] was considered as a pioneer for developing the theory of brittle fracture. The general principle is that the energy stored in a solid body is the sum of the elastic strain energy and the dissipated energy needed to create newly cracked surfaces. Furthermore, it is assumed that the crack propagates when the strain energy reaches a critical value. In the effort to simplify the structural calculation due to the consideration of the whole system energy with a known crack, Irwin [⁴] introduced the concept of stress intensity factor K, which considers the critical value of the stress in the vicinity of a crack tip as a function with respect to the crack length and the loading. Though these theories give criteria for crack growth, they are not capable of predicting the crack initiation.

Numerical method is a useful tool to simulate the fracture phenomena such as crack initiation, crack branching, dynamic cracking. However, there are challenges in capturing crack paths accurately, especially the curved crack paths. There are two common types of numerical model for fracture, namely discrete and diffusive approaches. The discrete approach is used to describe the sharp crack topology by incorporating a discontinuity in the displacement field. Typically successful methods of this approach includes the extended finite element method (XFEM) proposed by Moës et al. [⁵], meshfree methods developed by Rabczuk et al. [⁶–¹⁰], the cohesive zone method proposed by Remmers et al. [¹¹], the higher order gradient damage model developed by Thai et al. [¹²], extended isogeometric analysis (XIGA) by Ghorashi et al. [¹³], the interface finite element method suggested by Xu and Needleman [¹⁴], and the work of Miehe and Gürses [¹⁵] and Areias and Rabczuk [¹⁶] with a configurational-force driven model. Although they are able to successfully simulate the crack growth in two-dimensional cases, it has been shown that they suffer some drawbacks in three dimensions with crack branching where they need to capture the extent of the crack during computing process. When it comes to diffusive approach, e.g., phase field model, this issue is alleviated compared with the aforementioned approach in the sense that the discontinuities are not needed to be tracked during simulation. The fracture states, such as crack initiation and crack kinking, can be simulated in a straight-forward manner because its description is embedded in the variational form. For the diffusive approach, the sharp crack topology is approximated in the distributed sense by a so-called phase-field variable transforming from the unbroken to fully broken state of the material. This parameter is driven by a driving force, which is a function of strain energy. Moreover, with the appearance of this new quantity, the physical problem at hand has to be redefined in the form of a couple equation system with a momentum equation and a gradient-type evolution equation. The length scale parameter is introduced in the phase-field model to control the width of the regularized crack and to match the crack dissipation energy from experiment. It is assumed that when the length scale parameter closes to zero, the phase-field approach converges to Griffith-type fracture model. Another notable feature of this approach is that the evolution of the discontinuities is captured automatically by the evolution of a smooth scalar-valued phase-field without adaptive mesh strategies. Therefore, this capability avoids the difficulty in tracking complicated curved crack paths in three-dimensional cases. Clearly, this approach shows noticeable superiority over the discrete fracture models, which requires the implicit or explicit description of the discontinuities. For brittle fracture, phase-field formulation originated from the work of Francfort and Marigo [¹⁷], see also Ref. [¹⁸], which are based on Griffith-type energy minimization using variational methods. Besides that, the other variants of phase-field approach include the classical Ginzburg-Landau type evolution equation as introduced in the work of Hakim and Karma [¹⁹]. Until recently, this approach has been extended to various applications for both static [²⁰,²¹] and dynamic [²²] fracture problems in brittle material, for finite strain behaviour of fracture [²³–²⁹], structural mechanics as thin plates and shells [¹⁶,³⁰,³¹], nano-composite material [³²–³⁴], ductile material [³⁵,³⁶], hydraulic fracturing [³⁷], and contact problems [³⁸].

The goal of this paper is for the detailed implementation of a phase-field fracture model in solid mechanics setting, which is considered as a basic foundation that is necessary for various applications mentioned above and it is aimed at beginners for such diffusive fracture problems. We implement the phase-field problem on Kratos, an open sourced framework based on finite element method, which provides all necessary components for problems using finite element method such as solver, post- and pre-processing and assembly of system matrix. Therefore, we just transform the formulation at hand to a so-called custom-elements modul in Kratos. In particular, to demonstrate the concept of staggered scheme used to solve such a multi-field problem, a segment of Python code is shown in Appendix A of this paper.

A review of the phase-field model for brittle fracture

Following the work of Ref. [²⁰], this section provides detailed steps leading to discrete equations for the employed phase-field model. It starts with a review of the variational formulation based on the incremental variational approach. Quasi-static problem is assumed, in which the loading is subjected to act slowly, so the inertia term is neglected in the energy system. Accordingly, the weak form is derived by taking the derivative of energy functional with respect to time.

Problem statement

Let W⊂ R^nd be an arbitrary body in an nd-dimensional medium. G represents a discrete crack incorporated into the domain of the problem and ∂W is the external boundary of the body. Surface tractions t are imposed on the Neumann boundary ∂W^N and prescribed displacements u are described on the Dirichlet boundary ∂W^D. Assuming the small strain context, a coupled two-field problem describing a phase field fracture model is characterized by the displacement field,

(1)\mathit{u}=\{\begin{array}{l}\Omega \times T\to {\mathit{R}}^{nd}\\ (x,t)\mapsto \mathit{u}(x,t)\end{array},

and the phase-field which takes on the value of 1 corresponding to the intact state of materials and of 0 corresponding to the fully broken state of materials as depicted in Fig. 1.

(2)c:=\{\begin{array}{l}\Omega \times T\to [\mathrm{0},\mathrm{1}]\\ (x,t)\mapsto c(x,t)\end{array}\mathrm{.}

The symmetric strain tensor is expressed in terms of the gradient of the displacement as

(3)\epsilon ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\nabla \mathit{u}+\nabla {\mathit{u}}^T)\mathrm{.}

Energy functional

To model the process of crack initiation, propagation and branching, the minimization problem of an energy functional is established through the equation as

(4)\mathrm{\Pi }=ϵ+\mathcal{F}-{\mathcal{P}}_{\mathrm{e}\mathrm{x}\mathrm{t}},

where P is the total potential energy, \epsilon is the strain energy, ℱ is the fracture energy, and {\mathcal{P}}_{\mathrm{e}\mathrm{x}\mathrm{t}} is the external work.

In the following steps, each term in Eq. (4) will be explained thoroughly. Firstly, components of the fracture energy should be analysed, by which a scalar-valued phase field is introduced. Following the work of Ref. [³], a unit area is created by the work done due to fracture. This amount of energy dissipation is denoted as critical energy release rate G_c. The fracture energy can be postulated by

(5)ℱ={\displaystyle \underset{\Omega }{\int }{G}_c}d\Gamma \mathrm{.}

From the numerical point of view, it is non-trivial or not straightforward to find the minimization of Eq. (4) with constraint Eq. (5). This is due to the unpredictable character of the crack topology G. To further simplify the formulation, the surface integral is substituted by a volumetric integral, which is emanated from the idea of Ref. [¹⁸]:

(6){\displaystyle \underset{\Gamma }{\int }{G}_c}d\Gamma \approx {\displaystyle \underset{\Omega }{\int }{G}_c{\Gamma }_c}\mathrm{d\Omega },

with G_c considered a crack density function and introduced by Bourdin et al. [¹⁸] as

(7){\Gamma }_c(c,\nabla c)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{{\displaystyle l}}_{\mathrm{0}}}}{(c-\mathrm{1})}^{\mathrm{2}}+l_{\mathrm{0}}^{\mathrm{2}}{\left|\nabla c\right|}^{\mathrm{2}}\mathrm{.}

Equation (7) uses l₀ as a length scale parameter controlling the width of the transitional zone from the uncracked to cracked states of materials and depends on the phase-field variable and its gradient.

Concerning the strain energy\epsilon, it is defined as follows

(8)\epsilon ={\displaystyle \underset{\Omega }{\int }{\psi }_e}(\epsilon )d\Omega ,

with the strain energy density y_e for isotropic linear elasticity being expressed as

(9){\psi }_e(\epsilon )={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\lambda {(tr\epsilon )}^{\mathrm{2}}+\mu tr{\epsilon }^{\mathrm{2}}\mathrm{.}

A key aspect to note in fracture mechanics is that cracks only occur in tension, not in compression due to the fact that tensile stresses cannot be able to bring through a crack opening, but the compressive counterpart can. Hence, the strain energy should be split into tensile and compressive terms, which leads to the need of a split in the strain tensor as

(10)\epsilon ={\epsilon }^{+}+{\epsilon }^{-},

with {\epsilon }^{+} and {\epsilon }^{-} being computed through a spectral decomposition introduced by Miehe et al. [¹⁵] and De Souza et al. [³⁹] as

(11){\epsilon }^{\pm }(\epsilon )={\sum \langle {\epsilon }_i\rangle }^{\pm }{n}_i\otimes n_i,

where e_i is eigenvalues and n_i is eigenvector of the strain tensor. The ramp functions {\langle ·\rangle }^{+} and {\langle ·\rangle }^{-} are defined as

(12)

(13){\langle x\rangle }^{-}=\{\begin{array}{l}\mathrm{0}\hspace{0.17em}\begin{array}{ll}& \end{array}\mathrm{if}\hspace{0.17em}x>\mathrm{0}\\ x\hspace{0.17em}\begin{array}{ll}& \end{array}\mathrm{if}\hspace{0.17em}x\le \mathrm{0}\end{array}\mathrm{.}

With the split operator of strain tensor, the positive and negative terms of strain energy density and stress tensor can be expressed as follows

(14){\psi }_e{(\epsilon )}^{\pm }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\lambda {({\langle tr\epsilon \rangle }^{\pm })}^{\mathrm{2}}+{\mu tr\epsilon }^{{\pm }^{\mathrm{2}}},

(15){\sigma }^{\pm }=\lambda {\langle tr\epsilon \rangle }^{\pm }+{\mathrm{2}\mu \epsilon }^{\pm }\mathrm{.}

The tensile components are reduced by a degradation function g(c, k) = (1-k)c²+ k as follows

(16){\psi }_e=g(c,k){\psi }_e^{+}+{\psi }_e^{-},

(17)\sigma =g(c,k){\sigma }^{+}+{\sigma }^{+}\mathrm{.}

Finally, the external work P_ext in Eq. (4) takes the form of

(18){{\displaystyle P}}_{\mathrm{ext}}={\displaystyle \underset{\Omega }{\int }\overline{b}}\cdot ud\Omega +{\displaystyle \underset{\partial \Omega }{\int }\overline{t}}\cdot u,

where \bar{b} is a body force.

Variational formulation

With all components in the energy functional expressed above, we obtain the rate of the potential energy, from which the weak form of the problem is obtained, by taking the derivative of the Eq. (4) with respect to time t as

(19){\displaystyle \frac{\partial \Pi }{\partial t}}={\displaystyle \frac{\partial \epsilon }{\partial t}}+{\displaystyle \frac{\partial F}{\partial t}}+{\displaystyle \frac{\partial {{\displaystyle P}}_{\mathrm{ext}}}{\partial t}}\mathrm{.}

The rate of the first term on the right-hand side of Eq. (19) can be obtained by

(20)\begin{array}{l}\dot{\epsilon }={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\displaystyle \underset{\Omega }{\int }{\psi }_e}(\epsilon ,c)\\ ={\displaystyle \underset{\Omega }{\int }{\displaystyle \frac{\partial {\psi }_e}{\partial \epsilon }}}{\displaystyle \frac{\partial \epsilon }{\partial t}}\mathrm{d}\Omega +{\displaystyle \underset{\Omega }{\int }{\displaystyle \frac{\partial {\psi }_e}{\partial c}}}{\displaystyle \frac{\partial c}{\partial t}}\mathrm{d}\Omega \\ ={\displaystyle \underset{\Omega }{\int }{{\displaystyle \sigma \nabla }}_s}\dot{\mathit{u}}\mathrm{d}\Omega +{\displaystyle \underset{\Omega }{\int }\mathrm{2}c(\mathrm{1}-k){\psi }_e^{+}}\dot{c}\mathrm{d}\Omega \mathrm{.}\end{array}

Similarly, we obtain the derived form of the rate of the dissipated energy due to fracture by

(21)\begin{array}{l}\dot{ℱ}={\displaystyle \frac{d}{dt}}{\displaystyle \underset{\Omega }{\int }{\psi }_c}\mathrm{d\Omega }\\ ={\displaystyle \frac{d}{dt}}{\displaystyle \underset{\Omega }{\int }{G}_c{\Gamma }_c}(c,\nabla c)\mathrm{d\Omega }\\ ={\displaystyle \underset{\Omega }{\int }{G}_c}{\displaystyle \frac{\partial {\Gamma }_c}{\partial c}}\dot{c}\mathrm{d\Omega }+{\displaystyle \underset{\Omega }{\int }{G}_c}{\displaystyle \frac{\partial {\Gamma }_c}{\partial \nabla c}}\nabla \dot{c}\mathrm{d\Omega }\\ ={\displaystyle \underset{\Omega }{\int }{\displaystyle \frac{G_c}{\mathrm{2}{l}_{\mathrm{0}}}}}(c-\mathrm{1})\dot{c}\mathrm{d\Omega }+{\displaystyle \underset{\Omega }{\int }\mathrm{2}}{G}_c{l}_{\mathrm{0}}(\nabla c)\left(\nabla \dot{c}\right)\mathrm{d\Omega },\end{array}

and the external power due to the external load is simply expressed as

(22){{\displaystyle \dot{\mathit{P}}}}_{\mathrm{ext}}={\displaystyle \underset{\mathrm{\Omega }}{\int }\overline{\mathit{b}}\cdot \dot{\mathit{u}}\mathrm{d}\mathrm{\Omega }+{\displaystyle \underset{\partial \mathrm{\Omega }}{\int }\overline{\mathit{t}}}}\cdot \dot{\mathit{u}}\mathrm{.}

All terms of the rate of the potential energy Eq. (19) are rearranged in a form such that they only have a common term of the rate of the displacements \dot{u} or the rate of the phase field \dot{c}. This leads to the separation of the rate of the potential energy into an elastic contribution and a phase-field contribution. The elastic contribution takes the form of

(23){\displaystyle \underset{\mathrm{\Omega }}{\int }\mathit{\sigma }:{{\displaystyle \nabla }}_s}\dot{\mathit{u}}\mathrm{d\Omega }={\displaystyle \underset{\mathrm{\Omega }}{\int }\mathit{b}\cdot \dot{\mathit{u}}\mathrm{d\Omega }+{\displaystyle \underset{\partial {\mathrm{\Omega }}_t}{\int }\mathit{t}\cdot \dot{\mathit{u}}\mathrm{d}\partial \mathrm{\Omega }}},

whereσ is the stress tensor, b the body force. The phase field contribution is represented by

(24){\displaystyle \underset{\mathrm{\Omega }}{\int }\left({\displaystyle \frac{\mathrm{4}{l}_{\mathrm{0}}(\mathrm{1}-k){\psi }_e^{+}}{G_c}}+\mathrm{1}\right)}c\dot{c}\mathrm{d\Omega }+{\displaystyle \underset{\mathrm{\Omega }}{\int }\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}}(\nabla c)\left(\nabla \dot{c}\right)\mathrm{d\Omega }={\displaystyle \underset{\mathrm{\Omega }}{\int }\dot{c}\mathrm{d\Omega }}

Employing integration by parts to Eqs. (23) and (24), the strong form of the coupled two field problem and boundary conditions for the phase field fracture model are obtained accordingly. The first term of Eq. (23) can be written in a form of

(25){\displaystyle \underset{\mathrm{\Omega }}{\int }\mathit{\sigma }:{{\displaystyle \nabla }}_s}\dot{\mathit{u}}\mathrm{d\Omega }={\displaystyle \underset{\partial \mathrm{\Omega }}{\int }(\mathit{\sigma }\mathit{n})}\cdot \dot{\mathit{u}}\mathrm{d}\partial \mathrm{\Omega }-{\displaystyle \underset{\mathrm{\Omega }}{\int }\mathrm{div}\mathit{\sigma }\cdot }\dot{\mathit{u}}\mathrm{d\Omega },

and the second term of Eq. (24) can be obtained as

(26){\displaystyle \underset{\Omega }{\int }\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}}(\nabla c)\left(\nabla \dot{c}\right)\mathrm{d\Omega }={\displaystyle \underset{\Omega }{\int }\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}}(\nabla c\cdot n)\dot{c}\mathrm{d\Omega }-{\displaystyle \underset{\Omega }{\int }\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}}\Delta c\dot{c}\mathrm{d\Omega }\mathrm{.}

Substituting Eqs. (25) and (26) into Eqs. (23) and (24), respectively, we obtain

(27)-{\displaystyle \underset{\mathrm{\Omega }}{\int }(\mathrm{div}\mathit{\sigma }+\mathit{b})}\cdot \dot{u}\mathrm{d\Omega }+{\displaystyle \underset{\partial \mathrm{\Omega }}{\int }\left(\mathit{\sigma }\cdot \mathit{n}-\overline{\mathit{t}}\right)}\cdot \dot{\mathit{u}}\mathrm{d}\partial \mathrm{\Omega }=\mathrm{0},

and

(28){\displaystyle \underset{\mathrm{\Omega }}{\int }\left({\displaystyle \frac{\mathrm{4}{l}_{\mathrm{0}}(\mathrm{1}-\kappa ){\psi }_e^{+}}{G_c}}+\mathrm{1}\right)c\dot{c}\mathrm{d\Omega }}+{\displaystyle \underset{\mathrm{\Omega }}{\int }\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}(\nabla c\cdot \mathit{n})\dot{c}\mathrm{d\Omega }}-{\displaystyle \underset{\mathrm{\Omega }}{\int }\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}}\Delta c\dot{c}\mathrm{d\Omega }={\displaystyle \underset{\mathrm{\Omega }}{\int }\dot{c}\mathrm{d\Omega }}\mathrm{.}

Based on Eqs. (27) and (28), the strong form and boundary conditions of the phase field fracture model are given by

(29)\{\begin{array}{l}\mathrm{div}\mathit{\sigma }+\mathit{b}=\mathrm{0}\begin{array}{llll}& & & \begin{array}{llll}& & & \end{array}\end{array}\mathrm{in}\begin{array}{l}\end{array}\mathrm{\Omega }\\ \left({\displaystyle \frac{\mathrm{4}{l}_{\mathrm{0}}(\mathrm{1}-\kappa ){\psi }_{\mathrm{e}}^{+}}{G_{\mathrm{c}}}}+\mathrm{1}\right)c-\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}\Delta c=\mathrm{1}\begin{array}{l}\end{array}\mathrm{in}\begin{array}{l}\end{array}\mathrm{\Omega }\\ \mathit{\sigma }\cdot \mathit{n}=\overline{\mathit{t}}\begin{array}{llll}& & & \begin{array}{ll}& \end{array}\end{array}\mathrm{on}\begin{array}{l}\end{array}\partial {\mathrm{\Omega }}_{\mathrm{N}}\\ \mathit{u}=\overline{\mathit{u}}\hspace{0.17em}\begin{array}{llll}& & & \begin{array}{ll}& \end{array}\end{array}\mathrm{on}\begin{array}{l}\end{array}\partial {\mathrm{\Omega }}_{\mathrm{D}}\\ \nabla d\cdot \mathit{n}=\mathrm{0}\begin{array}{llll}& & & \end{array}\mathrm{on}\begin{array}{l}\end{array}\partial \mathrm{\Omega }\end{array}\mathrm{.}

In order to guarantee the irreversibility of the crack evolution, a history field is introduced, representing the maximum of the positive strain energy density during loading applied within the time instant t

(30)\mathscr{H}(x,t)=\underset{s\in [\mathrm{0},t]}{\max }{\psi }_{\mathrm{0}}^{+}(\epsilon (\mathit{x},s))\mathrm{.}

Equation (30) facilitates the irreversibility condition

(31)\dot{ℱ}\ge \mathrm{0},

leading the degradation of the fracture energy from the energy system.

Discretization

To find a solution for the strong form composing of the momentum equation and the phase field equation, function spaces are defined as

(32)\{\begin{array}{l}{\mathcal{V}}^u=\left\{\mathit{u}\in {\mathit{H}}^{\mathrm{1}}(\mathrm{\Omega })|\mathit{u}=\overline{\mathit{u}}\begin{array}{l}\end{array}\mathrm{on}\begin{array}{l}\end{array}\partial {\mathrm{\Omega }}^{\mathrm{D}}\right\}\\ {\mathcal{V}}^c=H^{\mathrm{1}}(\Omega )\end{array},

for the approximated solutions of the problem, and

(33)\{\begin{array}{l}{{\displaystyle \mathcal{V}}}_{\mathrm{0}}^u=\left\{\dot{\mathit{u}}\in {\mathit{H}}^{\mathrm{1}}(\mathrm{\Omega })|\dot{\mathit{u}}=\mathbf{0}\begin{array}{l}\end{array}\mathrm{on}\begin{array}{l}\end{array}\partial {\mathrm{\Omega }}^{\mathrm{D}}\right\}\\ {{\displaystyle \mathcal{V}}}_{\mathrm{0}}^c=\left\{\dot{c}\in {\mathit{H}}^{\mathrm{1}}(\mathrm{\Omega })|\dot{c}=\mathbf{0}\begin{array}{l}\end{array}\mathrm{on}\begin{array}{l}\end{array}\partial \mathrm{\Omega }\right\}\end{array},

for the trial functions of the problem. The displacements and phase field are discretized as follows

(34)u(x)=N^u(x)d^u,c(x)=N^c(x)d^c,

(35)\dot{u}(x)=N^u(x)w^u,c(x)=N^c(x)w^c\mathrm{.}

Phase-field part

Replacing the phase-field approximation and its corresponding test function into the weak form of the phase field part, we obtain

(36){\displaystyle \underset{\mathrm{\Omega }}{\int }\left[{\displaystyle \frac{\mathrm{4}{l}_{\mathrm{0}}(\mathrm{1}-\kappa )\mathscr{H}}{g_c}}+\mathrm{1}\right]}{({\mathit{N}}^c{\mathit{w}}^c)}^{T}c\mathrm{d\Omega }+{\underset{\mathrm{\Omega }}{\int }[\nabla ({\mathit{N}}^c{\mathit{w}}^c)]}^T\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}(\nabla c)\mathrm{d\Omega }={\displaystyle \underset{\mathrm{\Omega }}{\int }{\mathit{N}}^c{\mathit{w}}^c\mathrm{d\Omega }}\mathrm{.}

Then the residual form for phase-field part is expressed as follows

(37){\mathit{r}}^c={{\displaystyle \mathit{f}}}_{\mathrm{int}}^c-{{\displaystyle \mathit{f}}}_{\mathrm{ext}}^c,

where the internal and external forces take the forms

(38){{\displaystyle \mathit{f}}}_{\mathrm{int}}^c={\displaystyle \underset{\Omega }{\int }\left[{\displaystyle \frac{\mathrm{4}{l}_{\mathrm{0}}(\mathrm{1}-\kappa )\mathscr{H}}{g_c}}+\mathrm{1}\right]{({\mathit{N}}^c)}^{T}c\mathrm{d\Omega }}+{\displaystyle \underset{\Omega }{\int }\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}}(\nabla c)\cdot ({\nabla N}^c)\mathrm{d\Omega },

(39){{\displaystyle \mathit{f}}}_{\mathrm{ext}}^c={\displaystyle \underset{\Omega }{\int }{\mathit{N}}^c}\mathrm{d\Omega }\mathrm{.}

And the stiffness matrix is obtained by taking the derivative of r^c with respect to d^c

(40)\begin{array}{l}{\mathit{K}}^c={\displaystyle \frac{{\partial \mathit{r}}^c}{{\partial \mathit{d}}^c}}\\ ={\displaystyle \underset{\mathrm{\Omega }}{\int }\left[{\displaystyle \frac{\mathrm{4}{l}_{\mathrm{0}}(\mathrm{1}-\kappa )\mathscr{H}}{g_c}}+\mathrm{1}\right]{({\mathit{N}}^c)}^T}({\mathit{N}}^c{\mathit{d}}^c)\mathrm{d\Omega }+{\displaystyle \underset{\Omega }{\int }\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}{(\nabla {\mathit{N}}^c)}^T}(\nabla {\mathit{N}}^c{d}^c)\mathrm{d\Omega }\\ ={\displaystyle \underset{\mathrm{\Omega }}{\int }\left[{\displaystyle \frac{\mathrm{4}{l}_{\mathrm{0}}(\mathrm{1}-\kappa )\mathscr{H}}{g_c}}+\mathrm{1}\right]{({\mathit{N}}^c)}^T{\mathit{N}}^c}\mathrm{d\Omega }+{\displaystyle \underset{\Omega }{\int }\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}{(\nabla {\mathit{N}}^c)}^T}(\nabla {\mathit{N}}^c)\mathrm{d\Omega }.\end{array}

Elastic part

Analogously, substituting the displacement field and its corresponding test function into the weak form of the elastic part, we obtain

(41){\displaystyle \underset{\mathrm{\Omega }}{\int }\mathit{\sigma }:{{\displaystyle \nabla }}_s}({\mathit{N}}^u{\mathit{w}}^u)\mathrm{d\Omega }={\displaystyle \underset{\mathrm{\Omega }}{\int }\mathit{b}\cdot ({\mathit{N}}^u{\mathit{w}}^u)}\mathrm{d\Omega }+{\displaystyle \underset{\partial {\mathrm{\Omega }}_t}{\int }\mathit{t}\cdot }({\mathit{N}}^u{\mathit{w}}^u)d\partial \mathrm{\Omega }\mathrm{.}

The residual form for elastic part is represented by

(42){\mathit{r}}^u={{\displaystyle \mathit{f}}}_{\mathrm{int}}^u-{{\displaystyle \mathit{f}}}_{\mathrm{ext}}^u,

with the internal and external forces as follows:

(43){{\displaystyle \mathit{f}}}_{\mathrm{int}}^u={\displaystyle \underset{\mathrm{\Omega }}{\int }\mathit{\sigma }:({{\displaystyle \nabla }}_s{\mathit{N}}^u)}\mathrm{d\Omega },

(44){{\displaystyle \mathit{f}}}_{\mathrm{ext}}^u={\displaystyle \underset{\mathrm{\Omega }}{\int }\mathit{b}\cdot {\mathit{N}}^u\mathrm{d\Omega }+{\displaystyle \underset{\partial {\mathrm{\Omega }}_t}{\int }\mathit{t}\cdot }{\mathit{N}}^u}\mathrm{d}\partial \mathrm{\Omega }\mathrm{.}

The corresponding stiffness matrix is given by taking the derivative of Eq. (42) with respect to d^u

(45)\begin{array}{l}{\mathit{K}}^u={\displaystyle \frac{\partial r^u}{{\partial d}^u}}={\displaystyle \underset{\Omega }{\int }{({\nabla }_s{N}^u)}^T{\displaystyle \frac{\partial \sigma }{\partial \epsilon }}}{\displaystyle \frac{\partial \epsilon }{\partial d^u}}\mathrm{d\Omega }\\ ={\displaystyle \underset{\Omega }{\int }{({\nabla }_s{N}^u)}^T}C({\nabla }_s{N}^u)\mathrm{d\Omega },\end{array}

where s is the Cauchy stress and C is the material tangent.

Numerical examples

The objective of this section is to confirm the efficiency and robustness of the employed phase field fracture model using the finite element method by selected numerical benchmarks. Firstly, two classical benchmark problems, which are single edge notched tension and shear tests, are investigated. Subsequently, other three examples consisting of three-point bending test, asymmetric three-point bending test and notched plate with hole are presented.

Single edge notched tension test

The specimen of this test is described by a squared plate with a horizontal notch located at middle height of the left edge stretching to the center of the specimen. The geometry setup and applied deformation are depicted as in Fig. 2 with vertical displacements imposed on all points of the top edge. A mesh, which is refined in zones where crack is expected to grow, contains 25000 elements with an effective element size of h≈ 0.005 mm in the critical area. Material parameters are chosen as in Table 1. The computation is conducted under displacement control with a displacement increment of Du =1 × 10^-5 mm in the first 500 loading steps and of Du =1 × 10^-6 mm in the subsequent loading steps. The contour plots for the evolution of the crack path at several stages of imposed displacements are given in Fig. 3. Figure 4 shows the load deflection curve for this example.

Single edge notched pure shear test

We now consider the same specimen as in previous example, but horizontal displacement is imposed on the complete top edge instead of vertical displacement. Figure 5 illustrates the geometric set up and boundary conditions of the problem. A discretization consists of 25000 triangular elements and is refined in the expected crack propagation. An effective mesh of h≈ 0.004 mm which is one half of the length scale is chosen to track properly crack propagation. The simulation is performed under displacement control with displacement increment of Du =1 × 10^-5 mm during the process of applying loading. The material parameters are chosen as in Table 1. Figure 6 shows the load-deflection curve which is in good agreement with the curve obtained in Ref. [¹⁵]. Crack patterns of the specimen at different stages of displacement are depicted as in Fig. 7.

Symmetric three-point bending test

This test is a classical benchmark problem investigated in phase field literature. The geometric set up and boundary conditions of the problem are shown as in Fig. 8. A mesh containing 30000 triangular elements is refined a priori in areas where crack is forecasted to propagate. An effective element size of h≈ 0.005 mm in the interest zone is selected to be one half of the length scale parameter. The computation is performed under displacement control with displacement increment of Du =1 ×10^-4 mm during the entirely computational process. Table 2 gives material parameters chosen for this problem. The resulting crack ranges at several stages of displacement are shown as in Fig. 9. Figure 10 illustrates the global response of the specimen.

Asymmetric notched three-point bending test

In this example, an asymmetric notched beam with three holes and being supported at two points is considered as in Fig. 11. Material parameters are chosen as in Table 3. The discretization consists of 60000 quadrilateral elements and is refined a priori in expected crack propagation zones. The simulation is performed under fixed displacement control with fixed displacement increments of Du =1 × 10^-3 mm for the first 200 loading steps and Du =1 × 10^-4 mm in the subsequent steps. Fig. 12 shows the load-deflection curve of this problem, while the progression of the crack path is depicted as in Fig. 13.

Notched plate test with hole

Next, we test the employed phase field model with the example of the notched plate with a hole. The geometric set up and boundary conditions are given as in Fig. 14, in which external loading is applied at a top pin and a lower pin is fixed with a larger hole between the two holes to yield a mixed-mode fracture case. Material parameters are selected as in Table 4. The mesh comprises of 80000 triangular elements and is refined a priori in forecasted crack growth areas. The computation is conducted in a monotonic displacement-controlled context with fixed displacement increments Du =1 × 10^-3 mm during the simulation process. Crack patterns at different stages of deformation are shown as in Fig. 15. One observes that a curvilinear crack propagates from the notch from the left edge to the largest hole of the specimen, and then crack faces tend to penetrate the closure towards the right edge. The resulting crack path and the load-deflection curve of the problem as depicted in Fig. 16 are compared to the result obtained from an experiment by Ambati et al. [³⁵].

Conclusions

This paper provided a comprehensive overview for an implementation of the employed phase-field approach in brittle fracture. Starting with reviewing the potential energy governing the whole process of fracture including crack initiation, branching or merging, a discretization of system of equation is presented, in which the key aspect is that for the correctness of fracture phenomena, a split into tensile and compressive terms of the strain energy is performed, which allows crack to occur in tension, not in compression. For the purpose of confirming the accuracy and efficiency of the approach, we tested various numerical examples, which shows good agreement with those obtained from previous work. Particularly, this paper is aimed at beginners on this approach, who want to attain the results quickly and apply to other applications in the computational mechanics fields. Concerning disadvantages of the method, the length scale parameter and the step size of each displacement increment are needed to choose suitably, and a local refinement for the mesh in critical regions where cracks are expected to propagate is required. Moreover, due to primary variables are updated separately, the staggered scheme accumulates a small sum of errors, but this can be overcome by the methods in Refs. [⁴⁰,⁴¹].

References

[1] Budarapu P R, Gracie R, Yang S W, Zhuang X, Rabczuk T. Efficient coarse graining in multiscale modeling of fracture. Theoretical and Applied Fracture Mechanics, 2014, 69: 126–143

[2] Talebi H, Silani M, Bordas S P A, Kerfriden P, Rabczuk T. A computational library for multiscale modeling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[3] Griffith A A. The phenomena of rupture and flow in solids. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1921, 221(582–593): 163–198

[4] Irwin G R. Analysis of stresses and strains near the end of a crack traversing a plate. Journal of Applied Mechanics, 1957, 24: 361–364

[5] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131–150

[6] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[7] Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Computational Mechanics, 2007, 39(6): 743–760

[8] Rabczuk T, Zi G, Gerstenberger A, Wall W A. A new crack tip element for the phantom-node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

[9] Rabczuk T, Areias P M A, Belytschko T. A meshfree thin shell method for non-linear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[10] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[11] Remmers J J C, de Borst R, Needleman A. The simulation of dynamic crack propagation using the cohesive segments method. Journal of the Mechanics and Physics of Solids, 2008, 56(1): 70–92

[12] Thai T Q, Rabczuk T, Bazilevs Y, Meschke G. A higher-order stress-based gradient-enhanced damage model based on isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2016, 304: 584–604

[13] Ghorashi S S, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based XIGA for fracture analysis of orthotropic media. Computers & Structures, 2015, 147: 138–146

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

[15] Miehe C, Gürses E. A robust algorithm for configurational-force-driven brittle crack propagation with R-adaptive mesh alignment. International Journal for Numerical Methods in Engineering, 2007, 72(2): 127–155

[16] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotation. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

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

[18] Bourdin B, Francfort G A, Marigo J J. The variational approach to fracture. Journal of Elasticity, 1998, 91(1–3): 5–148

[19] Hakim V, Karma A. Laws of crack motion and phase-field models of fracture. Journal of the Mechanics and Physics of Solids, 2009, 57(2): 342–368

[20] Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45–48): 2765–2778

[21] Amiri F, Millan D, Arroyo M, Silani M, Rabczuk T. Fourth order phase-field model for local max-ent approximants applied to crack propagation. Computer Methods in Applied Mechanics and Engineering, 2016, 312(C): 254–275

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

[23] Areias P, Rabczuk T. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elements in Analysis and Design, 2017, 132: 27–41

[24] Areias P, Rabczuk T, de Sá J C. A novel two-stage discrete crack method based on the screened Poisson equation and local mesh refinement. Computational Mechanics, 2016, 58(6): 1003–1018

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

[26] Areias P, Reinoso J, Camanho P, Rabczuk T. A constitutive-based element-by-element crack propagation algorithm with local remeshing. Computational Mechanics, 2015, 56(2): 291–315

[27] Areias P M A, Rabczuk T, Camanho P P. Finite strain fracture of 2D problems with injected anisotropic softening elements. Theoretical and Applied Fracture Mechanics, 2014, 72: 50–63

[28] Areias P, Rabczuk T, Dias-da-Costa D. Element-wise fracture algorithm based on rotation of edges. Engineering Fracture Mechanics, 2013, 110: 113–137

[29] Areias P, Rabczuk T, Camanho P P. Initially rigid cohesive laws and fracture based on edge rotations. Computational Mechanics, 2013, 52(4): 931–947

[30] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69: 102–109

[31] Areias P, Rabczuk T, Msekh M. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312(C): 322–350

[32] Msekh M A, Nguyen-Cuong H, Zi G, Areias P, Zhuang X, Rabczuk T. Fracture properties prediction of clay/epoxy nanocomposites with interphase zones using a phase field model. Engineering Fracture Mechanics, 2017

[33] Msekh M A, Silani M, Jamshidian M, Areias P, Zhuang X, Zi G, He P, Rabczuk T. Predictions of J integral and tensile strength of clay/epoxy nanocomposites material using phase field model. Composites. Part B, Engineering, 2016, 93: 97–114

[34] Hamdia K, Msekh M A, Silani M, Vu-Bac N, Zhuang X, Nguyen-Thoi T, Rabczuk T. Uncertainty quantification of the fracture properties of polymeric nanocomposites based on phase field modeling. Composite Structures, 2015, 133: 1177–1190

[35] Ambati M, Gerasimov T, De Lorenzis L. Phase-field modeling of ductile fracture. Computer Methods in Applied Mechanics and Engineering, 2015, 55: 1017–1040

[36] Areias P, Dias-da-Costa D, Sargado J M, Rabczuk T. Element-wise algorithm for modeling ductile fracture with the Rousselier yield function. Computational Mechanics, 2013, 52(6): 1429–1443

[37] Mauthe S, Miehe C. Hydraulic fracture in poro-hydro-elastic media. Mechanics Research Communications, 2017, 80: 69–83

[38] Franke M, Hesch C, Dittmann M. Phase-field approach to fracture for finite deformation contact problems. Proceedings in Applied Mathematics and Mechanics, 2016, 16(1): 123–124

[39] de Souza Neto E A, Petric D, Owen D R J. Computational Methods for Plasticity: Theory and Applications. Chichester: Wiley, 2008

[40] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[41] Singh N, Verhoosel C, de Borst R, van Brummelen E. A fracture-controlled path-following technique for phase-field modeling of brittle fracture. Finite Elements in Analysis and Design, 2016, 113: 14–29

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature