Fracture Mechanics Online Class

Numerical Methods > LEFM

Simple Finite Element Simulations

Picture VI.1: Example of conforming mesh.

In order to model crack propagation, a simple method is a finite element simulation in which the crack is explicitly modeled with adequate boundary conditions (BCs). The mesh is therefore conforming with the crack lips. Crack propagation can then be evaluated and conducted with the following steps:

Despite its simplicity, this method is not used because:

Cohesive elements


Picture VI.2: Barenblatt's model.

The cohesive method is based on Barenblatt's model, which is an idealization of the brittle fracture mechanisms that involves the breaking of atomic bonding at crack tips (cleavage), see Picture VI.2. As long as they are not separated by a distance $\delta_t$, the crack lips are under attractive forces (see overview lecture).

In the cohesive zone lecture, it will be shown that \begin{equation}J = \int_{0}^{\delta_t} \sigma_y(\delta) d\delta, \end{equation} so the area below the $\sigma - \delta$ curve corresponds to $G_C$ if the crack grows straight ahead.

Picture VI.3: Cohesive law or $\sigma - \delta$ curve.

This model is defined trough a cohesive law, see Picture VI.3, which requires only 2 parameters:

the shape of the curves has no effect for brittle materials as long as it is monotonically decreasing.

Computational framework

Picture VI.4: Cohesive elements.

To integrate the cohesive law such as depicted in Picture VI.3, cohesive elements are inserted between two volume elements, see Picture VI.4. In a 3D context, the computation of the opening is calculated with the normal to the interface in the deformed configuration $\mathbf{N^-}$.

The opening, also called jump, between the two volume elements can be decomposed into two parts: the normal opening is equal to $ \delta_n = \max( [\![ \mathbf{u}]\!].\mathbf{N^-},0)$ while the sliding reads $\mathbf{\delta_s} = [\![ \mathbf{u}]\!] - [\![ \mathbf{u}]\!].\mathbf{N^-}\mathbf{N^-}$.

The resulting (or effective) opening is obtained as $ \delta = \sqrt{\delta_n^2 + \beta_c^2 \Vert\delta_s\Vert^2}$ with $\beta_c$ the ratio between the shear and opening critical stress.

A potential $\phi = \phi(\delta)$ can be used to derive the 3D cohesive law, which is also called traction separation law (TSL). The traction (in the deformed configuration) across the crack lips derives from this potential: \begin{equation} \mathbf{t} = \dfrac{\partial \phi}{\partial \mathbf{\delta}} = \dfrac{\partial \phi}{\partial \delta_n} \mathbf{N^-} + \dfrac{\partial \phi}{\partial \delta_s} \dfrac{\partial \mathbf{\delta_s}}{\delta_s}. \end{equation}

Crack insertion

Two approaches can be used to insert the cohesive element, see Picture VI.4, in the finite element mesh.

The first method uses an intrinsic cohesive law:

Picture VI.5: Intrinsic cohesive law.

The pre-failure response and the failure criterion are incorporated within the cohesive law, see Picture VI.5. The cohesive elements can then be inserted since the very beginning of the simulation as they are modeling the pre-failure response. When the maximum traction is reached, the TSL has an irreversible part (during crack opening) and a reversible part (during crack closing).

The main advantage of this method is its ease of implementation: one has only to insert cohesive elements in the initial finite element mesh.

However the method suffers from serious drawbacks:

This intrinsic method can however be used with an a priori knowledge of the crack path. Indeed, in this case the cohesive elements can be inserted in a limited zone, reducing the artificial compliance issue. For example, the case of delamination in composite laminates can be treated in this way.

The second method uses an extrinsic law:

Picture VI.6: Extrinsic cohesive law.

The extrinsic cohesive law is modeling the failure response of the crack tip only, see Picture VI.6. In this case, the simulation begins with a usual continuous finite element mesh, and when the failure criterion is met ($\sigma \geq \sigma_{max}$) [Ortiz & Pandolfi 1999], interface elements are dynamically inserted in the finite element mesh to integrate the extrinsic TSL.

This method does not suffer from the artifical compliance of the intrinsic method. However, the implementation in 3D is complex especially for parallel simulations, which are mandatory in fracture mechanics.

Improvement: The hybrid Discontinuous Galerkin/Extrinsic Cohesive Law method

Picture VI.7: Fragmentation of an impacted ceramic plate [Radovitzky et al.].
Picture VI.8: Crack propagation in a blasted aluminum alloy cylinder [Becker et al.].
Picture VI.9: Crack propagation in a detonated canister [Becker et al.].
Picture VI.10: Composite fiber debonding [Wu et al.].

