Phase field approximation of dynamic brittle fracture

Abstract

Numerical methods that are able to predict the failure of technical structures due to fracture are important in many engineering applications. One of these approaches, the so-called phase field method, represents cracks by means of an additional continuous field variable. This strategy avoids some of the main drawbacks of a sharp interface description of cracks. For example, it is not necessary to track or model crack faces explicitly, which allows a simple algorithmic treatment. The phase field model for brittle fracture presented in Kuhn and Müller (Eng Fract Mech 77(18):3625–3634, 2010) assumes quasi-static loading conditions. However dynamic effects have a great impact on the crack growth in many practical applications. Therefore this investigation presents an extension of the quasi-static phase field model for fracture from Kuhn and Müller (Eng Fract Mech 77(18):3625–3634, 2010) to the dynamic case. First of all Hamilton’s principle is applied to derive a coupled set of Euler-Lagrange equations that govern the mechanical behaviour of the body as well as the crack growth. Subsequently the model is implemented in a finite element scheme which allows to solve several test problems numerically. The numerical examples illustrate the capabilities of the developed approach to dynamic fracture in brittle materials.

Appendix: Finite element implementation

The set of coupled partial differential Eqs. (29) and (30) can be solved by standard multi-field finite element solvers. This section outlines the implementation of the present dynamic phase field model with the finite element code FEAP. The initial boundary value problem is solved simultaneously for the displacements $\varvec{u}$ and the crack field $s$. To start with, the weak forms of (29) and (30) are formulated

$$\begin{aligned}&\int \limits _{\varOmega }{\left( -\rho \ddot{\varvec{u}}\cdot \delta {\varvec{u}}-\varvec{\sigma }:\delta {{\varvec{\varepsilon }}}\right) }\ \mathrm{d}V\nonumber \\&\quad +{\int \limits _{\partial {\varOmega _t}}{\varvec{t}^*\cdot \delta \varvec{u}}\ \mathrm{d}A}=0 \end{aligned}$$

(80)

and

$$\begin{aligned}&\int \limits _{\varOmega }{\left( \left[ 2s\psi _{e,zs}-\mathcal{G}_c\frac{1-s}{2\epsilon }\right] \delta {s}\right. }\nonumber \\&\quad +{\left. 2\epsilon \mathcal{G}_c\nabla {s}\cdot \nabla {\delta {s}}\right) \mathrm{d}V}=0. \end{aligned}$$

(81)

Standard isoparametric, four node in 2D and eight node in 3D finite elements are used to discretize the body $\varOmega$. Therefore, the primary fields can be interpolated in the interior of an element by means of the nodal values $\underline{\varvec{u}}_I, {s}_I$ and the scalar bilinear or trilinear shape functions $N_I$

$$\begin{aligned} \begin{array}{ll} \underline{\varvec{u}}_h= \displaystyle {\sum _{I=1}^{N_e}}{N_I\underline{\varvec{u}}_I}&\text { and }{s}_h= \displaystyle {\sum _{I=1}^{N_e}}{N_I{{s}}_I} \end{array} \end{aligned}$$

(82)

where $N_e$ is the total number of nodes per element. The 2D, plane strain implementation is discussed in more detail in the following. However a 3D finite implementation is conceptually similar and straightforward. For further information about the finite element method the reader is referred to textbooks like Hughes [22]. In a two dimensional setting the spatial derivatives in (80) and (81) can be expressed by means of the matrices

$$\begin{aligned} \begin{array}{ll} \underline{\varvec{B}}_I^{\varvec{u}}= \left[ \begin{array}{ll} N_{I,1} &{} 0 \\ 0 &{} N_{I,2} \\ N_{I,2} &{} N_{I,1} \end{array} \right] \text { and } &{} \underline{\varvec{B}}_I^{s}= \left[ \begin{array}{l} N_{I,1} \\ N_{I,2} \end{array} \right] . \end{array} \end{aligned}$$

(83)

In these matrices $N_{I,1}$ and $N_{I,2}$ denote the differentiation of the shape functions with respect to the coordinates $x_1$ and $x_2$. This way, it is possible to express the gradients

$$\begin{aligned} \underline{{\varvec{\varepsilon }}}_h= \left[ \begin{array}{c} \varepsilon _{11_h}\\ \varepsilon _{22_h}\\ 2\varepsilon _{12_h} \end{array} \right] = \displaystyle {\sum _{I=1}^{N_e}}{\underline{\varvec{B}}_I^{\varvec{u}}\underline{\varvec{u}}_I} \end{aligned}$$

(84)

