An adaptive space-time phase field formulation for dynamic fracture of brittle shells based on LR NURBS

Abstract

We present an adaptive space-time phase field formulation for dynamic fracture of brittle shells. Their deformation is characterized by the Kirchhoff–Love thin shell theory using a curvilinear surface description. All kinematical objects are defined on the shell’s mid-plane. The evolution equation for the phase field is determined by the minimization of an energy functional based on Griffith’s theory of brittle fracture. Membrane and bending contributions to the fracture process are modeled separately and a thickness integration is established for the latter. The coupled system consists of two nonlinear fourth-order PDEs and all quantities are defined on an evolving two-dimensional manifold. Since the weak form requires $C^1$-continuity, isogeometric shape functions are used. The mesh is adaptively refined based on the phase field using Locally Refinable (LR) NURBS. Time is discretized based on a generalized-$\alpha$ method using adaptive time-stepping, and the discretized coupled system is solved with a monolithic Newton–Raphson scheme. The interaction between surface deformation and crack evolution is demonstrated by several numerical examples showing dynamic crack propagation and branching.

Appendices

Appendix

Time integration scheme

The system in Eq. (87) with intermediate quantities and the quantities at time step $n+1$

$$\begin{aligned} {\mathbf {x}}_{n+1}= & {} {\mathbf {x}}_n+\Delta t\, \dot{{\mathbf {x}}}_n+\big (\big (0.5 -\beta \big )\Delta t^2 \big )\ddot{{\mathbf {x}}}_n+\beta \Delta t^2 \ddot{{\mathbf {x}}}_{n+1}\,, \nonumber \\ {\dot{{\mathbf {x}}}}_{n+1}= & {} {\dot{{\mathbf {x}}}}_n+\big (\big (1 -\gamma \big )\Delta t \big )\ddot{{\mathbf {x}}}_n+\gamma \Delta t\ddot{{\mathbf {x}}}_{n+1}\,, \nonumber \\ {\mathbf {x}}_{n+\alpha _\mathrm {f}}= & {} \big (1-\alpha _\mathrm {f}\big ){\mathbf {x}}_n+\alpha _\mathrm {f}{\mathbf {x}}_{n+1}\,, \nonumber \\ {\dot{{\mathbf {x}}}}_{n+\alpha _\mathrm {f}}= & {} \big (1-\alpha _\mathrm {f}\big ){\dot{{\mathbf {x}}}}_n+\alpha _\mathrm {f}{\dot{{\mathbf {x}}}}_{n+1}\,, \nonumber \\ \ddot{\mathbf {x}}_{n+\alpha _\mathrm {m}}= & {} \big (1-\alpha _\mathrm {m}\big )\ddot{\mathbf {x}}_n+\alpha _\mathrm {m}\ddot{\mathbf {x}}_{n+1}\,, \end{aligned}$$

(98)

has to be solved. Here, $\Delta t = t_{n+1}-t_n$ refers to the time step size. Numerical dissipation is controlled by the parameters $\gamma$, $\beta$, $\alpha _\mathrm {f}$ and $\alpha _\mathrm {m}$. They are expressed in terms of $\rho _{\infty }\in [0,1]$, which resembles an algorithmic parameter that corresponds to the spectral radius of the amplification matrix as $\Delta t \rightarrow \infty$ (see Chung and Hulbert [16] for further details), i.e.

$$\begin{aligned} \alpha _\mathrm {f}= & {} \displaystyle \frac{1}{1+\rho _{\infty }}\,,\quad \alpha _\mathrm {m}= \displaystyle \frac{2-\rho _{\infty }}{1+\rho _{\infty }}\,,\nonumber \\ \gamma= & {} \displaystyle \frac{1}{2}+\alpha _\mathrm {m}-\alpha _\mathrm {f}\,,\quad \beta = \displaystyle \frac{1}{4}\,(1+\alpha _\mathrm {m}-\alpha _\mathrm {f})^2\,. \end{aligned}$$

(99)

We have found $\rho _\infty = 0.5$ to be a good choice and have used this in all computations. To solve the nonlinear system of equations in Eq. (87) using the Newton–Raphson procedure, it has to be linearized, i.e.