In order to recover the consistency when introducing the interface elements at the beginning of the simulation, the hybrid Discontinuous Galerkin/Extrinsic Cohesive Law method has recently been developed. Examples of simulations can be seen in Pictures VI.7, VI.8, VI.9, and VI.10.

Advantages and drawbacks of the cohesive zone method

Picture VI.11: Mesh dependency of the crack path.



The eXtended Finite Element Method

Key principles

The previously presented methods require the finite element mesh to be conforming with the crack. As a results the meshing operation can be complicated and requires fine discretizations. It would thus be advantageous to allow for discontinuities in the finite element model to be non conforming with the mesh as in Picture VI.12.

Picture VI.12: Discontinuity in a Finite Element (FE) discretization.

As a reminder, for a classical Finite Element (FE) discretization, the displacement field at coordinates $\xi^i$ is estimated by the sum on the nodes in the set $I$ (11 nodes here in Picture VI.12) of the product of the shape functions with the nodal displacements, i.e.:

\begin{equation} \mathbf{u}_h \left( \xi^i \right) = \sum_{a \in I} N^a \left( \xi^i \right) \mathbf{u}^a, \label{eq:FEinterpolation}\end{equation}


However, such interpolation (\ref{eq:FEinterpolation}) of the displacement field cannot account for a discontinuiy. New degrees of freedom are thus required in the model.

Picture VI.13: New nodes could be added to account for the discontinuity

The new degrees of freedom could be introduced by inserting new nodes near the crack tips (red crosses in Picture VI.13), but this would be inefficient as it corresponds to remeshing. Instead, in the eXtended Finite Element Method (XFEM), the new degrees of freedom are introduced through new shape functions in the displacement field interpolation (\ref{eq:FEinterpolation}). As the discontinuity does not affect all the finite elements, only the shape functions related to the nodes linked to the edges that intersect the crack are modified (blue dots in Picture VI.13). The displacement field (\ref{eq:FEinterpolation}) thus becomes:

\begin{equation} \mathbf{u}_h \left( \xi^i \right) = \sum_{a \in I} N^a \left( \xi^i \right) \mathbf{u}^a + \sum_{a \in J} N^a \left( \xi^i \right) F^a \mathbf{u}^{*a},\label{eq:FEModifiedInterpolation} \end{equation}


Picture VI.14: The signed distance of a domain from the crack surface is needed.

The new shape functions $F^a$ have now to be defined so that they can represent the discontinuity in the domain. A logical choice is the Heaviside's function defined such that their evaluation is $+1$ above and $-1$ below the crack. In order to know if we are above or below the crack, the signed-distance from the crack surface should be available for any domain point. This is achieved by defining the level set functions of the domain related to the crack. The displacement field (\ref{eq:FEModifiedInterpolation}) then becomes:

\begin{equation} \mathbf{u}_h \left( \xi^i \right) = \sum_{a \in I} N^a \left( \xi^i \right) \mathbf{u}^a + \sum_{a \in J} N^a \left( \xi^i \right) H(lsn(\xi^i,\xi^{*i})) \mathbf{u}^{*a},\label{eq:FEModifiedInterpolation2} \end{equation}


Picture VI.15: Fracture mechanics solution is introduced through the enrichment of the elements containing the crack tip.

At this point, a discontinuity can be introduced in the mesh and problems with a crack can be solved without the need for a conforming mesh. However, in the context of LEFM one can take advantage of the knowledge of the asymptotic solution expression to improve the accuracy in more efficient way than with a mesh refinement. Toward this end, a new enrichment is needed on the nodes of the elements embedding the crack tip (red dots in Picture VI.15). Indeed, as the asymptotic LEFM solution vanish far away from the crack tip, only closest nodes to the crack tip have to be enriched. The displacement field (\ref{eq:FEModifiedInterpolation2}) thus becomes:

\begin{equation}\begin{cases} \mathbf{u}_h \left( \xi^i \right) &=& \sum_{a \in I} N^a \left( \xi^i \right) \mathbf{u}^a + \sum_{a \in J} N^a \left( \xi^i \right) H(lsn(\xi^i,\xi^{*i})) \mathbf{u}^{*a} \\ &&+\sum_{a \in K} N^a \left( \xi^i \right) \sum_{b} \Psi_b(\xi^i) \mathbf{\psi}_b^a,\end{cases}\label{eq:FEModifiedInterpolation3} \end{equation}


Picture VI.16: Referential $x'$, $y'$, $z'$ linked to the crack tip direction.

What remain to be defined are the new shape functions $\Psi_b$. Toward this end, the asymptotic displacement field is expressed, in the referential $x'$, $y'$, $z'$ linked to the crack tip, see Picture VI.16, in terms of the three stress intensity factors:

\begin{equation} \begin{cases} \mathbf{u}_{x'} & = \frac{1+\nu}{E} \sqrt{\frac{r}{2\pi}} &\left\lbrace K_I \cos{\frac{ \theta}{2}} \left[\kappa-1+2 \sin^2{\frac{\theta}{2}}\right] +\right.\\ &&\left.K_{II} \sin{\frac{ \theta}{2}} \left[\kappa+1+2 \cos^2{\frac{\theta}{2}}\right] \right\rbrace, \\ \mathbf{u}_{y'} & = \frac{1+\nu}{E} \sqrt{\frac{r}{2\pi}} &\left\lbrace K_I \sin{\frac{ \theta}{2}} \left[\kappa+1-2 \cos^2{\frac{\theta}{2}}\right] + \right.\\ &&\left. K_{II} \cos{\frac{ \theta}{2}} \left[1-\kappa+2 \sin^2{\frac{\theta}{2}}\right] \right\rbrace, \\ \mathbf{u}_{z'} & = \frac{1+\nu}{E} \sqrt{\frac{r}{2\pi}} 2& K_{III} \sin{\frac{\theta}{2}}. \end{cases} \end{equation}

Basic trigonometric relations allow these fields to be rewritten

\begin{equation} \begin{cases} \mathbf{u}_{x'} &=& \frac{1+\nu}{E} \sqrt{\frac{r}{2\pi}} \left\lbrace K_I \cos{\frac{\theta}{2}} \left[ \kappa-1 \right] + K_I \sin{\frac{\theta}{2}} \sin{\theta} + K_{II} \sin{\frac{\theta}{2}} \left[ \kappa+1 \right] + K_{II} \cos{\frac{\theta}{2}} \sin{\theta} \right\rbrace,\\ \mathbf{u}_{y'} &=& \frac{1+\nu}{E} \sqrt{\frac{r}{2\pi}} \left\lbrace K_I \sin{\frac{\theta}{2}} \left[\kappa+1\right] - K_I \cos{\frac{\theta}{2}} \sin{\theta} + K_{II} \cos{\frac{\theta}{2}} \left[1-\kappa\right] + K_{II} \sin{\frac{\theta}{2}} \sin{\theta} \right\rbrace, \\ \mathbf{u}_{z'} & = & \frac{1+\nu}{E} \sqrt{\frac{r}{2\pi}} 2 K_{III} \sin{\frac{\theta}{2}}.\end{cases}\label{eq:manipulation1}\end{equation}

We have then defined four shape functions $\Psi_b$. In the global axes, the shape functions $\Psi_b$ remain the same as $\theta$ is linked to the crack tip direction and we can write

\begin{equation}\begin{cases} \mathbf{u}_x &=& \psi_x^1 \underbrace{ \sqrt{r} \sin{\frac{\theta}{2}} }_{\Psi_1} + \psi_x^2 \underbrace{ \sqrt{r} \cos{\frac{\theta}{2}} }_{\Psi_2} + \psi_x^3 \underbrace{ \sqrt{r} \sin{\frac{\theta}{2}} \sin{\theta}}_{\Psi_3} + \psi_x^4 \underbrace{ \sqrt{r} \cos{\frac{\theta}{2}} \sin{\theta} }_{\Psi_4},\\ \mathbf{u}_y &=& \psi_y ^1 \underbrace{ \sqrt{r} \sin{\frac{\theta}{2}} }_{\Psi_1} + \psi_y^2 \underbrace{ \sqrt{r} \cos{\frac{\theta}{2}} }_{\Psi_2} + \psi_y^3 \underbrace{ \sqrt{r} \sin{\frac{\theta}{2}} \sin{\theta}}_{\Psi_3} + \psi_y^4 \underbrace{ \sqrt{r} \cos{\frac{\theta}{2}} \sin{\theta} }_{\Psi_4}\,,\\ \mathbf{u}_z &=& \psi_z ^1 \underbrace{ \sqrt{r} \sin{\frac{\theta}{2}} }_{\Psi_1} + \psi_z^2 \underbrace{ \sqrt{r} \cos{\frac{\theta}{2}} }_{\Psi_2} + \psi_z^3 \underbrace{ \sqrt{r} \sin{\frac{\theta}{2}} \sin{\theta}}_{\Psi_3} + \psi_z^4 \underbrace{ \sqrt{r} \cos{\frac{\theta}{2}} \sin{\theta} }_{\Psi_4}. \end{cases}\label{eq:manipulation2} \end{equation}

We have thus 12 new degrees of freedom $\psi_b^a$ on the LEFM-enriched nodes. Note that as $\Psi_1$ is discontinuous, the Heaviside functions are not required anymore on the nodes in the subset $K$. The enriched displacement field (\ref{eq:FEModifiedInterpolation3}) finally becomes

\begin{equation} \begin{cases}\mathbf{u}_h \left( \xi^i \right) &=& \sum_{a \in I} N^a \left( \xi^i \right) \mathbf{u}^a + \sum_{a \in J} N^a \left( \xi^i \right) H(lsn(\xi^i,\xi^{*i})) \mathbf{u}^{*a} +\\ && \sum_{a \in K} N^a \left( \xi^i \right) \sum_{b=1}^4 \Psi_b\left(r(\xi^i), \theta(\xi^i)\right) \mathbf{\psi}_b^a.\end{cases}\label{eq:xFEMInterpolation} \end{equation}
Picture VI.17: The normal and tangent level set functions.

In order to evaluate the new shape functions $\Psi_b\left(r(\xi^i), \theta(\xi^i) \right)$, both the distance $r$ from the crack tip and the direction $\theta$ from the crack tip tangent have to be evaluated, see Picture VI.17. Toward this end, the normal and tangential signed distances from the crack tip have to be obtained through the level set functions. The values of $r$ and $\theta$ are then evaluated by:

\begin{equation} \begin{cases} r &=& \sqrt{lsn(\xi^i,xi^{**})^2 + lst(\xi^i,xi^{**})^2}, \\ \theta &=& \pm \arctan{ \frac{lsn(\xi^i,x^{**})} {lst(\xi^i,x^{**})}}, \end{cases} \end{equation}


With this framework, the displacement and stress fields of a domain embedding a crack can be computed using non-conforming mesh discretizations. The next step is now to evaluate whether the crack will propagate or not, and if so in which direction.

SIFs computation

In the LEFM context, the crack propagation can be assessed from the computation of the SIFs. Toward this end, the SIFs could be deduced from the new degrees of freedom $\psi_b^a$. Indeed, comparing the displacement fields (\ref{eq:manipulation1}), e.g. in the $x'$-direction linked to the crack tip,

\begin{equation} \mathbf{u}_{x'} = \frac{1+\nu}{E} \sqrt{\frac{r}{2\pi}} \left\lbrace K_I \cos{\frac{\theta}{2}} \left[ \kappa-1 \right] + K_I \sin{\frac{\theta}{2}} \sin{\theta} + K_{II} \sin{\frac{\theta}{2}} \left[ \kappa+1 \right] + K_{II} \cos{\frac{\theta}{2}} \sin{\theta} \right\rbrace, \end{equation}

with the formulation (\ref{eq:manipulation2}), e.g. in the $x'$-direction linked to the crack tip

\begin{equation} \mathbf{u}_{x'} = \psi_{x'}^1 \underbrace{ \sqrt{r} \sin{\frac{\theta}{2}} }_{\Psi_1} + \psi_{x'}^2 \underbrace{ \sqrt{r} \cos{\frac{\theta}{2}} }_{\Psi_2} + \psi_{x'}^3 \underbrace{ \sqrt{r} \sin{\frac{\theta}{2}} \sin{\theta}}_{\Psi_3} + \psi_{x'}^4 \underbrace{ \sqrt{r} \cos{\frac{\theta}{2}} \sin{\theta} }_{\Psi_4}, \end{equation}

shows the relationship between the degrees of freedom $\psi_{b'}^a$ expressed in the referential linked to the crack tip and the SIFs $K_i$:

\begin{equation} \begin{cases} K_I&=&\frac{\psi_{x'}^3\sqrt{2\pi}E}{1+\nu},\\K_{II}&=&\frac{\psi_{x'}^4\sqrt{2\pi}E}{1+\nu},\\K_{III}&=&\frac{\psi_{z'}^1\sqrt{2\pi}E}{1+\nu}.\end{cases}\end{equation}

However, this way of proceeding is not accurate as it is obtained from the displacement field interpolation, whose accuracy strongly depends on the mesh discretization. Therefore the SIFs are usually deduced from the $J$-integral. Toward this end, we apply the methodology described previously, which is summarized by

Crack propagation criterion

With the developed tools, the xFEM method is an efficient tool to simulate crack propagation. This is achieved as follows


Numerical examples

Resection of a brain tumor [L. Vigneron et al.]

Picture VI.18: Resection of a brain tumor [L. Vigneron et al.]

As current neuronavigation systems cannot adapt to changing intraoperative conditions over time an experimental end-to-end system capable of updating 3D preoperative images in the presence of brain shift and successive resections was developed using the XFEM technology. The biomechanical model is deformed using the xFEM in order to model tumor resection without recourse to a conforming mesh, see Picture VI.18.

Advantages and drawbacks