Three Dimensional CS-FEM Phase-Field Modeling Technique for Brittle Fracture in Elastic Solids

Abstract

The cell based smoothed finite element method (CS-FEM) was integrated with the phase-field technique to model brittle fracture in 3D elastic solids. The CS-FEM was used to model the mechanics behavior and the phase-field method was used for diffuse fracture modeling technique where the damage in a system was quantified by a scalar variable. The integrated CS-FEM phase-field approach provides an efficient technique to model complex crack topologies in three dimensions. The detailed formulation of our combined method is provided. It was implemented in the commercial software ABAQUS using its user-element (UEL) and user-material (UMAT) subroutines. The coupled system of equations were solved in a staggered fashion using the in-built non-linear Newton–Raphson solver in ABAQUS. Eight node hexahedral (H8) elements with eight smoothing domains were coded in CS-FEM. Several representative numerical examples are presented to demonstrate the capability of the method. We also discuss some of its limitations.

1. Introduction

In the past decade, Liu et al. [1,2] generalized the gradient smoothing approach in meshfree method [3] and proposed the smoothed finite element method(S-FEM) to overcome some of the inherent shortcomings of the classical finite element method(FEM) such as overly stiff behavior, sensitivity to mesh distortions, and stress inaccuracy. The SFEM combines the FEM with the traditional meshfree methods in producing more accurate results with higher efficiency [1,2]. It uses the base mesh of the FEM and reconstructs the strain field using the gradient smoothing technique. Thus, in contrast to the weak formulations of the FEM [4], the weakened weak formulation used here [1,2] further softens the model [5], giving it a closer to exact stiffness [6]. The gradient smoothing operation facilitates creation of various smoothing domains based on elements, nodes, edges (2D), and faces (3D) over which the stress/strain is evaluated [1]. This, in turn, produces a wide variety of results and gives the analyst much needed freedom to design models as per the requirement. For example, the node based smoothing domains (NS-FEM) produce upper bound solutions for force driven problems [7]. The edge/face based smoothing domains (ES/FS-FEM) produce very accurate results using even triangular and tetrahedral mesh [8]. A combination of the NS-FEM and the ES-FEM produces very accurate close to exact solutions in the error norm [6,9]. Amongst other applications [10,11,12,13,14], the S-FEM has been effectively applied to fracture problems for quasi-static crack propagation [15], anisotropic materials [16], dynamic fracture [1] and for singular geometries of arbitrary order [17,18]. Although these studies have proven to be quite accurate and efficient when compared to the FEM, they treat the crack as an discrete entity, a formulation which has its own inherent problems.

Fracture is the primary cause of failure in the majority of engineering structures. An initially existing small crack often leads to catastrophic failures by propagation under external loads. Over the past few decades, several theories and computational techniques have been developed to accurately model fracture and predict crack propagation in engineering systems. Classically, cracks are dealt with discretely in theoretical/computational fracture mechanics. The discontinuities have been separately modeled and special methodologies have been developed to accurately model the singular stress field near the crack tip. The collapsed quadrilateral or triangular element in the FEM [4] is the most widely used, however several other methods such as the boundary element method (BEM) [19], mesh-free methods [3,20], the extended finite element method [21], and SFEM [18,22,23] have proven to be equally efficient. BEM uses a mesh only on the boundaries, thereby reducing the computational complexity by one order, but is unable to treat material non-linearities. S-FEM uses a base mesh of linear elements with an enrichment only around the crack-tip, thereby mitigating the computational cost of using a quadratic mesh throughout the domain [18]. However, in S-FEM, as in the classical FEM, the crack path is mesh dependent; cracks can only propagate along edges of elements and the tracking of crack front for three-dimensional problems is computationally very expensive. The X-FEM overcomes this mesh dependency but uses additional degrees of freedom and the integration in the weak form for the elements containing the crack becomes particularly very complex in three-dimensional fracture problems.

An alternate way of modeling fracture has gained popularity in the computational community in the past decade or two. In this method, the crack is modeled as a diffused entity and represented by a continuous scalar variable called the phase-field [24]. This variable differentiates between the broken and intact material phases and, contrary to the discrete approach, provides a smooth transition between them [25,26,27,28]. This method is based on the variational theory of fracture [24] addressing the shortcomings of the original Griffith’s theory, which could not predict crack nucleation or complex crack paths. The elegance of this method lies in the fact that it can successfully predict complex crack behaviors such as crack branching, curved crack paths, and crack merging even in three dimensions without any ad-hoc criterion. This gives it an immediate advantage over the traditional methods where the prediction of such phenomenons is quite complex. Apart from the primary applications of brittle fracture [28], the phase-field approach has been developed and applied to many complex fracture mechanics problems including large deformations [29], plasticity [30], multiphysics [31,32], and dynamics [33,34]. It has also been previously implemented in the commercial code ABAQUS by using its user subroutine features [27,28]. The phase-field model, however, is not without its disadvantages, including having a very refined mesh in the expected crack propagation region. Advanced fracture algorithms such as the screened Poisson’s equations [35,36,37] for crack propagation are developed to overcome such shortcomings. However, here we limit ourselves to the application of the phase-field method and discuss its useful features as well as disadvantages.