$$\begin{aligned} \displaystyle \begin{bmatrix} {\mathbf {K}}_\mathrm {x}&{} {\mathbf {K}}_\phi \\ {{\bar{{\mathbf {K}}}}}_\mathrm {x}&{} {{\bar{{\mathbf {K}}}}}_\phi \end{bmatrix} \begin{bmatrix} \Delta {\mathbf {x}}_{n+1} \\ \Delta \varvec{\phi }_{n+1} \end{bmatrix}= & {} \displaystyle - \begin{bmatrix} {\mathbf {f}}\left( {\mathbf {x}}_{n+\alpha _\mathrm {f}},\ddot{{\mathbf {x}}}_{n+\alpha _\mathrm {m}},\varvec{\phi }_{n+1} \right) \\ {{\bar{{\mathbf {f}}}}}\left( {\mathbf {x}}_{n+\alpha _\mathrm {f}},\varvec{\phi }_{n+1} \right) \end{bmatrix}, \nonumber \\ \end{aligned}$$

(100)

where the tangent matrix blocks are computed from

$$\begin{aligned} {\mathbf {K}}_\mathrm {x}= & {} \displaystyle \frac{\partial {\mathbf {f}}}{\partial {\mathbf {x}}_{n+1}} = \displaystyle \alpha _\mathrm {f}\frac{\partial {\mathbf {f}}}{\partial {\mathbf {x}}_{n+\alpha _\mathrm {f}}} + \frac{\alpha _\mathrm {m}}{\beta \Delta t^2} \frac{\partial {\mathbf {f}}}{\partial \ddot{\mathbf {x}}_{n+\alpha _\mathrm {f}}}\,, \nonumber \\ {\mathbf {K}}_\phi= & {} \displaystyle \frac{\partial {\mathbf {f}}}{\partial \varvec{\phi }_{n+1}}\,,\nonumber \\ {{\bar{{\mathbf {K}}}}}_\mathrm {x}= & {} \displaystyle \frac{\partial {{\bar{{\mathbf {f}}}}}}{\partial {\mathbf {x}}_{n+1}} = \alpha _\mathrm {f}\frac{\partial {{\bar{{\mathbf {f}}}}}}{\partial {\mathbf {x}}_{n+\alpha _\mathrm {f}}}\,,\nonumber \\ {{\bar{{\mathbf {K}}}}}_\phi= & {} \displaystyle \frac{\partial {{\bar{{\mathbf {f}}}}}}{\partial \varvec{\phi }_{n+1}}\,. \end{aligned}$$

(101)

The required linearizations of the force vectors are shown in Appendix B:. The initial guess for the Newton–Raphson iteration is set to

$$\begin{aligned}&{\mathbf {x}}_{n+1}^0 = {\mathbf {x}}_n+\Delta t\, \dot{{\mathbf {x}}}_n+\big (\big (0.5 -\beta \big )\Delta t^2 \big )\ddot{{\mathbf {x}}}_n+\big (\beta \Delta t^2 \big )\ddot{{\mathbf {x}}}^0_{n+1}\,, \nonumber \\&{\dot{{\mathbf {x}}}}_{n+1}^0 = {\dot{{\mathbf {x}}}}_n\,,\nonumber \\&\ddot{\mathbf {x}}_{n+1}^0 = \ddot{\mathbf {x}}_n\displaystyle \frac{\gamma -1}{\gamma }\,, \nonumber \\&\varvec{\phi }_{n+1}^0 = \varvec{\phi }_n\,, \end{aligned}$$

(102)

and then updated from iteration step $i\rightarrow i+1$ by

$$\begin{aligned} {\mathbf {x}}_{n+1}^{i+1}= & {} {\mathbf {x}}_{n+1}^{i} + \Delta {\mathbf {x}}_{n+1}^{i+1}~, \nonumber \\ {\dot{{\mathbf {x}}}}_{n+1}^{i+1}= & {} {\dot{{\mathbf {x}}}}_{n+1}^{i} + \Delta {\mathbf {x}}_{n+1}^{i+1} \displaystyle \frac{1}{\gamma \,\Delta t}\,, \nonumber \\ \ddot{\mathbf {x}}_{n+1}^{i+1}= & {} \ddot{\mathbf {x}}_{n+1}^{i} + \Delta {\mathbf {x}}_{n+1}^{i+1} \displaystyle \frac{1}{\beta \,\Delta t^2}\,, \nonumber \\ \varvec{\phi }_{n+1}^{i+1}= & {} \varvec{\phi }_{n+1}^{i}+\Delta \varvec{\phi }_{n+1}^{i+1}\,, \end{aligned}$$

(103)

until convergence is achieved. At iteration i, we check for the two convergence criteria

