A FEniCS implementation of the phase field method for quasi-static brittle fracture

Abstract

In the recent years, the phase field method for simulating fracture problems has received considerable attention. This is due to the salient features of the method: 1) it can be incorporated into any conventional finite element software; 2) has a scalar damage variable is used to represent the discontinuous surface implicitly and 3) the crack initiation and subsequent propagation and branching are treated with less complexity. Within this framework, the linear momentum equations are coupled with the diffusion type equation, which describes the evolution of the damage variable. The coupled nonlinear system of partial differential equations are solved in a ‘staggered’ approach. The present work discusses the implementation of the phase field method for brittle fracture within the open-source finite element software, FEniCS. The FEniCS provides a framework for the automated solutions of the partial differential equations. The details of the implementation which forms the core of the analysis are presented. The implementation is validated by solving a few benchmark problems and comparing the results with the open literature.

Keywords

phase field method / FEniCS / brittle fracture / crack propagation / variational theory of fracture

Introduction

Fracture is one of the most predominant modes of failures in the materials. The estimation of various fracture parameters, viz., the failure strength, the crack initiation and subsequent crack propagation and its trajectory, the crack branching and coalescence are of prime importance for a reliable design and for the safe operation of the structure. Hence, understanding the failure of materials due to the externally applied load in the presence of discontinuities is of great importance to the practicing engineers and scientists. Experimental investigations of the fracture parameters are cumbersome: 1) due to the complex nature of the fracture process and/or 2) due to the presence of multi-physical interactions.

To aid in the understanding and to simulate fracture process, researchers rely on numerical methods. The numerical methods could be boundary based, such as the boundary element method [¹,²] and the scaled boundary finite element method (SBFEM) [³] or domain based. The domain-based methods can be further classified as mesh based, for example, the finite element method or meshless methods [⁴,⁵]. Of these, the FEM is industry accepted, widely used method. This can be attributed to its robustness and the flexibility it offers to treat complex geometries and boundary conditions. Within the framework of the FEM, there are two approaches to represent the discontinuity: discrete [⁶–⁹] and diffuse or smeared crack models [¹⁰–¹²]. In the discrete model, the crack is modeled as a geometric discontinuity. The discrete model can further be classified into implicit or explicit depending on the choice of representation of the crack within the domain. In explicit discrete model [⁶], the crack is modeled via the displacement-discontinuity in an interface element that separates the two solid elements. The accuracy of the explicit discrete models could further be enhanced by combining it with strain smoothing technique [¹³–¹⁸]. Although robust, the explicit model requires the crack topology to conform to the underlying discretization. This increases the computational burden when the discontinuity evolves in time. On another front, within the implicit discrete model [⁷,⁸], the crack or other discontinuities can be represented independent of the underlying FE discretization, thanks to the partition of unity method. One of the robust methods, which also has been implemented now in commercial software such as Abaqus [¹⁹] is the eXtended Finite Element Method (XFEM) [²⁰,²¹]. Meshfree methods are also a better alternative to the explicit discrete models, which can be combined with cohesive zone model for the crack initiation in the concrete failure [²²].

In the smeared crack model, the cracked domain is considered to be a continuum and the description of the discontinuity is represented in terms of the degradation of the stress-strain relations. Recently, the variational approach to fracture has gained popularity [²³–²⁵], wherein the fracture process is modeled as the minimizer of a global energy functional. The method avoids explicit tracking of the interfaces by making them spatially diffuse by introducing a scalar field variable. The scalar field variable has two roles: one it distinguishes different states of the material, viz., the damaged and the undamaged state and second, it aids in automatically tracking the crack trajectories by an evolution equation [²⁶]. The Phase Field Methods (PFM) was originally developed for phase transformations in the material [²⁷–³⁰] and in the last decade, the PFM has been applied to fracture problems starting from the seminal work of Francfort and Marigo [³¹]. Bourdin et al. [²⁶], regularized the formulation to approximate the strong discontinuities by considering a scalar field variable. Since then, the PFM has gained notable attention and has been applied to wide variety of fracture problems, viz., dynamic brittle fracture [³²,³³], fracture in biological tissues [³⁴], ductile fracture [³⁵–³⁷], thermal fracture [³⁸,³⁹] to name a few. In the recent years, an emphasis has been given to implement this method in commercial software Abaqus/Standard through UEL (User ELement) [⁴⁰–⁴²].