In this work, we integrate the cell based S-FEM with phase-field model to simulate 3-D crack propagation. This integration is done on the commercial platform of ABAQUS. ABAQUS is one of the most widely used commercial codes and its excellent inbuilt solvers and sophisticated visualization tools are particularly attractive for implementing user developed element formulations. A CS-FEM formulation for the eight-node hexahedral elements (3D-CS-FEM-H8) is used, similar to the one used in Xuan et al. [38]. We consider each node to have an additional degree of freedom (phase-field) in addition to the standard displacements. The phase-field and the displacement variables are solved in a staggered fashion to attenuate instability [27]. The 3D CS-FEM-H8 has already been proven to have faster convergence rate and better accuracy than the standard FEM, thus its assimilation with the phase-field method which computes complex crack topologies without any ad-hoc criterion, produces another efficient technique to model crack propagation. One of the significant disadvantages of implementing this method in ABAQUS is that the software does not let users access data for any surrounding elements, thus we deviate from the classical way of calculating strains/stresses in CS-FEM [1,38], and develop another novel, simple yet quite effective approach. It is noteworthy that this is also the reason for our inability to implement S-FEM models of higher accuracy into the software. However, with all these shortcomings, as demonstrated by a number of examples, this method proves to be quite an effective tool to predict complex crack behaviors.

This paper is outlined as follows: In Section 2, we outline the phase-field model of fracture and the 3D-CS-FEM-H8, and discuss the implementation details of the combined method into ABAQUS. In Section 3, we provide numerical examples to validate our method. In Section 4, we conclude the paper with a summary of our findings and proposed future works.

2. Methods

2.1. Governing Equations of a 3D Elastic Solid with Discontinuity

Consider a three-dimensional (3D) arbitrary, homogeneous, linear elastic domain $\mathbf{\Omega }$ bounded by $\mathbf{\Gamma }$ such that $\mathbf{\Gamma }={\mathbf{\Gamma }}_{\mathit{u}}\cup {\mathbf{\Gamma }}_{\mathit{t}}$, ${\mathbf{\Gamma }}_{\mathit{u}}\cap {\mathbf{\Gamma }}_{\mathit{t}}=\mathbf{0}$ and an internal traction free crack ${\mathbf{\Gamma }}_{\mathit{c}}$, as shown in Figure 1, where ${\mathbf{\Gamma }}_{\mathit{u}}$ and ${\mathbf{\Gamma }}_{\mathit{t}}$ are the displacement and traction boundary surfaces, respectively. The equilibrium equations are given as [1,4]:

$$\nabla ·\mathit{\sigma }+{\mathit{f}}^{\mathit{b}}=0$$

(1)

where ∇ is the divergence operator, $\mathit{\sigma }$ is the Cauchy stress tensor and ${\mathit{f}}^{\mathit{b}}$ is the body force.

The Dirichlet and Neumann boundary conditions are given as:

$$\mathit{u}(\mathit{x},t)=\overline{\mathit{u}}(\mathit{x},t)\hspace{1em}\mathrm{on}\hspace{1em}{\mathbf{\Gamma }}_{\mathit{u}}$$

(2)

$$\mathit{\sigma }·\mathit{n}={\mathit{f}}^{\mathit{t}}\hspace{1em}\mathrm{on}\hspace{1em}{\mathbf{\Gamma }}_{\mathit{t}}$$

(3)

$$\mathit{\sigma }·\mathit{n}=\mathbf{0}\hspace{1em}\mathrm{on}\hspace{1em}{\mathbf{\Gamma }}_{\mathit{c}}$$

(4)

where $\mathit{n}$ is the outward unit normal vector on the boundary area $\mathbf{\Gamma }$ and $\overline{\mathit{u}}$ is the prescribed displacement on the boundary ${\mathbf{\Gamma }}_{\mathit{u}}$.

The stress–strain relation is given by the constitutive equation

$$\mathit{\sigma }=\mathit{C}·\mathit{ϵ}$$

(5)

where C is the matrix of elastic constants.

Assuming small displacements and strain, we can define the compatibility relation between strain and displacement as:

$$\mathit{ϵ}={\nabla }^s\mathit{u}\left(\mathit{x}\right)$$

(6)

The elastic strain energy density is given by