and

$$\begin{aligned} \underline{\nabla {s}}_h= \left[ \begin{array}{l} s_h,_1\\ s_h,_2 \end{array}\right] = \displaystyle {\sum _{I=1}^{N_e}}{\underline{\varvec{B}}_I^{s}{s}_I} \end{aligned}$$

(85)

in matrix notation. Furthermore, two matrices are introduced that represent the deviatoric part of $\underline{{\varvec{\varepsilon }}}_h$

$$\begin{aligned} \underline{{\varvec{\varepsilon }}_D}_h^*= \left[ \begin{array}{c} \frac{1}{2}(\varepsilon _{11_h}-\varepsilon _{22_h})\\ \frac{1}{2}(\varepsilon _{22_h}-\varepsilon _{11_h})\\ 2\varepsilon _{12_h} \end{array} \right] \end{aligned}$$

(86)

and

$$\begin{aligned} \underline{{\varvec{\varepsilon }}_{D}}_h= \left[ \begin{array}{c} \frac{1}{2}(\varepsilon _{11_h}-\varepsilon _{22_h})\\ \frac{1}{2}(\varepsilon _{22_h}-\varepsilon _{11_h})\\ \varepsilon _{12_h} \end{array} \right] . \end{aligned}$$

(87)

By means of the above terms, it is possible to express the elastic energy density as

$$\begin{aligned}&\psi _{\text {e}_h} = \frac{K_n}{2}{\langle \varepsilon _{11_h}+\varepsilon _{22_h}}\rangle _{-}^2\nonumber \\&\quad +({s}_h^2+\eta )\underbrace{\left[ \frac{K_n}{2}{\langle \varepsilon _{11_h} +\varepsilon _{22_h}}\rangle _{+}^2+\mu \left( {\underline{{\varvec{\varepsilon }}_D}_h^*}^T \underline{{\varvec{\varepsilon }}_{D}}_h\right) \right] }_{\displaystyle \psi _{\text {e,ts}_{h}}}, \end{aligned}$$

(88)

where the bracket terms $\langle \cdot \rangle _{+/-}$ are treated according to (12) and (13). Consequently the stress tensor is

$$\begin{aligned} \underline{\varvec{\sigma }}_h=&\dfrac{\partial \psi _{e_h}}{\partial {\varvec{\varepsilon }}_h}\nonumber \\ =&K_n\text {tr}^{-}(\underline{{\varvec{\varepsilon }}}_h)\underline{\varvec{1}}\nonumber \\&+(s_h^2+\eta )\left[ K_n\text {tr}^{+}(\underline{{\varvec{\varepsilon }}}_h)\underline{\varvec{1}}+2\mu \underline{{\varvec{\varepsilon }}_{D}}_h\right] , \end{aligned}$$

(89)

where

$$\begin{aligned} \underline{\varvec{1}}&= \left[ \begin{array}{l@{\quad }l@{\quad }l} 1&1&0 \end{array} \right] ^T \end{aligned}$$

(90)

is the second-order identity tensor in matrix notation. By means of these expressions one can reformulate the governing equations (80) and (81) in discretized form

$$\begin{aligned} \displaystyle \bigcup _{e=1}^{n_e}\displaystyle \sum _{I=1}^{{N_e}} {\delta \underline{\varvec{u}}_I^T\underline{\varvec{R}}_{I,e}^{\varvec{u}}}={{0}},\end{aligned}$$

(91)

$$\begin{aligned} \displaystyle \bigcup _{e=1}^{n_e}\sum _{I=1}^{N_e}\delta {s}_I{R}_{I,e}^{s}={{0}}, \end{aligned}$$

(92)

where $n_e$ and $N_e$ are the total number of elements and the number of nodes per element, respectively. The symbol $\bigcup$ denotes the assembly operator. The contribution of element $e$ at node $I$ to the global residuum are

$$\begin{aligned} \begin{array}{ll} \underline{\varvec{R}}_{I,e}^{\varvec{u}}&{}=\underbrace{\int \limits _{\varOmega _e} {(-{\left[ \underline{\varvec{B}}_I^{\varvec{u}}\right] }^T\underline{\varvec{\sigma }}_h-\rho {N_I}\underline{\ddot{\varvec{u}}}_h)\mathrm{d}V}}_{\displaystyle -\underline{\varvec{P}}_{I,e}^{\varvec{u}}}\\ &{}\quad +\underbrace{\int \limits _{\partial {\varOmega }_e^t\cap \partial {\varOmega ^t}} {N_I\underline{\varvec{t}}_h^*\mathrm{d}A}}_{\displaystyle \underline{\varvec{F}}_{I,e}^{\varvec{u}}} \end{array} \end{aligned}$$