The salient feature of the phase field approach for fracture stems from the fact that it does not require additional criteria for the crack initiation, propagation, branching, and coalescence. Moreover, it can handle complex crack trajectories even in three dimensions [⁴³] without additional modifications. The main objective of the present work is (i) to implement the two variants of the PFM, viz., the isotropic model and the hybrid model [⁴⁴] in the open source FEniCS (Finite Element Computational Software) [⁴⁵] to simulate the quasi-static brittle fracture in an elastic material in two and three dimensions, (ii) to implement PFM for crack growth/deflection near the material interfaces.

FEniCS is a framework for the automated solutions of the partial differential equations (PDEs) based on the FEM [⁴⁶]. The user has to specify the finite element variational form with the corresponding geometry and the mesh information to the software. FEniCS software takes care of the element stiffness matrix calculations and assembly to obtain the global stiffness matrix. Another advantage of the FEniCS is the relative ease in extending the two-dimensional analysis to three-dimensional analysis. This feature will be exploited for solving three-dimensional fracture problems.

The paper is structured as follows: Section 2 presents an overview and the governing equations for the phase field method applied to linear elastic fracture mechanics. The corresponding weak form, the spatial discretization and the FEniCS implementation details are presented in Section 3. Section 4 presents the solution of a few boundary value problems to validate the implementation of the phase field method. Major conclusions are discussed in the last section.

Overview of the phase field method for brittle fracture

Within the framework of the phase field method, the strong discontinuity, viz., crack is modeled by a diffused field variable \varphi [0,1]. This field variable differentiates between multiple phases through a smooth transition, with \varphi [0,1] representing the undamaged and the fully damaged state of the material. The choice of the function \varphi is arbitrary, however, it should satisfy the following conditions: 1) symmetric along the crack face and 2) monotonically decay to zero from the crack face to the boundary. A typical choice of the function in case of one dimension is shown in Fig. 1. In this section, important equations pertaining to the modeling of the strong discontinuity using the PFM are presented. For more detailed discussion, interested readers can refer to the literature [⁴⁰,⁴⁴,⁴⁷].

Consider linear elasto-static body with a discontinuity occupying the domain Ω \subset {\mathrm{ℝ}}^d, where d = 2,3. The boundary is considered to admit the decomposition of the boundary \mathit{\Gamma }with outward normal n into three disjoint sets, i.e., \mathit{\Gamma }={\mathit{\Gamma }}_D\cup {\mathit{\Gamma }}_N\cup {\mathit{\Gamma }}_C and {\mathit{\Gamma }}_D\cap {\mathit{\Gamma }}_N=\varnothing. The closure of the domain is \overline{\mathbf{\Omega }}\equiv \mathbf{\Omega }\cup \mathit{\Gamma }. Let (\mathit{u},\mathit{\varphi }):\mathbf{\Omega }\to {\mathrm{ℝ}}^d be the displacement field and the scalar damage variable at a point x of the elastic body when subjected to external tractions \bar{t} on {\mathit{\Gamma }}_N and imposed essential boundary conditions (\widehat{\mathit{u}},\widehat{\mathit{\varphi }}) on {\mathit{\Gamma }}_D. The boundary value problem in the absence of body force is then: find (\mathit{u},\mathit{\varphi }):\mathbf{\Omega }\to {\mathrm{ℝ}}^dsuch that

(1a)[{(\mathrm{1}-\varphi )}^{\mathrm{2}}+k]\nabla \cdot \mathit{\sigma }=0\mathrm{in}\mathrm{\Omega ，}

(1b)-G_c{\ell }_o{\nabla }^{\mathrm{2}}\varphi +\left[{\displaystyle \frac{G_c}{{\ell }_o}}+\mathrm{2}H\right]\varphi =\mathrm{2}H\mathrm{in}\mathrm{\Omega },

with the following boundary conditions:

[{(\mathrm{1}-\varphi )}^{\mathrm{2}}+k]\nabla \cdot \mathit{\sigma }=\overline{\mathit{t}}\mathrm{on}{{\displaystyle \mathit{\Gamma }}}_{\mathrm{N}},

\mathit{u}=\overline{\mathit{u}}\mathrm{o}\mathrm{n}{{\displaystyle \mathit{\Gamma }}}_{\mathrm{D}},

(2)\nabla \varphi \cdot \mathit{n}=\mathrm{0}\mathrm{on}{{\displaystyle \mathit{\Gamma }}}_{\mathrm{N}},