$${\psi }_e\left(\mathit{ϵ}\right)=\frac{1}{2}\lambda {\mathit{ϵ}}_{\mathit{ii}}{\mathit{ϵ}}_{\mathit{jj}}+\mu {\mathit{ϵ}}_{\mathit{ij}}{\mathit{ϵ}}_{\mathit{ij}}$$

(7)

where $\lambda$ and $\mu$ are the Lamé constants.

2.2. Review of Phase-Field Model for Brittle Fracture

Following the derivations in [25,27,39], we discuss briefly about the formulations of phase-field approximations for diffuse fracture modeling.

According to the variational theory of fracture, the crack propagates in such a way the total energy of a system is always minimized. Thus, we approach by minimizing the energy functional. The total energy of a discontinuous system is given by

$$\Pi ={\psi }_e\left(\mathit{ϵ}\right)+{\psi }_f-{\psi }_{ext}\left(\mathit{u}\right)$$

(8)

where ${\psi }_e\left(\mathit{ϵ}\right)$ is the elastic strain energy as given in Equation (7) integrated over the entire volume $\Omega$, ${\psi }_f$ is the fracture surface energy and ${\psi }_{ext}\left(\mathit{u}\right)$ is the external potential energy.

The fracture surface energy ${\psi }_f$ is given by

$${\psi }_f=G_c{\int }_{\Omega }\gamma \left(c\right)d\mathbf{\Omega }$$

(9)

where $G_c$ is the critical energy release rate proposed by Griffith, c is the phase-field parameter and $\gamma \left(c\right)$ is the density function which was calculated to be [28]

$$\gamma \left(c\right)=\frac{1}{2}[\frac{1}{l_c}{c}^2+\frac{l_c}{2}{|\nabla c|}^2]$$

(10)

where $l_c$ is the length scale parameter.

The external potential energy is given by the summation of the body force $\mathit{b}$ and the surface traction ${\mathit{t}}_{\mathbf{\Gamma }}$:

$${\Psi }_{ext}\left(\mathit{u}\right)={\int }_{\Omega }\mathit{b}·\mathit{u}d\mathit{V}+{\int }_{\Gamma }{\mathit{t}}_{\mathbf{\Gamma }}·\mathit{u}d\mathit{A}$$

(11)

The total internal potential energy of a system can be written as the summation of the bulk energy and the energy required for the formation of crack [24,28,39]

$$\Psi (\mathit{u},c)={\int }_{\Omega }[{(1-c)}^2+d]\psi \left(\mathit{ϵ}\right)d\mathbf{\Omega }+{\int }_{\Omega }\frac{G_c}{2}[l_c\nabla c·\nabla c+\frac{1}{l_c}{c}^2]d\mathbf{\Omega }$$

(12)

where d is a numerical stabilization parameter.

By variation of the external and internal energy potentials (using Equations (8)–(12)) and thereby imposing the principal of virtual displacements we obtain the governing equations of the model:

$$[{(1-c)}^2+d]\frac{\partial {\mathit{\sigma }}_{\mathit{ij}}}{\partial {\mathit{x}}_{\mathit{i}}}+{\mathit{b}}_{\mathit{j}}=0\hspace{1em}\mathrm{in}\hspace{1em}\mathbf{\Omega }$$

(13)

$$[{(1-c)}^2+d]{\mathit{n}}_{\mathit{i}}{\mathit{\sigma }}_{\mathit{ij}}={\mathit{t}}_{\mathit{j}}\hspace{1em}\mathrm{in}\hspace{1em}{\mathbf{\Gamma }}_{\mathit{t}}$$

(14)

$${\mathit{u}}_{\mathit{j}}={\overline{\mathit{u}}}_{\mathit{j}}\hspace{1em}\mathrm{in}\hspace{1em}{\mathbf{\Gamma }}_{\mathit{u}}$$

(15)

$$-G_c{l}_c\frac{{\partial }^{2}c}{\partial x_i\partial x_i}+[\frac{G_c}{l_c}+2\psi \left(ϵ\right)]c=2\psi \left(ϵ\right)\hspace{1em}\mathrm{in}\hspace{1em}\mathbf{\Omega }$$

(16)

$$\frac{\partial c}{\partial x_i}{n}_i=0\hspace{1em}\mathrm{in}\hspace{1em}\mathbf{\Gamma }$$

(17)

2.3. Three-Dimensional Cell Based Smoothed-Finite Elements

