Phase Field Fracture

Formulation

Minimization of Energy Functional

In the phase field framework, fracture processes can be determined via the minimization of the following energy functional over displacement field $\mathit{u}$ and phase field variable $d$:

$$\begin{array}{r}\mathrm{\Pi }(\mathit{u},d)={\int }_{\mathrm{\Omega }}\psi (\mathit{\epsilon }(\mathit{u}),d)\text{ d}\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}{g}_c\gamma (d,\mathrm{\nabla }d)\text{ d}\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}\mathit{b}\cdot \mathit{u}\text{ d}\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_N}\mathit{t}\cdot \mathit{u}\text{ d}\mathrm{\Gamma },\end{array}$$

where $g_c$ is the Griffith-type critical energy release rate and $\gamma$ is the regularized crack surface density function (per volume):

$$\begin{array}{r}\gamma (d,\mathrm{\nabla }d)=\frac{1}{2l}{d}^2+\frac{l}{2}|\mathrm{\nabla }d{|}^2,\end{array}$$

with $l$ being the length-scale parameter. The bulk elastic energy $\psi$ is assumed to take the following form:

$$\begin{array}{r}\psi (\mathit{\epsilon },d)=g(d){\psi }_{+}(\mathit{\epsilon })+{\psi }_{-}(\mathit{\epsilon }),\end{array}$$

with $g(d)=(1-d{)}^2$ being the degradation function that models partial loss of stiffness due to the presence of cracks.

Miehe's model [1] of spectral decomposition of the strain tensor is adopted:

$$\begin{array}{rl}{\psi }_{\pm }(\mathit{\epsilon })& =\frac{\lambda }{2}⟨\text{tr}(\mathit{\epsilon }){⟩}_{\pm }^2+\mu {\mathit{\epsilon }}_{\pm }:{\mathit{\epsilon }}_{\pm },\end{array}$$

where ${\mathit{\epsilon }}_{\pm }:={\mathrm{\Sigma }}_{a=1}^n⟨{ϵ}_a{⟩}_{\pm }{\mathit{n}}_a\otimes {\mathit{n}}_a$, with $\{{ϵ}_a{\}}_{a=1}^n$ and $\{{\mathit{n}}_a{\}}_{a=1}^n$ being the principal strains and principal directions of $\mathit{\epsilon }$, respectively. We have $⟨x{⟩}_{\pm }:=\frac{1}{2}(x\pm |x|)$ being the bracket operator. Miehe's model prevents crack generation in compression.

Strong Form

The governing equations and boundary conditions obtained by minimizing the total energy functional are:

$$\begin{array}{rlrl}\mathrm{\nabla }\cdot \mathit{\sigma }+\mathit{b}& =\mathbf{0}& & \text{in}\mathrm{\Omega },\\ \frac{g_c}{l}(d-l^2\mathrm{\Delta }d)& =2(1-d)\mathcal{H}& & \text{in}\mathrm{\Omega },\\ \mathit{u}& ={\mathit{u}}_D& & \text{on}{\mathrm{\Gamma }}_D,\\ \mathit{\sigma }\cdot \mathit{n}& =\mathit{t}& & \text{on}{\mathrm{\Gamma }}_N,\\ \mathrm{\nabla }d\cdot \mathit{n}& =0& & \text{on}\mathrm{\Gamma },\end{array}$$

where $\mathit{\sigma }=\frac{\partial \psi }{\partial \mathit{\epsilon }}=g(d){\mathit{\sigma }}^{+}+{\mathit{\sigma }}^{-}$, and we have

$$\begin{array}{r}{\mathit{\sigma }}^{\pm }=\frac{\partial {\psi }^{\pm }}{\partial \mathit{\epsilon }}=\lambda ⟨\text{tr}(\mathit{\epsilon }){⟩}_{\pm }\mathit{I}+2\mu {\mathit{\epsilon }}_{\pm }.\end{array}$$

The history variable prevents crack from healing itself during unloading:

$$\begin{array}{r}\mathcal{H}(\mathit{x},t)=\underset{s\in [0,t]}{max}{\psi }_{+}(\mathit{\epsilon }(\mathit{x},s)).\end{array}$$

Weak Form

For arbitrary test function $(\delta \mathit{u},\delta d)$, we pose the following variational problem