where \mathit{\sigma }={\displaystyle \frac{\partial \psi (\epsilon )}{\partial \epsilon }} is the Cauchy stress tensor and \psi (\epsilon )={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\lambda {(\mathrm{tr(}\epsilon ))}^{\mathrm{2}}+\mu \mathrm{tr(}{\epsilon }^{\mathrm{2}}) is the elastic energy density with \mathit{\mu } and \lambda the Lamé constants, \epsilon ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}[\nabla u+\nabla u^{\mathrm{T}}] is the small strain tensor and the history variable H is defined as:

(3)

where H_n is the strain energy computed at load step n. The above model is referred to as an isotropic model in the literature [⁴⁴]. Although this model captures the discontinuity propagation, it fails to distinguish between the growth of the discontinuity when the local stress state is tension and compression. In order to prevent the propagation of the discontinuity in compressive region, various approaches have been proposed in the literature: 1) additive decomposition of the elastic energy density function, viz., decomposition based on volumetric and deviatoric [⁴⁸] and spectral decomposition [⁴⁷] and 2) hybrid formulation [⁴⁴]. The first approach leads to a non-linear system of equations for the elasticity as a result of the decomposition of the strain tensor. This could be computationally expensive. In order to alleviate this, Ambati et al. [⁴⁴], introduced a hybrid approach. In the hybrid approach, the linear momentum equation is retained same as the isotropic model (see Eq. (1a)) and the evolution of the phase field variable is controlled by the tensile elastic energy {\psi }^{+}. This ensures that the crack surfaces do not grow in the regions of compression and moreover the linear momentum equation is linear which reduces the computational complexity. The modified phase field equation is given by:

(4)-G_c{\ell }_o{\nabla }^{\mathrm{2}}\varphi +\left[{\displaystyle \frac{G_c}{{\ell }_o}}+\mathrm{2}{H}^{+}\right]\varphi =\mathrm{2}{H}^{+}\mathrm{in}\mathrm{\Omega },

where H^{+}={}_{\tau \in [\mathrm{0},t]}{}^{\max }\psi {}^{+}(\mathit{\epsilon }(\mathit{x},\tau )). Eq. (4) is supplemented with the following constraint to prevent crack face inter penetration:

(5)

and

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

with

(7){\mathit{\epsilon }}_{\pm }={\displaystyle \sum _{I=\mathrm{1}}^{\mathrm{3}}{\langle {\epsilon }_I\rangle }_{\pm }{\mathit{n}}_I\otimes {\mathit{n}}_I}\mathrm{and}\epsilon ={\displaystyle \sum _{I=\mathrm{1}}^{\mathrm{3}}\langle {\epsilon }_I\rangle {\mathit{n}}_I\otimes {\mathit{n}}_I},

where {\epsilon }_I and {\mathit{n}}_I are the principle strains and principle strain directions, respectively.

The weak form and an overview of FEniCS implementation

In this section, we develop the weak form of the governing equations (see Eq. (1)). Let us define the infinite dimensional trial (\mathcal{U},\mathcal{P}) and test spaces ({\mathcal{V}}^{\mathrm{0}},\mathcal{Q}). Let \mathcal{W}(\mathrm{\Omega })include the linear displacement field and phase field variable:

(\mathcal{P},{\mathcal{Q}}^0):=\left\{(\phi ,q)\in {\left[C^0\left(\mathrm{\Omega }\right)\right]}^d:\right(\phi ,q)\in {\left[\mathcal{W}\left(\mathrm{\Omega }\right)\right]}^d\subseteq {\left[H^1\left(\mathrm{\Omega }\right)\right]}^d,(\phi ,q)=0\mathrm{o}\mathrm{n}{\mathit{\Gamma }}_{\phi }\}\mathrm{.}

Upon applying the standard Bubnov-Galerkin procedure, we obtain the following weak form: Find \mathit{u}\in \mathcal{U} & \phi \in \mathcal{P}such that, for all \mathit{v}\in {\mathcal{V}}^{\mathrm{0}}& \theta \in \mathcal{Q},

(8a)a(\mathit{u},\mathit{v})=l(\mathit{v})，

(8b)b(\varphi ,\theta )=l(\theta ),

where

a(\mathit{u},\mathit{v})={\displaystyle {\int }_{\mathrm{\Omega }}[{(\mathrm{1}-\varphi )}^{\mathrm{2}}+k]\mathit{\sigma }(\mathit{u}):\mathit{\epsilon }(\mathit{v})}\mathrm{d}\mathrm{\Omega }

l(\mathit{v})={\displaystyle {\int }_{\mathrm{\Omega }}\mathit{b}\cdot \mathit{v}\mathrm{d}\mathrm{\Omega }+{\displaystyle {\int }_{{\Gamma }_t}\overline{\mathit{t}}\cdot \mathit{v}\mathrm{d}\Gamma ,}}