$$\begin{aligned} \max \displaystyle \left\{ \frac{\Vert {\mathbf {f}}^i_{n+1}\Vert }{\Vert {\mathbf {f}}^0_{n+1}\Vert },\frac{\Vert {{\bar{{\mathbf {f}}}}}^i_{n+1}\Vert }{\Vert {{\bar{{\mathbf {f}}}}}^0_{n+1}\Vert } \right\} \le \text {tol}^\mathrm {dyn}\,, \end{aligned}$$

(104)

with $\Vert ...\Vert$ denoting the Euclidean norm and $\text {tol}^\mathrm {dyn}=10^{-4}$ and

$$\begin{aligned} \begin{bmatrix}{\mathbf {f}}\\ {{\bar{{\mathbf {f}}}}} \end{bmatrix} \cdot \begin{bmatrix}\Delta {\mathbf {x}}\\ \Delta \varvec{\phi }\end{bmatrix} \le \text {tol}^\mathrm {nrg}\,, \end{aligned}$$

(105)

with $\text {tol}^\mathrm {nrg}=10^{-25}$.

Linearization

This section presents the respective elemental contributions for the tangent blocks in Eq. (101). The linearization of the mechanical force vector ${\mathbf {f}}^e:={\mathbf {f}}^e_\mathrm {kin}+{\mathbf {f}}^e_\mathrm {int}-{\mathbf {f}}^e_\mathrm {ext}$ of finite element $\Omega ^e$ with respect to the respective nodal positions ${\mathbf {x}}_e$ can be found in the work of Duong et al. [19]. Since we model the pressure as a function of the phase field variable, we need to linearize the external force vector with respect to $\phi$. This linearization of the pressure part ${\mathbf {f}}_{\mathrm {ext}p}^e$ of the external elemental force vector reads

$$\begin{aligned} \Delta _\phi \,{\mathbf {f}}_{\mathrm {ext}p}^e := \displaystyle \int _{\Omega ^e}{\mathbf {N}}^\mathrm{T}\,{{\bar{p}}}\,{\varvec{n}}^h\,{{\bar{{\mathbf {N}}}}}\,\mathrm {d}a\,\Delta \varvec{\phi }_e\,. \end{aligned}$$

(106)

For the linearization of the internal force vector, the four material tangents

$$\begin{aligned} \begin{aligned} c^{\alpha \beta \gamma \delta }&:= 2\displaystyle \frac{\partial {\tau ^{\alpha \beta }}}{\partial {a_{\gamma \delta }}}\,,\quad&d^{\alpha \beta \gamma \delta }&:= \displaystyle \frac{\partial {\tau ^{\alpha \beta }}}{\partial {b_{\gamma \delta }}}\,,\\ e^{\alpha \beta \gamma \delta }&:= 2\displaystyle \frac{\partial {M_0^{\alpha \beta }}}{\partial {a_{\gamma \delta }}}\,,\quad&f^{\alpha \beta \gamma \delta }&:= \displaystyle \frac{\partial {M_0^{\alpha \beta }}}{\partial {b_{\gamma \delta }}}\,,\\ \end{aligned} \end{aligned}$$

(107)

have to be defined. Since we assume the constitutive in-plane response to be fully decoupled from the out-of-plane response, it follows that $d^{\alpha \beta \gamma \delta }=e^{\alpha \beta \gamma \delta }=0$. According to Eqs. (49) and (50), the first tangent matrix can be computed based on the contributions

$$\begin{aligned} \displaystyle \frac{\partial {\tau _\mathrm {dil}^{\alpha \beta }}}{\partial {a_{\gamma \delta }}}= & {} \dfrac{K}{2}\Bigr (J^2a^{\alpha \beta }a^{\gamma \delta }+\bigl (J^2-1\bigr )\,a^{\alpha \beta \gamma \delta }\Bigr )\,,\nonumber \\ \displaystyle \frac{\partial {\tau _\mathrm {dev}^{\alpha \beta }}}{\partial {a_{\gamma \delta }}}= & {} \dfrac{G}{2J}\left( \dfrac{I_1}{2}a^{\alpha \beta }a^{\gamma \delta }-I_1 a^{{\alpha \beta }{\gamma \delta }}-a^{\alpha \beta }A^{\gamma \delta }- A^{\alpha \beta }a^{\gamma \delta }\right) .\nonumber \\ \end{aligned}$$

(108)

Based on Eqs. (52) and (53), the tangent matrix $f^{{\alpha \beta }{\gamma \delta }}$ can be computed from the contribution

$$\begin{aligned} \displaystyle \frac{\partial ^2{{\tilde{\Psi }}_\mathrm {bend}(\xi )}}{\partial {b_{\alpha \beta }}\,\partial {b_{\gamma \delta }}}=\xi ^2\dfrac{12}{T^3}\,c\,A^{\alpha \gamma }A^{ \beta \delta }\,. \end{aligned}$$