(93)

and

$$\begin{aligned} \begin{array}{ll} {R}_{I,e}^{s}=&{}-{P}_{I,e}^s=\displaystyle \int \limits _{\varOmega _e}\left( -[\underline{\varvec{B}}_I^{s}]^T2\epsilon \mathcal{G}_c\underline{\nabla {s}}_h\right. \\ &{}\left. -N_I\left( 2{s}_h{\psi _{e,ts}}_h -\frac{\mathcal{G}_c}{2\epsilon }(1-{s}_h)\right) \right) \mathrm{d}V. \end{array} \end{aligned}$$

(94)

Terms that can be attributed to inner forces are denoted $\underline{\varvec{P}}_{I,e}^{\varvec{u}}$ and ${P}_{I,e}^{s}$ respectively, whereas external mechanical loads are labelled $\underline{\varvec{F}}_{I,e}^{\varvec{u}}$. The residuals $\underline{\varvec{R}}_{I,e}^{\varvec{u}}$ and ${R}_{I,e}^{s}$ have to be computed on element level and are subsequently assembled to a global residuum $\underline{\varvec{R}}$. The discretized equations (91) and (92) define a global system of equations

$$\begin{aligned} \delta \underline{\varvec{d}}^T\underline{\varvec{R}}=0\quad \Leftrightarrow \quad \underline{\varvec{R}}=\underline{\varvec{0}} \end{aligned}$$

(95)

where

$$\begin{aligned} \underline{\varvec{R}} = \bigcup _{e=1}^{n_e}{\underline{\varvec{R}}_{e}}= \left[ \begin{array}{lll} \left( \underline{\varvec{R}}_{1}\right) ^T&...&\left( \underline{\varvec{R}}_{N}\right) ^T \end{array} \right] ^T, \end{aligned}$$

(96)

and

$$\begin{aligned} \begin{array}{ll} \delta \underline{\varvec{d}}= \left[ \begin{array}{lll} \left( \delta \underline{\varvec{d}}_1\right) ^T &{}...&{}\left( \delta \underline{\varvec{d}}_N\right) ^T \end{array} \right] ^T,&{} \\ \delta \underline{\varvec{d}}_I= \left[ \begin{array}{ll} \left( \delta \underline{\varvec{u}}_I\right) ^T &{}\delta {s}_I \end{array} \right] ^T \end{array} \end{aligned}$$

(97)

are the global vectors of nodal residuals and nodal values of the test functions respectively.

1.1 Appendix B: Time discretization and linearization

By means of the internal and external forces $\underline{\varvec{P}}$ and $\underline{\varvec{F}}$ the global system of Eq. (95) reads

$$\begin{aligned} \underline{\varvec{R}}=\underline{\varvec{F}}-\underline{\varvec{P}}\left( {{\underline{\varvec{d}}}}, \dot{{{\underline{\varvec{d}}}}},\ddot{{{\underline{\varvec{d}}}}}\right) =\underline{\varvec{0}} \end{aligned}$$

(98)

or in time discretized form ($t_{n+1}=t_n+\Delta {t}$)

$$\begin{aligned} \underline{\varvec{R}}_{n+1}=\underline{\varvec{F}}_{n+1}-\underline{\varvec{P}} \left( {{\underline{\varvec{d}}}}_{n+1},\dot{{{\underline{\varvec{d}}}}}_{n+1},\ddot{{{\underline{\varvec{d}}}}}_{n+1}\right) =\underline{\varvec{0}}. \end{aligned}$$

(99)

Herein, the subscript $n+1$ denotes terms that are evaluated at time $t_{n+1}$. In order to approximate $\dot{{{\underline{\varvec{d}}}}}_{n+1}$ and $\ddot{{{\underline{\varvec{d}}}}}_{n+1}$ the implicit Newmark method

$$\begin{aligned} \dot{{\underline{\varvec{d}}}}_{n+1}&= \dot{{\underline{\varvec{d}}}}_n + \Delta t\ \ddot{{\underline{\varvec{d}}}}_\gamma , \text { where }\nonumber \\ \ddot{{\underline{\varvec{d}}}}_\gamma&= (1-\gamma )\ \ddot{{\underline{\varvec{d}}}}_n+\gamma \ \ddot{{\underline{\varvec{d}}}}_{n+1}\,,\end{aligned}$$

(100)