b(\varphi ,\theta )={\displaystyle \underset{\mathrm{\Omega }}{\int }\left\{\nabla \theta G_c{\ell }_o\nabla \varphi +\theta \left[{\displaystyle \frac{G_c}{{\ell }_o}}+\mathrm{2}{H}^{+}\right]\varphi \right\}\mathrm{d}\mathrm{\Omega }},

(9)l(\theta )={\displaystyle \underset{\mathrm{\Omega }}{\int }\mathrm{2}{H}^{+}\theta \mathrm{d\Omega +}{\displaystyle \underset{\Gamma }{\int }\nabla \varphi \cdot \mathit{n}\theta \mathrm{d}\Gamma }}.

In this work, Eq. (8) is solved using a staggered solver, wherein, we solve for the phase field (\varphi) using the displacement fields (u). Then, the updated phase field is used to solve for the displacement field. This staggered process is depicted in Fig. 2. The salient feature of FEniCS is that it takes the weak form (see Eq. (8)) directly. The finite element discretization, the assembly and the solution of the resulting system of equations are done automatically. The weak form and the corresponding finite element discretization are specified using a domain-specific language, called UFL (Unified Form Language) embedded in Python. The domain is partitioned into non-overlapping elements and the discrete trial and test spaces are defined by constructing finite element shape functions over the union of all elements in \mathbf{\Omega }. This is done by the following instruction to FEniCS:

V\mathrm{=\; FunctionSpace(mesh,\text{'}}\mathit{P}\mathrm{\text{'},}\mathit{m}\mathrm{),}

where ‘mesh’ consists of cells which in two dimensions are triangles and tetrahedra in three dimensions, 'P'defines the Lagrange family of elements and ‘m’ denotes the order of the polynomial. Listings 1‒3 (please see the Appendix) gives the Python code for the staggered solution scheme to be used within FEniCS.

Numerical examples

In this section, we present the results concerning the FEniCS implementation of the PFM for brittle fracture through various boundary value problems available in the literature. These fracture problems are chosen to show various characteristics of cracks, such as crack propagation in straight and tracing curvilinear crack paths, crack kinking, crack merging towards a geometric singularity (hole), crack initiation at a notch and bi-material interface. The characteristic length {\ell }_oin all the problems is chosen as {\ell }_o=\mathrm{2}h or otherwise, stated, where h is the size of the smallest element in the domain. Plane strain condition has been considered in all the two-dimensional examples.

Benchmark problem (uni-axial tension applied to a plane strain domain): domain discretized with two triangular elements

To validate the FEniCS implementation of PFM for brittle fracture, we start with a simplest example where the analytical solutions are available. We consider a two-dimensional plain strain domain subjected to uni-axial tension as shown in Fig. 3(a). Since FEniCS supports only triangular elements (the quad element is under development), the plate is discretized with two triangular elements as shown in Fig. 3(b). Stress state in the chosen problem is uniform and hence linear triangular element predicts the exact solution. The material properties and phase field parameters are chosen as: Young’s modulus E = 210 kN/mm², Poisson’s ratio n = 0.3, critical energy release rate G_c = 5×10⁻³ kN/mm and characteristic length {\ell }_o= 0.1 mm. The displacement is applied at the top edge with the increment of Du=1×10⁻³ mm up to u= 0.1 mm. The analytical solutions of axial stress, elastic strain energy, and phase-field parameter are given as [⁴²]:

(10a){\sigma }_{y,\mathrm{0}}={\displaystyle \frac{E(\mathrm{1}-\nu )}{(\mathrm{1}+\nu )(\mathrm{1}-\mathrm{2}\nu )}}{\epsilon }_y,

(10b)\psi ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{E(\mathrm{1}-\nu )}{(\mathrm{1}+\nu )(\mathrm{1}-\mathrm{2}\nu )}}{\epsilon }_y^{\mathrm{2}},

(10c)\varphi ={\displaystyle \frac{\mathrm{2}\Psi }{{\displaystyle \frac{G_c}{{\ell }_o}}+\mathrm{2}\Psi }}={\displaystyle \frac{{\displaystyle \frac{E(\mathrm{1}-\nu )}{(\mathrm{1}+\nu )(\mathrm{1}-\mathrm{2}\nu )}}{\epsilon }_y^{\mathrm{2}}}{{\displaystyle \frac{G_c}{{\ell }_o}}+{\displaystyle \frac{E(\mathrm{1}-\nu )}{(\mathrm{1}+\nu )(\mathrm{1}-\mathrm{2}\nu )}}{\epsilon }_y^{\mathrm{2}}}},