$$\begin{array}{rl}a((\mathit{u},d),(\delta \mathit{u},\delta d))& =F((\delta \mathit{u},\delta d)),\end{array}$$

where we have

$$\begin{array}{rl}a((\mathit{u},d),(\delta \mathit{u},\delta d))& ={\int }_{\mathrm{\Omega }}\mathit{\sigma }(\mathit{\epsilon },d):\mathrm{\nabla }\delta \mathit{u}\text{ d}\mathrm{\Omega }+g_c{\int }_{\mathrm{\Omega }}(\frac{d}{l}\delta d+l\mathrm{\nabla }d\cdot \mathrm{\nabla }\delta d)\text{ d}\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}2\mathcal{H}d\delta d\text{ d}\mathrm{\Omega },\\ F((\delta \mathit{u},\delta d))& ={\int }_{{\mathrm{\Gamma }}_N}\mathit{t}\cdot \delta \mathit{u}\text{ d}\mathrm{\Gamma }+{\int }_{\mathrm{\Omega }}\mathit{b}\cdot \delta \mathit{u}\text{ d}\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}2\mathcal{H}\delta d\text{ d}\mathrm{\Omega }.\end{array}$$

There are two common schemes to solve this coupled nonlinear problem: monolithic and staggered schemes. In this example, we use the staggered scheme. More details can be found in our paper [2].

Eigenvalue/Eigenvector Derivative with Degeneracy

Before we move to the implementation section, caveats on computing the derivative of eigenvalues and eigenvectors (especially with degenerate eigenvalues) are briefly discussed here. One may tend to fully rely on JAX automatic differentiation to compute the derivative of eigenvalues and eigenvectors. However, when repeated eigenvalues occur, JAX native jax.grad may fail and return np.nan, as discussed in this post. The issue has its own complexity, and is not resolved yet.

One workaround is to add a small random noise to the matrix so that it always has distinct eigenvalues. This approach proves to be effective in our implementation of the phase field method.

The second approach is to define custom derivative rules with knowledge to handle repeated eigenvalues. In our example, JAX-FEM needs to computes $\frac{\partial \mathit{\sigma }}{\partial \mathit{\epsilon }}$, which further requires to compute $\frac{\partial {\mathit{\epsilon }}^{+}}{\partial \mathit{\epsilon }}$ and $\frac{\partial {\mathit{\epsilon }}^{-}}{\partial \mathit{\epsilon }}$. More generally, if a second order tensor $\mathit{A}$ decompose as $\mathit{A}={\mathrm{\Sigma }}_{a=1}^n{\lambda }_a{\mathit{n}}_a\otimes {\mathit{n}}_a$ and we define tensor map $\mathit{F}(\mathit{A}):={\mathrm{\Sigma }}_{a=1}^{n}f({\lambda }_a){\mathit{n}}_a\otimes {\mathit{n}}_a$, then we are interested in computing $\frac{\partial \mathit{F}}{\partial \mathit{A}}$. The procedures are well presented in Miehe's paper [3], in particular, Eq. (19) is what we are concerned about. We implemented the algorithms in the file eigen.py. In this file, you will see how native AD of JAX fails on repeated eigenvalues, but once custom derivative rules are specified, the issues is resolved.

Finally, make sure your JAX version is up-to-date, since we have observed some possible unexpected behavior of the function np.linalg.eigh in older versions of JAX, e.g., 0.3.x version.

Execution

Run

python -m demos.phase_field_fracture.example

from the jax-fem/ directory.

Results

Results can be visualized with ParaWiew.

Deformation (x10)

Loading history and tensile force

References

[1] Miehe, Christian, Martina Hofacker, and Fabian Welschinger. "A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits." Computer Methods in Applied Mechanics and Engineering 199.45-48 (2010): 2765-2778.

[2] Xue, Tianju, Sigrid Adriaenssens, and Sheng Mao. "Mapped phase field method for brittle fracture." Computer Methods in Applied Mechanics and Engineering 385 (2021): 114046.

[3] Miehe, Christian, and Matthias Lambrecht. "Algorithms for computation of stresses and elasticity moduli in terms of Seth–Hill's family of generalized strain tensors." Communications in numerical methods in engineering 17.5 (2001): 337-353.