$$\begin{aligned} {{\underline{\varvec{d}}}}_{n+1}&= {{\underline{\varvec{d}}}}_n+\Delta t\ \dot{{\underline{\varvec{d}}}}_n+\frac{1}{2}\Delta t^2\ \ddot{{\underline{\varvec{d}}}}_\beta , \text { where }\nonumber \\ \ddot{{\underline{\varvec{d}}}}_\beta&= (1-2\beta )\ \ddot{{\underline{\varvec{d}}}}_n+2\beta \ \ddot{{\underline{\varvec{d}}}}_{n+1}\, \end{aligned}$$

(101)

with parameters $\gamma =\beta =0.5$ is utilized. The time step size $\Delta t$ is chosen by an automatic time step size control that adapts the time step size according to the number of Newton iterations needed for convergence in the previous time step. Substituting the relations (100) and (101) in (99) yields a system of nonlinear equations for the unknowns $\varvec{\underline{d}}_{n+1}$ that has to be solved in every time step

$$\begin{aligned} \underline{\varvec{R}}_{n+1}&=\widehat{\underline{\varvec{R}}}\left( \underline{\varvec{F}}_{n+1}, {{\underline{\varvec{d}}}}_{n},\dot{{{\underline{\varvec{d}}}}}_{n},\ddot{{{\underline{\varvec{d}}}}}_{n},{{\underline{\varvec{d}}}}_{n+1}\right) \nonumber \\&=\underline{\varvec{F}}_{n+1}-\widehat{\underline{\varvec{P}}} \left( {{\underline{\varvec{d}}}}_{n},\dot{{{\underline{\varvec{d}}}}}_{n},\ddot{{{\underline{\varvec{d}}}}}_{n},{{\underline{\varvec{d}}}}_{n+1}\right) =\underline{\varvec{0}}. \end{aligned}$$

(102)

For this purpose, Newton’s method is applied

$$\begin{aligned} \underline{\varvec{R}}_{n+1}^{(k+1)}&=\underline{\varvec{R}}_{n+1}^{(k)}+\Delta {\underline{\varvec{R}}_{n+1}^{(k)}}\nonumber \\&\approx \underline{\varvec{R}}_{n+1}^{(k)}+\underbrace{D{\underline{\varvec{R}}_{n+1}^{(k)}}}_{=-\underline{\varvec{S}}_{n+1}^{(k)}} \Delta {{{\underline{\varvec{d}}}}_{n+1}}\nonumber \\&= \underline{\varvec{R}}_{n+1}^{(k)}-\underline{\varvec{S}}_{n+1}^{(k)}\Delta {{{\underline{\varvec{d}}}}_{n+1}^{(k)}}=\underline{\varvec{0}} \end{aligned}$$

(103)

where the index $k$ indicates the current Newton iteration and $\underline{\varvec{S}}_{n+1}^{(k)}$ is the overall tangent matrix of the problem. Equation (103) defines an iteration scheme, in which a system of linear equations (SLE) has to be solved for $\Delta {{{\underline{\varvec{d}}}}_{n+1}^{(k)}}$ in each iteration. Subsequently the vector of unknowns is updated according to

$$\begin{aligned} {{{\underline{\varvec{d}}}}_{n+1}^{(k+1)}} = {{{\underline{\varvec{d}}}}_{n+1}^{(k)}} +\Delta {{{\underline{\varvec{d}}}}_{n+1}^{(k)}} \end{aligned}$$

(104)

and the norm $\left\| \underline{\varvec{R}}\left( {{\underline{\varvec{d}}}}_{n+1}^{(k+1)}\right) \right\|$ is compared to a predefined tolerance $\delta$ to check for convergence. The overall tangent matrix is computed from

$$\begin{aligned} \underline{\varvec{S}}_{n+1}^{(k)}&= \frac{\partial \widehat{\underline{\varvec{P}}}}{\partial {{\underline{\varvec{d}}}}_{n+1}^{(k)}}\left( {{\underline{\varvec{d}}}}_{n},\dot{{{\underline{\varvec{d}}}}}_{n},\ddot{{{\underline{\varvec{d}}}}}_{n},{{\underline{\varvec{d}}}}_{n+1}^{(k)}\right) \nonumber \\&= \frac{\mathrm{d}\underline{\varvec{P}}}{\mathrm{d}{{\underline{\varvec{d}}}}_{n+1}^{(k)}}\left( {{\underline{\varvec{d}}}}_{n+1}^{(k)},\dot{{{\underline{\varvec{d}}}}}_{n+1}^{(k)},\ddot{{{\underline{\varvec{d}}}}}_{n+1}^{(k)}\right) \nonumber \\&= \underline{\varvec{K}}^{(k)}+\underline{\varvec{D}}^{(k)}\frac{\partial \dot{{{\underline{\varvec{d}}}}}_{n+1}^{(k)}}{\partial {{\underline{\varvec{d}}}}_{n+1}^{(k)}} +\underline{\varvec{M}}^{(k)}\frac{\partial \ddot{{{\underline{\varvec{d}}}}}_{n+1}^{(k)}}{\partial {{\underline{\varvec{d}}}}_{n+1}^{(k)}}\nonumber \\&= \underline{\varvec{K}}^{(k)}+\frac{\gamma }{\beta \Delta {t}}\underline{\varvec{D}}^{(k)}+\frac{1}{\beta \Delta {t}^2}\underline{\varvec{M}}^{(k)}. \end{aligned}$$