(10d){\sigma }_y={\sigma }_{y,\mathrm{0}}{(\mathrm{1}-\varphi )}^{\mathrm{2}},

where {\epsilon }_yis the axial strain and {\sigma }_y is the axial stress with considering damage.

Figure 4(a), shows the evolution of the phase-field parameter as a function of applied axial strain. Phase-field parameter \varphi is evolving from zero initially (i.e., no damage) with the applied axial strain and reaches a maximum value of one (i.e., fully damaged state).

Figure 4(b), shows the evolution of axial stress as a function of applied axial strain. It can be observed that the load-bearing capacity of the material (stiffness of the material) degrades with the evolution of the phase-field parameter (\varphi). The results show an excellent agreement between the analytical and the FEniCS solutions.

Three-point bend test

Here, we consider a simply supported three-point bend test to verify the implementation with the recent literature [⁴⁴]. Figure 5(a) show the domain and the boundary conditions for the boundary value problem chosen. The material properties are chosen as l = 12 kN/mm², m = 8 kN/mm² and G_c = 5.4×10⁻⁴ kN/mm. The mesh is refined in the region where the crack is expected to propagate (see shown in Fig. 5(b)). The displacement at the top is applied incrementally, with an increment Du = 1×10⁻⁴ mm up to u = 4×10⁻² mm and subsequently Du = 1×10⁻⁵ mm till the complete failure of the specimen. Figure 6 shows the crack initiation at the notch and further propagation at the various stages of loading. The red line shows the complete failure of the material (\varphi=1) and blue region (\varphi=0) shows the intact state of the material. Note that there is a smooth transition from blue to red color, which describes the diffuse crack as mentioned in Section 2. The crack initiates at the notch and propagates towards the maximum tensile stress regime. The load-displacement curve (Fig. 7(a)) shows a sudden drop after the peak stress. Since the simulation uses a displacement controlled approach it is not possible to predict the snap-back after the peak stress as observed in Ref. [⁴⁹]. The choice of length scale parameter {\ell }_oslightly influences the maximum reaction force as shown in Fig. 7(b).

Plate with an edge crack subjected to tension and shear load

We consider a square plate with a straight edge crack. Firstly, the specimen is subjected to displacement at the top in y-direction to simulate a pure mode-I fracture. The geometry and boundary conditions are shown in Fig. 8(a). We present the results for three different domain discretizations: 32,768 ({\ell }_o= 0.022), 131,072 ({\ell }_o= 0.011), and 524,288 ({\ell }_o= 0.0055) triangular elements, respectively without any localized refinement. The material properties are chosen as l = 121.15 kN/mm², m = 80 kN/mm², G_c = 2.7 kN/mm. Displacement is applied with an increment of Du = 1×10⁻⁵ mm up to u = 5×10⁻³ mm and Du = 1×10⁻⁶ mm up to the failure of the specimen. Figure 9 shows the crack propagation at successive stages of applied displacement. The crack propagates perpendicular to the applied external displacement. Figure 10(a) shows the load-displacement curves with different characteristic length {\ell }_o. It is inferred that the maximum reaction approaches to that of Ref. [⁴⁴] as the characteristic length scale {\ell }_odecreases as shown in Fig. 10(b).

Now the same specimen is loaded in shear mode. The domain and boundary conditions are shown in Fig. 8(b). The domain is discretized with 131,072 ({\ell }_o= 0.011) triangular elements uniformly and the load increment of Du = 1×10⁻⁵ mm is applied. In the shear mode, portion above the crack is in a compressive state of stress whereas portion below the crack is in the tensile state of stress. Hence, crack is expected to propagate through the portion below the crack, i.e., in the tensile regime but two symmetric branches of crack start to propagate from the crack tip as shown in Fig. 11(a) while using without energy decomposition model. This unrealistic crack propagation is not physical (see also reported in Refs. [⁴⁷,⁵⁰]). Model without tension-compression split does not distinguish between positive energy and negative energy (as discussed in Section 3). In this case, the magnitude of positive and negative energy is same and hence two symmetric branches of crack propagate. Unphysical crack propagation in a compressive state of stress is avoided using the hybrid model as shown in Fig. 11(b), where crack propagates only due to tensile part of the energy. Figure 12 shows the load-displacement curve for both isotropic and hybrid model and compared with [⁴⁴] hybrid model.

Cruciform plate with a corner crack