In S-FEM, the idea of strain smoothing is combined with the standard underlying mesh of FEA by dividing each element into a number of sub-cells [1]. We use the eight-noded hexahedral elements with eight smoothing cells (Figure 2). They have already been proven to have a better accuracy than the standard eight-noded hexahedral element used in FEA with eight Gaussian points [38] using the same set of elements. Classical FEM involves isoparametric transformation for every element, making the problem very much mesh dependent because it involves inverting the Jacobian matrix. The solutions tend to vary with minor mesh distortions and for better modeling of certain topological features one needs mesh refinement or higher order elements, which makes the problem computationally expensive. The S-FEM can solve majority of computational problems using a base of linear triangular/tetrahedral mesh, which can be generated automatically. Since the shape functions are created on the basis of radial point interpolation, creating special elements to treat special topological features [18] is of minor hassle and does not need any transition or patch elements. In S-FEM, a smoothing operation is performed on the gradient of the displacement field for each smoothing cell in an element. Subsequently, the interior integration on each of the smoothing cell is transferred to the boundary surface area using the divergence theorem. Herein lies one of the biggest advantages of the S-FEM. This feature makes the solution pretty impervious to mesh distortion since we can avoid the isoparametric transformation and also contributes in saving computational cost. In this case, the gradient smoothing is the spatial average of the strain, over a smoothing hexahedral cell. A smoothing operation is performed to the gradient of displacement for each smoothing cell in an element:

$$\nabla u^h\left({\mathit{x}}_C\right)={\int }_{\Omega }\nabla u^h\left(\mathit{x}\right)\Phi (\mathit{x}-{\mathit{x}}_C)d\mathbf{\Omega }$$

(18)

where ${\mathit{x}}_C$ and $\mathbf{\Omega }$ represent the center and volume of the smoothing domain, respectively, and $\Phi$ is a smoothing function.

Integration by parts on Equation (18) yields:

$$\nabla u^h\left({\mathit{x}}_C\right)={\int }_{\Gamma }{u}^h\left(\mathit{x}\right)\mathit{n}\left(\mathit{x}\right)\Phi (\mathit{x}-{\mathit{x}}_C)d\mathbf{\Gamma }-{\int }_{\Omega }{u}^h\left(\mathit{x}\right)\nabla \Phi (\mathit{x}-{\mathit{x}}_C)d\mathbf{\Omega }$$

(19)

For simplicity, a piecewise constant smoothing function is applied, which is assumed to be constant within ${\Omega }_C$ and vanish everywhere else,