(105)

Similar to the residuals, the stiffness matrix

$$\begin{aligned} \underline{\varvec{K}}^{(k)}=\frac{\partial \underline{\varvec{P}}}{\partial {{\underline{\varvec{d}}}}_{n+1}^{(k)}}, \end{aligned}$$

(106)

the damping matrix

$$\begin{aligned} \underline{\varvec{D}}^{(k)}=\frac{\partial \underline{\varvec{P}}}{\partial \dot{{{\underline{\varvec{d}}}}}_{n+1}^{(k)}}=\varvec{\underline{0}}, \end{aligned}$$

(107)

and the mass matrix

$$\begin{aligned} \underline{\varvec{M}}^{(k)}=\frac{\partial \underline{\varvec{P}}}{\partial \ddot{{{\underline{\varvec{d}}}}}_{n+1}^{(k)}}. \end{aligned}$$

(108)

are computed for each finite element separately and are assembled to the global matrices. The stiffness matrix of each four node element can be written as

$$\begin{aligned} \underline{\varvec{K}}_e^{(k)}= \left[ \begin{array}{l@{\quad }l@{\quad }l@{\quad }l} \underline{\varvec{K}}_{11,e}^{(k)} &{}\underline{\varvec{K}}_{12,e}^{(k)} &{}\underline{\varvec{K}}_{13,e}^{(k)} &{}\underline{\varvec{K}}_{14,e}^{(k)}\\ \underline{\varvec{K}}_{21,e}^{(k)} &{}\underline{\varvec{K}}_{22,e}^{(k)} &{}\underline{\varvec{K}}_{23,e}^{(k)} &{}\underline{\varvec{K}}_{24,e}^{(k)}\\ \underline{\varvec{K}}_{31,e}^{(k)} &{}\underline{\varvec{K}}_{32,e}^{(k)} &{}\underline{\varvec{K}}_{33,e}^{(k)} &{}\underline{\varvec{K}}_{34,e}^{(k)}\\ \underline{\varvec{K}}_{41,e}^{(k)} &{}\underline{\varvec{K}}_{42,e}^{(k)} &{}\underline{\varvec{K}}_{43,e}^{(k)} &{}\underline{\varvec{K}}_{44,e}^{(k)} \end{array} \right] , \end{aligned}$$

(109)

where the submatrices $\underline{\varvec{K}}_{IJ,e}^{(k)}$ have the structure

$$\begin{aligned} \underline{\varvec{K}}_{IJ,e}^{(k)}= \left[ \begin{array}{ll} \underline{\varvec{K}}_{IJ,e}^{\varvec{u}\varvec{u}} &{}\underline{\varvec{K}}_{IJ,e}^{\varvec{u}{s}}\\ \underline{\varvec{K}}_{IJ,e}^{s\varvec{u}} &{} {K}_{IJ,e}^{ss} \end{array} \right] . \end{aligned}$$

(110)

Note that the iteration index is omitted for readability. The submatrices follow from

$$\begin{aligned} \underline{\varvec{K}}_{IJ,e}^{\varvec{u}\varvec{u}}&=\frac{\partial \underline{\varvec{P}}_{I,e}^{\varvec{u}}}{\partial \underline{\varvec{u}}_J}\nonumber \\&=\int \limits _{\varOmega _e}{\left[ \underline{\varvec{B}}_I^{\varvec{u}}\right] ^T\frac{\partial \underline{\varvec{\sigma }}_h}{\partial \underline{{\varvec{\varepsilon }}}_h}\left[ \frac{\partial \underline{{\varvec{\varepsilon }}}_h}{\partial \underline{\varvec{u}}_J}\right] \mathrm{d}V}\nonumber \\&=\int \limits _{\varOmega _e}{\left[ \underline{\varvec{B}}_I^{\varvec{u}}\right] ^T\frac{\partial \underline{\varvec{\sigma }}_h}{\partial \underline{{\varvec{\varepsilon }}}_h}\left[ \underline{\varvec{B}}_J^{\varvec{u}}\right] \mathrm{d}V},\end{aligned}$$