(109)

Since we consider the fully linearized system in Eq. (100), we also need to linearize the mechanical force vector with respect to the phase field, i.e.

$$\begin{aligned} \Delta _\phi {\mathbf {f}}^e = \big [{{\mathbf {k}}}^e_{\sigma \phi }+{{\mathbf {k}}}^e_{M\phi }\big ]\,\Delta \varvec{\phi }_e\,, \end{aligned}$$

(110)

with

$$\begin{aligned} {{\mathbf {k}}}^e_{\sigma \phi }:= & {} \displaystyle \int _{\Omega _0^e}g'(\phi )\,\tau ^{\alpha \beta }_+\,{\mathbf {N}}^\mathrm{T}_{\!,\alpha }\,{\varvec{a}}^h_\beta \,{{\bar{{\mathbf {N}}}}}\,\mathrm {d}A\,, \nonumber \\ {{\mathbf {k}}}^e_{M\phi }:= & {} \displaystyle \int _{\Omega _0^e}g'(\phi )\,M_{0,+}^{\alpha \beta }\,{\mathbf {N}}^\mathrm{T}_{\!;\alpha \beta }\,{\varvec{n}}^h\,{{\bar{{\mathbf {N}}}}}\,\mathrm {d}A\,, \end{aligned}$$

(111)

where $\tau ^{\alpha \beta }:=J\sigma ^{\alpha \beta }$ and $M^{\alpha \beta }_0:=JM^{\alpha \beta }$ has been used to map the integrals to the element domain in the reference configuration. According to Eq. (86), the linearization of ${{\bar{{\mathbf {f}}}}}^e$ with respect to the respective nodal positions ${\mathbf {x}}_e$ yields

$$\begin{aligned} \Delta _\mathrm {x}{{\bar{{\mathbf {f}}}}}^e_\mathrm {el} := \displaystyle \int _{\Omega _0^e}{{\bar{{\mathbf {N}}}}}^\mathrm{T}\displaystyle \frac{2\ell _0}{{\mathcal {G}}_\mathrm {c}}g'(\phi )\,\Delta _\mathrm {x}{\mathcal {H}}\,\mathrm {d}A\,\Delta {\mathbf {x}}_\mathrm {e}\,, \end{aligned}$$

(112)

with

$$\begin{aligned} \Delta _\mathrm {x}{\mathcal {H}}:= \Delta _\mathrm {x}\max \limits _{\tau \in [0,t]}\Psi _\mathrm {el}^{+}({\varvec{x}},\tau )\,, \end{aligned}$$

(113)

and

$$\begin{aligned} \Delta _\mathrm {x}\Psi _\mathrm {el}^{+}:= {\tau ^{\alpha \beta }_{\mathrm {el},+}}\,{\varvec{a}}_\alpha \cdot {\mathbf {N}}_{\!,\beta } + {M^{\alpha \beta }_{0,+}}\,{\varvec{n}}\cdot {\mathbf {N}}_{\!;\alpha \beta } \,. \end{aligned}$$

(114)

The linearization of ${\bar{{\mathbf {f}}}}_\mathrm {int}^e$ with respect to the phase field variables of $\Omega ^e$ reads

$$\begin{aligned} \Delta _\phi {{\bar{{\mathbf {f}}}}}^e_\mathrm {int} := \Bigl [{{\bar{{\mathbf {k}}}}}^e_0 + {{\bar{{\mathbf {k}}}}}^e_\mathrm {el}\Bigr ]\,\Delta \varvec{\phi }_\mathrm {e}\,, \end{aligned}$$

(115)

with

$$\begin{aligned} {{\bar{{\mathbf {k}}}}}^e_\mathrm {el} := \displaystyle \int _{\Omega _0^e}{{\bar{{\mathbf {N}}}}}^\mathrm{T}\bigg (\displaystyle \frac{2\ell _0}{{\mathcal {G}}_\mathrm {c}}g''(\phi ){\mathcal {H}}\bigg ){{\bar{{\mathbf {N}}}}}\,\,\mathrm {d}A\,. \end{aligned}$$

(116)

The matrices ${{\bar{{\mathbf {k}}}}}^e_0$ and ${{\bar{{\mathbf {k}}}}}^e_\mathrm {el}$ both contribute to the tangent block ${{\bar{{\mathbf {K}}}}}_\phi$ in Eq. (100).