The geometry and boundary conditions are shown in Fig. 13(a). The ratio of a/L is taken as 0.2 and the crack is inclined 45° from the horizontal axis [⁵¹]. The material properties are chosen as: E = 21.84 kN/mm², n = 0.3, G_c = 2.7×10⁻³ kN/mm. The domain is discretized with 143,606 triangular elements with {\ell }_o= 0.0083. The displacement is applied at the top with an increment of Du = 1×10⁻³ mm up to u= 0.012 mm and to show the effect of maximum reaction force on loading Du = 1×10⁻⁵ and 1×10⁻⁶ is applied till the failure of the specimen. The crack propagation path has been shown in Fig. 13(b), which matches the predictions of [⁵¹,⁵²]. Figure 14 shows the load-displacement curve with different magnitude of Du applied after u = 0.012. The maximum reaction force decreases for the slow loading, i.e., Du = 1×10⁻⁶.

Notched plate with a hole

Now, we consider a notched plate with a hole. The center of the hole is kept offset from the median axis of the plate to induce mixed-mode fracture. The domain and applied boundary conditions are shown in Fig. 15. The domain is discretized with 100,395 triangular elements with characteristic length {\ell }_o= 0.45. The mesh is refined where the crack is expected to propagate (i.e., near the hole and notch). The material parameters are chosen as: l = 1.94 kN/mm², m = 2.45 kN/mm², G_c = 2.28×10⁻³ kN/mm. The displacement is applied with the increment of Du = 1×10⁻³ mm up to 0.38 mm and Du = 1×10⁻⁴ mm till complete failure of the specimen. The crack propagation with successive intervals is shown in Fig. 16, which shows that initially a straight crack propagates from the notch and due to the mode-mixity, the crack kinks and traverses along a curved path prior to merging into the hole. A secondary crack is then initiated from the hole and propagates straight till the edge of the specimen. The numerical simulation captures experimentally observed results shown in Refs. [²⁰,⁴⁴]. Figure 17 shows the load-displacement curve, in which the kink in the softening regime indicates the delayed fracture due to the merger of the crack into the hole.

Interaction of a crack with an elastic inclusion

In this example, we consider a bi-material plate with a crack to study the crack-inclusion interactions [⁵³]. The geometry and boundary conditions are shown in Fig. 18(a). Two different ratios of Young’s modulus (E_r = E_matrix/E_inclusion) (a) E_r = 0.1 and (b) E_r = 10, are chosen to study the crack growth in the presence of hard and soft inclusions. The domain is discretized with 137,804 triangular elements with {\ell }_o= 0.0278. The mesh is refined near the inclusion and notch. Crack propagation direction from the notch is highly dependent on the type of inclusion. Crack is expected to propagate straight from the crack for the case of homogeneous material (pure mode-I). On contrary here hard inclusion deflects the crack as shown in Fig. 18(b) whereas the soft inclusion attracts the crack as shown in Fig. 18(c) from its homogeneous material crack. Figure 19 shows the load-displacement curve for the hard and soft inclusions. The maximum reaction force for the soft inclusion is more than hard inclusion.

Three-dimensional single notched plate

To demonstrate the simplicity of the extension from the two-dimensional models to three-dimensional model in FEniCS, we consider a pure mode-I tension test in three dimensions. The domain and boundary conditions are shown in Fig. 20. The domain is discretized with 27,066 ({\ell }_o= 0.537), 41,245 ({\ell }_o= 0.342), 134,567 ({\ell }_o= 0.157) tetrahedron elements. The mesh is refined in the region where the crack is expected to propagate. The material parameters are chosen as E = 20.8 kN/mm², n = 0.3 and G_c = 5×10⁻⁴ kN/mm. The displacement is applied with an increment of Du =1×10⁻⁴ mm up to u = 0.25 mm and Du = 1×10⁻⁵ mm up to failure. The crack starts to propagate in the direction perpendicular to the displacement applied. The growth of the crack is shown in Fig. 21. The iso-surface of the phase-field parameter (\varphi) is shown for the visualization of the crack growth. Figure 22 shows the load-displacement curve with different characteristic length {\ell }_o. The maximum reaction force shows the similar behavior as we have observed in the two-dimensional tension test Section 4.3.

Cube with a rigid spherical inclusion