(111)

$$\begin{aligned} \underline{\varvec{K}}_{IJ,e}^{\varvec{u}{s}}&=\frac{\partial \underline{\varvec{P}}_{I,e}^{\varvec{u}}}{\partial {s}_J}\nonumber \\&=\int \limits _{\varOmega _e}{\left[ \underline{\varvec{B}}_I^{\varvec{u}}\right] ^T\frac{\partial \underline{\varvec{\sigma }}_h}{\partial {s_h}}\frac{\partial {s_h}}{\partial {s_J}}\mathrm{d}V}\nonumber \\&=\int \limits _{\varOmega _e}{\left[ \underline{\varvec{B}}_I^{\varvec{u}}\right] ^T\frac{\partial \underline{\varvec{\sigma }}_h}{\partial {s_h}}N_J^s\mathrm{d}V},\end{aligned}$$

(112)

$$\begin{aligned} \underline{\varvec{K}}_{IJ,e}^{{s}\varvec{u}}&=\frac{\partial {{P}}_{I,e}^{s}}{\partial {\underline{\varvec{u}}}_J}\nonumber \\&=\int \limits _{\varOmega _e}{N_I^{{s}}2\underline{s}_h\left[ \frac{\partial {\psi _{\text {e,ts}}}_h}{\partial {\underline{{\varvec{\varepsilon }}}_h}}\right] ^T\left[ \frac{\partial \underline{{\varvec{\varepsilon }}}_h}{\partial \underline{\varvec{u}}_J}\right] \mathrm{d}V}\nonumber \\&=\int \limits _{\varOmega _e}{N_I^{{s}}2\underline{s}_h\left[ \frac{\partial {\psi _{\text {e,ts}}}_h}{\partial {\underline{{\varvec{\varepsilon }}}_h}}\right] ^T\left[ \underline{\varvec{B}}_J^{\varvec{u}}\right] \mathrm{d}V},\end{aligned}$$

(113)

$$\begin{aligned} {K}_{IJ,e}^{{ss}} =&\,\,\frac{\partial {{P}}_{I,e}^{s}}{\partial {s}_J}\nonumber \\ =&\int \limits _{\varOmega _e}\left( \left[ \underline{\varvec{B}}_I^{s}\right] ^T\left[ 2\mathcal{G}_c\epsilon \right] \frac{\partial \underline{\nabla {s}}_h}{\partial {s_J}}\right. \nonumber \\&+\left. N_I\left( 2{\psi _{\text {e,ts}}}_h+\frac{\mathcal{G}_c}{2\epsilon }\right) \frac{\partial \underline{s}_h}{\partial {s_J}}\right) \mathrm{d}V\nonumber \\ =&\int \limits _{\varOmega _e}\left( \left[ \underline{\varvec{B}}_I^{s}\right] ^T\left[ 2\mathcal{G}_c\epsilon \right] \underline{\varvec{B}}_J^{s}\right. \nonumber \\&+\left. N_I\left( 2{\psi _{\text {e,ts}}}_h+\frac{\mathcal{G}_c}{2\epsilon }\right) N_J^s\right) \mathrm{d}V. \end{aligned}$$

(114)

The expressions

$$\begin{aligned} \frac{\partial \underline{\varvec{\sigma }}_h}{\partial \underline{{\varvec{\varepsilon }}}_h}={\varvec{\underline{\mathbb {C}}}}= \left[ \begin{array}{l@{\quad }l@{\quad }l} \mathbb {C}_{11} &{}\mathbb {C}_{12} &{}\mathbb {C}_{13}\\ \mathbb {C}_{21} &{}\mathbb {C}_{22} &{}\mathbb {C}_{23}\\ \mathbb {C}_{31} &{}\mathbb {C}_{32} &{}\mathbb {C}_{33} \end{array} \right] \end{aligned}$$

(115)

with