$$\Phi (\mathit{x}-{\mathit{x}}_C)=\left\{\begin{array}{c}{\displaystyle \frac{1}{V_c}}\hspace{1em}\mathrm{for}\hspace{1em}\mathit{x}\in {\mathbf{\Omega }}_{\mathit{c}}\\ 0\hspace{1em}\mathrm{for}\hspace{1em}\mathit{x}\notin {\mathbf{\Omega }}_{\mathit{c}}\end{array}\right.$$

(20)

where $V_C={\int }_{{\mathbf{\Omega }}_{\mathit{c}}}^d\mathbf{\Omega }$ and ${\mathbf{\Omega }}_{\mathit{c}}$ is the smoothing cell.

Substituting Equation (20) into Equation (18), we can get the smoothed gradient or strain field

$$\tilde{\nabla }{u}^h\left({\mathit{x}}_C\right)={\int }_{{\Gamma }_C}{u}^h\left(\mathit{x}\right)\mathit{n}\left(\mathit{x}\right)\Phi (\mathit{x}-{\mathit{x}}_C)d\mathbf{\Gamma }=\frac{1}{V_c}{\int }_{{\Gamma }_c}{u}^h\left(\mathit{x}\right)\mathit{n}\left(\mathit{x}\right)d\mathbf{\Gamma }$$

(21)

with $d{\mathbf{\Gamma }}_{\mathit{c}}$ denoting the boundary of the smoothing cell. It is noteworthy that the choice of the smoothing function makes the second term of Equation (19) vanish and the area integration is thus converted to line integration along the boundary of the smoothing cell.

The smoothed strain can be obtained as

$${\widetilde{ϵ}}^h\left({\mathit{x}}_C\right)=\sum _{I=1}^n{\tilde{B}}_I\left({\mathit{x}}_C\right){\mathit{d}}_{\mathit{I}}$$

(22)

where ${\tilde{B}}_I$ is the smoothed strain gradient matrix given by

$${\tilde{B}}_I=\left[\begin{array}{ccc}{\tilde{b}}_{I1}\left({\mathit{x}}_C\right)& 0& 0\\ 0& {\tilde{b}}_{I2}\left({\mathit{x}}_C\right)& 0\\ 0& 0& {\tilde{b}}_{I3}\left({\mathit{x}}_C\right)\\ {\tilde{b}}_{I2}\left({\mathit{x}}_C\right)& {\tilde{b}}_{I1}\left({\mathit{x}}_C\right)& 0\\ 0& {\tilde{b}}_{I2}\left({\mathit{x}}_C\right)& {\tilde{b}}_{I3}\left({\mathit{x}}_C\right)\\ {\tilde{b}}_{I1}\left({\mathit{x}}_C\right)& 0& {\tilde{b}}_{I3}\left({\mathit{x}}_C\right)\end{array}\right]$$

(23)

where

$${\tilde{b}}_{I1}\left({\mathit{x}}_C\right)=\frac{1}{V_C}{\int }_{{\Gamma }_C}{N}_I\left(\mathit{x}\right)n_k\left(\mathit{x}\right)d\mathbf{\Gamma }$$

(24)

$N_I$ is the regular shape functions for a eight-node hexahedral finite element. Since we are using eight-noded hexahedral elements, one Gauss point at the center is sufficient to integrate along each surface boundary. The above equation then reduces to

$${\tilde{b}}_{I1}\left({\mathit{x}}_C\right)=\frac{1}{V_C}\sum _{I=1}^M{N}_I\left({\mathit{x}}_i^{GP}\right)n_{ki}^C{A}_i^C$$

(25)

where $x_i$ is the midpoint or the intersection of the diagonals of the boundary segment and ${\mathit{n}}_{\mathit{i}}$ and $A_i$ are the corresponding outward unit normal and surface area, respectively.

The smoothed stiffness matrix of an element is thus given by the assembly of the stiffness of each smoothing cell

$${\tilde{K}}_e=\sum _C{\tilde{B}}_C^{T}C{\tilde{B}}_C{V}_C$$

(26)

The overall smoothed stiffness matrix of the domain is calculated by the summation of sub-stiffness matrices for nodes I in relation to node J as in FEM, except that the summation here is performed over the smoothing domains, not elements:

$${\tilde{K}}_{IJ}^{CS-FEM}=\sum _{i=1}^{N_e}\sum _{m=1}^{N_c}{\int }_{{\Omega }_{i,m}^s}{\tilde{B}}_I^{T}C{\tilde{B}}_{J}d\mathbf{\Omega }=\sum _{i=1}^{N_e}\sum _{m=1}^{N_c}{\tilde{B}}_I^{T}C{\tilde{B}}_J{V}_{i,m}^s=\sum _{k=1}^{N_e}{\tilde{B}}_I^{T}C{\tilde{B}}_J{V}_k^s$$

(27)

Thus, it is an algebraic summation over the stiffness of each element and the equations required to calculate the stiffness matrix of an element are the Equations (23) and (24) where we do not need the spatial derivative of the shape functions, thereby eliminating the Jacobian inversion problem. Only the values of shape function at certain points (Gauss point/center of each surface of a smoothing cell) are needed (Figure 2).

The calculation of displacements or force vector, based on the nature of the problem is similar to classical FEM, by solving the following algebraic equation

$${\tilde{\mathit{K}}}^{CS-FEM}\overline{\mathit{d}}=\overline{\mathit{f}}$$

(28)

Once the displacements are calculated, we deviate from the standard stress/strain calculations used in [1,2,38]. This is because ABAQUS does not let users obtain information from surrounding elements; the user element (UEL) is written such that it only specifies the formulation of a single element. In a typical CS-FEM setting [1,2], the stress/strains at a node is calculated by the weighted average of the values of the surrounding smoothing domains of the node. In this case, due to our inability to access data from surrounding elements and because of the way the UEL is formulated, we calculate the elemental strains as

$$\mathit{ϵ}=\tilde{B}\mathit{U}$$

(29)

where $\tilde{B}$ is the smoothed strain displacement matrix obtained at the center of the smoothing cell. The corresponding elemental stress is given by

$${\mathit{\sigma }}_{\mathit{ij}}={\mathit{C}}_{\mathit{ijkl}}{\mathit{ϵ}}_{\mathit{kl}}$$

(30)

ABAQUS then allows its internal algorithms to calculate the nodal stress/strain values. Similarly, the elemental displacement values are given as

$${\mathit{u}}_{\mathit{i}}=N_{ij}{\mathit{U}}_{\mathit{j}}$$

(31)

where $U_j$ is the $24\times 1$ nodal displacement vector and $u_i$ is the $3\times 1$ elemental displacement vector.

2.4. Implementation in ABAQUS UEL

We implemented the PhaseField-CSFEM in ABAQUS (ABAQUS 2017, Dassault Systemes, Providence, RI, USA) and solved it using a staggered scheme. Equations (13)–(17) contain the coupled phase and displacement fields which can be solved using a monolithic approach [26,28], however we decoupled and solved them separately using the SFEM in ABAQUS UEL. In the literature, the staggered scheme has presented quite a few stability advantages over the monolithic solution [25,27].

The decoupled governing equations are:

$$r_c={\int }_{\Omega }([\frac{g_c}{l_c}c-2(1-c)H]{\left(N^c\right)}^T+g_c{l}_c{\left(B^c\right)}^T\nabla c)\mathit{d}\mathbf{\Omega }$$

(32)

where $r_c$ is the residual force vector and $N^c$ and $B^c$ are the shape function and the strain gradient matrix for 3D eight-node hexahedral elements, respectively.

The phase-field component in the decoupled stiffness matrix is given by:

$$K_c={\int }_{\Omega }([\frac{G_c}{l_c}+2H]{\left(N^c\right)}^T\left(N^c\right)+G_c{l}_c{\left(B^c\right)}^T\left(B^c\right))\mathit{d}\mathbf{\Omega }$$

(33)

Similarly, the displacement-field contribution in the stiffness matrix is given by $B^u$ matrix calculated from Equation (23),

$$K_u={\int }_{\Omega }[{(1-c)}^2+d]{\left(B^u\right)}^T\mathit{C}\left(B^u\right)d\mathbf{\Omega }$$

(34)

and the internal force(residual) vector is

$$r_u={\int }_{\Omega }[{(1-c)}^2+d]{\left(B^u\right)}^T{\mathit{\sigma }}_{\mathbf{0}}d\mathbf{\Omega }$$

(35)

where ${\mathit{\sigma }}_{\mathbf{0}}$ is the stress calculated by Equation (30).

Thus, the following equation is solved using the modified Newton–Raphson scheme

$${\left[\begin{array}{c}u\\ c\end{array}\right]}_{t+\delta t}={\left[\begin{array}{c}u\\ c\end{array}\right]}_t-{\left[\begin{array}{cc}{K}_u& 0\\ 0& K_c\end{array}\right]}_t{\left[\begin{array}{c}{r}_u\\ r_c\end{array}\right]}_t$$

(36)

The term H in Equations (32) and (33) is a history field which is equal to the strain energy from the previous step. In the staggered scheme, in the first iteration of every load step, the history variable is updated via $H_{n+1}={\psi }_n$ and the displacement field (Equations (34) and (35)) is solved by updating the value of the phase-field from the previous load step ($c_n$). The history variable satisfies certain properties [25] and ensures that no penalty term is necessary to enforce the irreversibility of the crack field.

In this formulation, the $B^u$ matrix in Equations (34) and (35) is calculated based on the smoothing cells, using CS-FEM-H8.

The non-linear system of equations (Equation (36)) is solved via an incremental iterative approach using a Newton–Raphson scheme. The continuum mechanics equations are solved by the more accurate CS-FEM [2] and the phase-field equations are solved by the standard FEM. In the UEL, the user needs to specify the tangent stiffness matrix $AMATRX$ and the internal force vector $RHS$, which the solver calls for each element. The properties are calculated at the center of the smoothing domain for each subcell. In the staggered scheme of solutions [27], we use the $common$ block in the UEL to facilitate the transfer of variables. Through the $common$ block, the UEL allows a user to write formulations of multiple elements in a single code. We define three elements, the first two being the phase-field and the displacement elements and the third a dummy element written to facilitate post-processing in the ABAQUS viewer. ABAQUS is unaware of the inherent shape functions used in the element formulation(UEL) and thus is unable to extrapolate the results on its own. We write a user material (UMAT), where the results stored as solution variables $STATEV$ are transferred from the UEL. This is possible because the element connectivity and the shape function of the elements are exactly the same as the C3D8 elements in ABAQUS library. Thus, we imply that there is a dummy mesh of H8 elements whose material properties, chosen such that there is no resistance to strain, are provided to the UMAT. The corresponding variables are stored as $SVARS$ in the UEL and transferred to the UMAT via the $common$ statement.

3. Numerical Modeling

We present several numerical examples in this section to validate the accuracy of the 3D-CS-FEM-phase field method for linear elastic brittle fracture and also to discuss some shortcomings of the implementation in ABAQUS. First, three benchmark examples were tested and subsequently we simulated comparatively complex models. All models were meshed using eight-node hexahedral elements and eight smoothing domains were used for CS-FEM formulations. Since the phase-field represents a diffused fracture representation, no special crack-tip element was necessary for capturing the crack path. The mesh size was significantly reduced to successfully resolve the length scale parameter near the expected crack propagation zone. As discussed by Borden et al. [34], the approximate element size should be less than half of the corresponding length scale parameter. Contour plots representing the phase-field variable which signifies the crack path are presented for all examples.

3.1. Single Edge Notched Tensile Sample

We simulated crack propagation in a rectangular bar, with a finite opening, under far-field tension. The bottom surface was fixed and a tensile pull was applied on the top surface of the rectangular block (Figure 3). The crack was located at an edge and the material parameters were: E = 210 kN/mm${}^2$, $\nu$ = 0.3, and $G_c$ = 5 × ${10}^{-4}$ kN/mm. The length scale parameter was $l_c$ = 0.08 mm. We used a uniform mesh of 45,000 hexahedral elements, which was refined near the expected crack propagation zone (edge side of 0.02 mm) to successfully resolve the length scale parameter for better reproduction of the ultimate strength of the sample [25,39]. Since this was a diffused fracture representation, we did not need to use any special crack-tip elements; the length scale parameter provided the transition from the damaged to intact zone. As expected, we saw the contour plot of crack evolution; it propagated along the straight path (Figure 4). We also studied the influence of the Griffith’s energy release rate parameter $G_c$ on the crack pattern. We tested the model for $G_c$ = 5 × ${10}^{-4}$, 4 × ${10}^{-4}$, 2 × ${10}^{-4}$, and 1 × ${10}^{-4}$ kN/mm, and observd that this parameter had no influence on the ultimate crack path. The resultant load–displacement curve is presented in Figure 5. It shows that the slope of the curve is same for small displacements for different energy release rates. However, as displacement increased, the force reduced before it reached the ultimate load, which itself increased with $G_c$. There was a sudden drop of the load after the crack initiation, which signified decrease in strength and ultimate failure. We also performed experiments with further reduction of mesh size near the expected crack path, however that did not affect crack pattern and the force–displacement curve, thereby proving that the damage pattern was convergent.

3.2. Single Edge Notched Shear Sample

We next considered a similar example, to simulate crack propagation in a rectangular bar, with a finite opening, under shear loading. The crack was located at an edge and the material and length scale parameters were: E = 210 kN/mm${}^2$, $\nu$ = 0.2, $l_c$ = 0.2 mm, and $G_c$ = 2.7 × ${10}^{-3}$ kN/mm. We used a hexahedral mesh of 49,000 elements, refined with elements of edge size about 0.03 mm, near the expected crack path (Figure 6). The simulation used displacement controlled boundary condition, with the bottom surface fixed and a horizontal displacement acting on the top surface. We noticed a hint of damage originating from the side, after the deformation due to the crack beginning to propagate. The observed crack path was similar to the results obtained in Zeng et al. [22] who used ES-FEM-T3 to simulate crack propagation. We performed a similar mesh size reduction study as in Section 3.1, where our findings corroborated that the predicted crack patterns were convergent.

3.3. Double Edge Notched Tensile Sample

A doubly notched symmetric rectangular bar, subjected to uniaxial tension was tested in this example (Figure 7). The material and length scale parameters were: E = 210 kN/mm${}^2$, $\nu$ = 0.3, $l_c$ = 0.1 mm, and $G_c$ = 2.7 × ${10}^{-3}$ kN/mm. The model had approximately 35,000 H8 elements with reduced element size near the expected crack path and the boundary conditions were the same as presented in Section 3.1. The crack paths obtained by this code are presented in Figure 8, and were in excellent agreement with Msekh et al. [40].

3.4. 3D Specimen with Notch and Three Openings

In this example, we tested a relatively difficult problem to solve with traditional discrete fracture modeling techniques. In coordination with [40], we had a geometry with three openings and an initial notch under tensile pull (Figure 9). Initially, an incremental displacement of $\Delta u$ = ${10}^{-3}$ mm was applied. The relative positions of the openings and the notch were altered to see the difference in fracture pattern. The material and length scale parameters were: E = 210 kN/mm${}^2$, $\nu$ = 0.2, $l_c$ = 0.02 mm, and $G_c$ = 5.0 × ${10}^{-3}$ kN/mm. Here, we observed two different crack paths: when the openings were aligned, in the same line, the fracture zone from the notch merged with the one originating from the nearest opening and the crack then propagated through the openings. This is basically a phenomenon of crack arrest which comes into existence due to the merging of damaged zones. Although the stress field near the notch was singular, the field near the openings also suffered from the effect of stress concentration [4], which further led to damage. The beauty of this method lies in the fact that it can also elegantly capture the damage initiation from the openings without any ad-hoc criterion. Our results conformed with the findings of Msekh [40] (Figure 10). However, when the openings are not aligned symmetrically, we saw a difference in crack path behavior. The usual crack arrest phenomenon occurred as in the previous case, but, unlike in [40] (Figure 11), where the crack propagates through the first two holes and merges with the third hole and then continues, we saw that the crack propagation continued along the same line, with a slight deflection towards the damaged zone of the third opening. Additionally, another crack propagation started from the damaged zone of the third opening. We believe this is just a by-product of using different displacement increments for the simulation. When the step size of our displacement increments were reduced, $\Delta u$ = 10${}^{-5}$ mm, we observed a different crack path, as shown, in Figure 12 which was much closer to the reference results. This indicates a drawback of the method, where the size of each displacement increment has to be suitably chosen to simulate the correct phenomenon.

3.5. 3D Bi-Material Notched Specimen

This example highlighted one of the advantages of the phase-field method, where we reproduced the phenomenon of crack branching. The geometry and material parameters were chosen similar to those in [27], where we had a notch in the the softer material region with specifications given in Figure 13. The surface on the right was fixed while an incremental displacement pull was applied on the left. The plate had a thickness of 5 mm. We observed that the damage initiation and crack propagation occurred as usual in any edge-notched specimen but, as soon as the crack reached the stiffer material, it branched. This is because the energy required to propagate through the stiff material is higher than the energy required to branch and continue along the material boundaries. The resultant crack path was in pretty good agreement with the available result in [27].

3.6. Crack Propagation in Thick Walled Cylinder

We performed another numerical example to demonstrate the advantages of using the phase-field method to simulate crack propagation. We modeled two parallel cracks at different elevations in a thick walled cylinder under uni-axial tension. The dimensions were: t = 0.1 mm, r = 1 mm, a = 0.2 mm, and h = 1 mm (Figure 14). The material and length scale parameters were: E = 210 kN/mm${}^2$, $\nu$ = 0.2, $l_c$ = 0.075 mm, and $G_c$ = 2.7 × ${10}^{-3}$ kN/mm. We observed the crack merging phenomenon in this case. The cracks initiated from the respective notches but, when the damage zones were influenced by each other, the two cracks merged (Figure 15). This phenomenon was simulated based on the physics of the system without any ad-hoc criterion set up for the crack propagation path, thereby restating that the phase-field technique can be used to simulate many such cases that would traditionally be a bit more complex using discrete treatment.

4. Discussion

The primary objective of this study was to predict realistic crack paths in three-dimensional solids. The CS-FEM-H8 element is implemented in the commercial platform of ABAQUS. The element was already proven to have a very good accuracy in displacement norm [38]; here, we combined it with the phase-field model for diffuse fracture to simulate crack propagation in 3D linear elastic brittle solids. ABAQUS user elements (UEL) has its inherent limitations in not letting the user access data from surrounding elements, thus we developed a much simpler technique to calculate stress/strain for smoothing domains instead of the traditional weighted average approach. We sacrificed the stress convergence rate due to this new formulation, however still obtained accurate results that predict complex crack topologies such as crack branching and curved crack paths. Several drawbacks of the method were observed during the simulations: the phase-field variable propagation is very much mesh dependent, thus one needs to have a very fine mesh (approximate characteristic element size of $(1/2)$ the length scale parameter, or less [24,25,28]) near the expected crack propagation zone to efficiently resolve the length-scale parameter. This is not an issue for geometries where experimentally obtained crack paths are already known, however, for complex geometries and unknown crack paths, we needed to have a very fine mesh throughout the domain to facilitate independent crack propagation which increases the computational cost. Once the length scale parameter has been successfully resolved, further reduction of the mesh size does not alter the crack path. We also observed that the selection of displacement increment step size affects the resultant crack path. The final fracture pattern is not dependent on the Griffith’s energy release rate parameter, however the ultimate strength of the sample decreases with decreasing $G_c$. Implementation in ABAQUS, as already discussed, has its own drawbacks in not letting the user implement more accurate S-FEM models such as ES-FEM/$\alpha$-FEM [18]. This formulation was developed based on the H8 element with eight smoothing cells, which was not generated automatically for complex geometries. Analysts probably need to resort to the use of tetrahedral (T4) elements in cases of complex topologies. This in turn affects the accuracy because the formulations of the T4 elements are the same for both CS-FEM and FEM [1]. However, this opens up an avenue for future research in this domain, where in-house codes can be used to solve the same problem using a much more accurate ES-FEM/FS-FEM and a base mesh of T3/T4 elements, which are very easily generated for complex geometries. Moreover, adaptive meshing schemes designed based on values of the phase-field parameter can be utilized to reduce the computational cost of having very fine mesh. In 2D cases, CS-FEM is slightly more efficient than FEM [39]. However, in this 3D setting, the computational cost is similar (quantitative difference of around 1%) for both FEM and CS-FEM formulations. This is due to the assembly process of the stiffness matrix for each smoothing domain which involves integration over each surface of the cell. However, the better accuracy in the displacement norm and the comparative “soft” property in stiffness matrix with respect to FEM gives CS-FEM an advantage.

Thus, it can be concluded that this method is another efficient and novel computational technique to model crack propagation in 3D structures.