In this final example, we consider a cube with a rigid spherical inclusion as shown in Fig. 23. Tension is applied on top (ABCD) and bottom (EFGH) faces. The main objective of this problem is to simulate the crack initiation at the bimaterial interface in the absence of any geometric discontinuities such as a crack. The material parameters are kept the same as in Section 4.7. Taking the advantage of the symmetry, only one quarter of the problem is analyzed numerically. The domain is discretized with 140,142 tetrahedron elements with {\ell }_o= 0.0387. The mesh is refined near the bimaterial interface. The displacement is applied with an increment of Du =1×10⁻⁴ mm. Once the stress state reaches a critical value, the crack initiates at the top of the inclusion surface as shown in Fig. 24(a). Then, the crack continues to propagate horizontally until the complete failure of the material as shown in Fig. 24(b). Figure 25 shows the corresponding load-displacement curve.

Conclusions

In this work, we presented the procedure to implement the staggered phase field model within an open-source finite element package FEniCS for the two-dimensional as well as three-dimensional brittle fracture problems. The implementation is validated with the results available in the open literature. The effect of characteristic length scale ({\ell }_o) on the maximum reaction force is also investigated in the two-dimensional as well as three-dimensional fracture problems. The results show that the choice of the length scale and the rate of loading are important to obtain the accurate results.

The FEniCS codes of various benchmark problems discussed in this work can be downloaded from https: //home.iitm.ac.in/ratna/codes/phasefield. These codes can be a starting point for the material scientists and engineers to simulate various fracture problems of interest in their domain.

References

[1] Trevelyan J. Boundary Elements for Engineers: Theory and Applications, vol. 1. Southampton: Computational Mechanics, 1994

[2] Aliabadi M H. A new generation of boundary element methods in fracture mechanics. International Journal of Fracture, 1997, 86(1–2): 91–125

[3] Ooi E T, Song C, Tin-Loi F, Yang Z. Polygon scaled boundary finite elements for crack propagation modelling. International Journal for Numerical Methods in Engineering, 2012, 91(3): 319–342

[4] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

[5] Nguyen V P, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3): 763–813

[6] Dechaumphai P, Phongthanapanich S, Sricharoenchai T. Combined Delaunay triangulation and adaptive finite element method for crack growth analysis. Acta Mechanica Sinica, 2003, 19(2): 162–171

[7] 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

[8] Wang C, Xu X. An extended phantom node method study of crack propagation of composites under fatigue loading. Composite Structures, 2016, 154: 410–418

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

[10] Aranson I S, Kalatsky V A, Vinokur V M. Continuum field description of crack propagation. Physical Review Letters, 2000, 85(1): 118–121

[11] Marigo J J, Maurini C, Pham K. An overview of the modelling of fracture by gradient damage models. Meccanica, 2016, 51(12): 3107–3128

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

[13] Bordas S P A, Rabczuk T, Hung N X, Nguyen V P, Natarajan S, Bog T, Quan D M, Hiep N V. Strain smoothing in FEM and XFEM. Computers & Structures, 2010, 88(23‒24): 1419–1443

[14] Nguyen-Xuan H, Liu G R, Nourbakhshnia N, Chen L. A novel singular ES-FEM for crack growth simulation. Engineering Fracture Mechanics, 2012, 84: 41–66

[15] Nguyen-Xuan H, Liu G R, Bordas S, Natarajan S, Rabczuk T.An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order. Computer Methods Applied Mechanics and Engineering, 2013, 253: 252–273

[16] Natarajan S, Bordas S P A, Ooi E T. Virtual and smoothed finite elements: a connection and its application to polygonal/polyhedral finite element methods. International Journal for Numerical Methods in Engineering, 2015, 104(13): 1173–1199

[17] Francis A, Ortiz-Bernardin A, Bordas S P A, Natarajan S. Linear smoothed polygonal and polyhedral finite elements. International Journal for Numerical Methods in Engineering, 2017, 109(9): 1263–1288

[18] Surendran M, Natarajan S, Bordas S P A, Palani G S. Linear smoothed extended finite element method. International Journal for Numerical Methods in Engineering, 2017, 112(12): 1733–1749

[19] Abaqus. Abaqus documentation. Dassault Systmes Simulia Corp. , Provid. RI, USA, 2012

[20] Giner E, Sukumar N, Tarancón J E, Fuenmayor F J. An Abaqus implementation of the extended finite element method. Engineering Fracture Mechanics, 2009, 76(3): 347–368

[21] Shi J, Chopp D, Lua J, Sukumar N, Belytschko T. Abaqus implementation of extended finite element method using a level set representation for three-dimensional fatigue crack growth and life predictions. Engineering Fracture Mechanics, 2010, 77(14): 2840–2863

[22] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A geometrically non-linear three-dimensional cohesive crack method for reinforced concrete structures. Engineering Fracture Mechanics, 2008, 75(16): 4740–4758