$$\begin{aligned} \mathbb {C}_{11}&=\mathbb {C}_{22}=\frac{\partial \text {tr}^{-}(\underline{{\varvec{\varepsilon }}}_h)}{\partial \text {tr}(\underline{{\varvec{\varepsilon }}}_h)}K_n\nonumber \\&\quad +(s_h^2+\eta )\left( \frac{\partial \text {tr}^{+}(\underline{{\varvec{\varepsilon }}}_h)}{\partial \text {tr}(\underline{{\varvec{\varepsilon }}}_h)}K_n+\mu \right) ,\nonumber \\ \mathbb {C}_{12}&=\mathbb {C}_{21}=\frac{\partial \text {tr}^{-}(\underline{{\varvec{\varepsilon }}}_h)}{\partial \text {tr}(\underline{{\varvec{\varepsilon }}}_h)}K_n\nonumber \\&\quad +(s_h^2+\eta )\left( \frac{\partial \text {tr}^{+}(\underline{{\varvec{\varepsilon }}}_h)}{\partial \text {tr}(\underline{{\varvec{\varepsilon }}}_h)}K_n-\mu \right) ,\nonumber \\ \mathbb {C}_{33}&=(s_h^2+\eta )\mu ,\nonumber \\ \mathbb {C}_{23}&=\mathbb {C}_{32}=\mathbb {C}_{13}=\mathbb {C}_{31}=0 \end{aligned}$$

(116)

is the material tangent for plane strain, whereas

$$\begin{aligned} \frac{\partial {\psi _{\text {e,ts}}}_h}{\partial {\underline{{\varvec{\varepsilon }}}_h}} =&K_n\text {tr}^{+}(\underline{{\varvec{\varepsilon }}}_h)\frac{\partial \text {tr}^{+}(\underline{{\varvec{\varepsilon }}}_h)}{\partial \underline{{\varvec{\varepsilon }}}_h}\nonumber \\&+\mu \frac{\partial \left( \underline{{\varvec{\varepsilon }}_D}_h^*\cdot \underline{{\varvec{\varepsilon }}_{D}}_h\right) }{\partial \underline{{\varvec{\varepsilon }}}_h}\nonumber \\ =&K_n\text {tr}^{+}(\underline{{\varvec{\varepsilon }}}_h)\underline{\varvec{1}}+2\mu \underline{{\varvec{\varepsilon }}_{D}}_h\end{aligned}$$

(117)

is the tensile/deviatoric part of the stress tensor. The derivative of the stresses with respect to the phase field reads

$$\begin{aligned} \frac{\partial \underline{\varvec{\sigma }}_h}{\partial {s_h}}=2s_h\left( K_n\text {tr}^{+}(\underline{{\varvec{\varepsilon }}}_h)\underline{\varvec{1}}+2\mu \underline{{\varvec{\varepsilon }}_{D}}_h\right) . \end{aligned}$$

(118)

Similar to the element stiffness matrix it is possible to write the mass matrix as

$$\begin{aligned} \underline{\varvec{M}}_e= \left[ \begin{array}{l@{\quad }l@{\quad }l@{\quad }l} \underline{\varvec{M}}_{11,e} &{}\underline{\varvec{M}}_{12,e} &{}\underline{\varvec{M}}_{13,e} &{}\underline{\varvec{M}}_{14,e}\\ \underline{\varvec{M}}_{21,e} &{}\underline{\varvec{M}}_{22,e} &{}\underline{\varvec{M}}_{23,e} &{}\underline{\varvec{M}}_{24,e}\\ \underline{\varvec{M}}_{31,e} &{}\underline{\varvec{M}}_{32,e} &{}\underline{\varvec{M}}_{33,e} &{}\underline{\varvec{M}}_{34,e}\\ \underline{\varvec{M}}_{41,e} &{}\underline{\varvec{M}}_{42,e} &{}\underline{\varvec{M}}_{43,e} &{}\underline{\varvec{M}}_{44,e} \end{array} \right] , \end{aligned}$$

(119)

with submatrices

$$\begin{aligned} \underline{\varvec{M}}_{IJ,e}= \left[ \begin{array}{l@{\quad }l} \underline{\varvec{M}}_{IJ,e}^{\varvec{u}\varvec{u}} &{}\underline{\varvec{M}}_{IJ,e}^{\varvec{u}{s}} \\ \underline{\varvec{M}}_{IJ,e}^{s\varvec{u}} &{}\underline{M}_{IJ,e}^{ss} \end{array} \right] . \end{aligned}$$

(120)

The components of the submatrices follow from