[23] Bourdin B, Francfort G A, Marigo J J. The Variational Approach to Fracture. Springer Netherlands, 2008

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

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

[26] Bourdin B, Francfort G A, Marigo J J. Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids, 2000, 48(4): 797–826

[27] Cahn J W, Hilliard J E. Free energy of a nonuniform system. I. Interfacial free energy. Journal of Chemical Physics, 1958, 28(2): 258–267

[28] Chen L Q, Khachaturyan A G. Computer simulation of structural transformations during precipitation of an ordered intermetallic phase. Acta Metallurgica et Materialia, 1991, 39(11): 2533–2551

[29] Galenko P K, Herlach D M, Funke O, Gandham P. Phase-field modeling of dendritic solidification: verification for the model predictions with latest experimental data. In: Herlach D M, ed. Solidification and Crystallization. Wiley-VCH, 2005, 52–60

[30] Moelans N, Blanpain B, Wollants P. An introduction to phase-field modeling of microstructure evolution. Calphad, 2008, 32(2): 268–294

[31] 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

[32] Borden M J, Verhoosel C 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: 77–95

[33] Schlüter A, Willenbücher A, Kuhn C, Müller R. Phase field approximation of dynamic brittle fracture. Computational Mechanics, 2014, 54(5): 1141–1161

[34] Raina A, Miehe C. A phase-field model for fracture in biological tissues. Biomechanics and Modeling in Mechanobiology, 2016, 15(3): 479–496

[35] Miehe C, Aldakheel F, Raina A. Phase field modeling of ductile fracture at finite strains: a variational gradient-extended plasticity-damage theory. International Journal of Plasticity, 2016, 84: 1–32

[36] Ambati M, Gerasimov T, De Lorenzis L. Phase-field modeling of ductile fracture. Computational Mechanics, 2015, 55(5): 1017–1040

[37] Ambati M, Kruse R, De Lorenzis L. A phase-field model for ductile fracture at finite strains and its experimental verification. Computational Mechanics, 2016, 57(1): 149–167

[38] Kuhn C, Müller R. Phase field simulation of thermomechanical fracture. Proceedings in Applied Mathematics and Mechanics, 2009, 9(1): 191–192

[39] Schlüter A, Kuhn C, Müller R, Tomut M, Trautmann C, Weick H, Plate C. Phase field modelling of dynamic thermal fracture in the context of irradiation damage. Continuum Mechanics and Thermodynamics, 2017, 29(4): 977–988

[40] Msekh M A, Sargado J M, Jamshidian M, Areias P M, Rabczuk T. Abaqus implementation of phase-field model for brittle fracture. Computational Materials Science, 2015, 96: 472–484

[41] Liu G, Li Q, Msekh M A, Zuo Z. Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field mode. Computational Materials Science, 2016, 121: 35–47

[42] Molnár G, Gravouil A. 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture. Finite Elements in Analysis and Design, 2017, 130: 27–38

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

[44] Ambati M, Gerasimov T, De Lorenzis L. A review on phase-field models of brittle fracture and a new fast hybrid formulation. Computational Mechanics, 2015, 55(2): 383–405

[45] Alnæs M S, Blechta J, Hake H, Hoansson A, Kehlet B, Logg A, Richardson C, Ring J, Rognes M E, Wells G N. The FEniCS Project Version 1.5. Archive of Numerical Software, 2015, 3(100): 9–23

[46] Logg A, Mardal K A, Wells G N. Automated Solution of Differential Equations by the Finite Element Method. Berlin: Springer, 2012

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

[48] Amor H, Marigo J J, Maurini C. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. Journal of the Mechanics and Physics of Solids, 2009, 57(8): 1209–1229

[49] Moës N, Stolz C, Bernard P E, Chevaugeon N. A level set based model for damage growth: the thick level set approach. International Journal for Numerical Methods in Engineering, 2011, 86(3): 358–380

[50] 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): 2765–2778

[51] Prasad N N V, Aliabadi M H, Rooke D P. Incremental crack growth in thermoelastic problems. International Journal of Fracture, 1994, 66(3): R45–R50

[52] Duflot M. The extended finite element method in thermoelastic fracture mechanics. International Journal for Numerical Methods in Engineering, 2008, 74(5): 827–847

[53] Bouchard P O, Bay F, Chastel Y. Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria. Computer Methods in Applied Mechanics and Engineering, 2003, 192(35–36): 3887–3908

RIGHTS & PERMISSIONS

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