$$\begin{aligned} \underline{\varvec{M}}_{IJ,e}^{\varvec{u}\varvec{u}}&= \frac{\partial \underline{\varvec{P}}_{I,e}^{\varvec{u}}}{\partial \underline{\ddot{\varvec{u}}}_J}\nonumber \\&=\int \limits _{\varOmega _e}{N_I\rho \frac{\partial \ddot{\underline{\varvec{u}}}_h}{\partial \underline{\ddot{\varvec{u}}}_J}\mathrm{d}V\varvec{1}}\nonumber \\&=\int \limits _{\varOmega _e}{N_I\rho {N_J}\mathrm{d}V}\varvec{1},\end{aligned}$$

(121)

$$\begin{aligned} \underline{\varvec{M}}_{IJ,e}^{\varvec{u}{s}}&=\frac{\partial \underline{\varvec{P}}_{I,e}^{\varvec{u}}}{\partial {\ddot{s}_J}}=\underline{\varvec{0}},\end{aligned}$$

(122)

$$\begin{aligned} \underline{\varvec{M}}_{IJ,e}^{s\varvec{u}}&=\frac{\partial {{P}}_{I,e}^{s}}{\partial \underline{\ddot{\varvec{u}}}_J}=\varvec{\underline{0}},\end{aligned}$$

(123)

$$\begin{aligned} \underline{\varvec{M}}_{IJ,e}^{ss}&=\frac{\partial {{P}}_{I,e}^{s}}{\partial \ddot{s}_J}=0. \end{aligned}$$

(124)

Finally the element tangent matrix can be formulated

$$\begin{aligned} \underline{\varvec{S}}_e^{(k)}=\underline{\varvec{K}}_e^{(k)}+\frac{\gamma }{\beta \Delta {t}}\underline{\varvec{D}}_e^{(k)}+\frac{1}{\beta \Delta {t}^2}\underline{\varvec{M}}_e^{(k)} \end{aligned}$$

(125)

and the assembly of all element tangent matrix yields the global tangent matrix

$$\begin{aligned} \underline{\varvec{S}}_{n+1}^{(k)} = \bigcup _{e=1}^{n_e}{\underline{\varvec{S}}_e^{(k)}}. \end{aligned}$$

(126)

1.2 Irreversibility

As proposed in Kuhn [29], the irreversibility of the crack growth is modelled by defining homogeneous Dirichlet boundary conditions for the phase field $s(\varvec{x},t)$ once the value $s=0$ is reached for the first time i.e.

$$\begin{aligned} s_{I,n+1}=0 \ \ \text {if} \ \ s_{I,n}=0. \end{aligned}$$

(127)

This is accomplished by a reformulation of the residual and the tangent matrix on element level rather than changing global boundary conditions in order to achieve

$$\begin{aligned} \Delta {s}_{I,n+1}^{(k)}=-{s}_{I,n+1}^{(k)}. \end{aligned}$$

(128)

This causes the phase field to remain zero in subsequent time steps

$$\begin{aligned} s_{I,n+1}^{(k+1)}=s_{I,n+1}^{(k)}+\Delta {s}_{I,n+1}^{(k)} =s_{I,n+1}^{(k)}-{s}_{I,n+1}^{(k)}=0. \end{aligned}$$

(129)

The result (128) is enforced by reformulation of (103)

$$\begin{aligned} \underline{\varvec{R}}_{e,n+1}^{(k)}=\underline{\varvec{S}}_{e,n+1}^{(k)}\Delta \underline{\varvec{d}}_{e,n+1}^{(k)} \end{aligned}$$

(130)

on the element level.

• In a first step, the column ${\underline{\varvec{S}}_{e,n+1}^{(k)}}^{[:,s_I]}$ of the element tangent matrix is multiplied by ${s}_{I,n+1}^{(k)}$ and added to the residual $\underline{\varvec{R}}_{e,n+1}^{(k)}\leftarrow \underline{\varvec{R}}_{e,n+1}^{(k)}+{s}_{I,n+1}^{(k)}{\underline{\varvec{S}}_{e,n+1}^{(k)}}^{[:,s_I]}.$

• Then the corresponding residual entry to $s_I$ is replaced by $-{s}_{I,n+1}^{(k)}$, i.e. ${{R}_{e,n+1}^{(k)}}^{[s_I]}~=~-~{s}_{I,n+1}^{(k)}$.

• Finally, the row ${\underline{\varvec{S}}_{e,n+1}^{(k)}}^{[s_I,:]}$ and the column ${\underline{\varvec{S}}_{e,n+1}^{(k)}}^{[:,s_I]}$ are set to zero and the entry ${\underline{\varvec{S}}_{e,n+1}^{(k)}}^{[s_I,s_I]}$ is set to ${\underline{\varvec{S}}_{e,n+1}^{(k)}}^{[s_I,s_I]} = 1.$