Phase-field modelling of brittle fracture in thin shell elements based on the MITC4+ approach

Abstract

We present a phase field based MITC4+ shell element formulation to simulate fracture propagation in thin shell structures. The employed MITC4+ approach renders the element shear- and membrane- locking free, hence providing high-fidelity fracture simulations in planar and curved topologies. To capture the mechanical response under bending-dominated fracture, a crack-driving force description based on the maximum strain energy density through the shell-thickness is considered. Several numerical examples simulating fracture in flat and curved shell structures are presented, and the accuracy of the proposed formulation is examined by comparing the predicted critical fracture loads against analytical estimates.

Appendices

Appendix A: Jacobian for coordinate transformation

The Jacobian $[{\mathbf {J}}]$ for coordinate transformation mapping in a Reissner–Mindlin shell element and its first column are defined as in Eqs. (61) and (62). Eq. (62) can be subsequently used to derive expressions for second and third column in a similar manner.

$$\begin{aligned}{}[\mathbf{J}] = \begin{bmatrix} x_{,\xi } &{}\quad y_{,\xi } &{}\quad z_{,\xi } \\ x_{,\eta } &{}\quad y_{,\eta } &{}\quad z_{,\eta } \\ x_{,\zeta } &{}\quad y_{,\zeta } &{}\quad z_{,\zeta } \end{bmatrix} \end{aligned}$$

(61)

where

$$\begin{aligned} \begin{bmatrix} x_{,\xi } \\ x_{,\eta } \\ x_{,\zeta } \end{bmatrix} = \begin{bmatrix} \sum N_{i,\xi } \Big (x_i+\dfrac{\zeta t_i l_{3i}}{2}\Big ) \\ \sum N_{i,\eta } \Big (x_i+\dfrac{\zeta t_i l_{3i}}{2}\Big ) \\ \sum N_i \Big (\dfrac{t_i l_{3i}}{2}\Big ) \end{bmatrix} \end{aligned}$$

(62)

where ${\mathbf {x}}=[x,y,z]$ is the position vector of any arbitrary point within the shell element, $\{\xi , \eta , \zeta \}$ are the shell parametric coordinates, $t_i$ is the shell thickness and $\{l_{3i}, m_{3i}, n_{3i}\}$ are the direction cosines of normal vector ${\mathbf {V}}_{3i}$ to the shell mid-surface at any node i.

Appendix B: Coordinate-transformation matrix for rotation of strain tensors

The strains can be rotated from any one coordinate system (say $C_1$ with normalized basis vectors $\bar{e}$) to another coordinate system ($C_2$ with normalized basis vectors ${\hat{e}}$) by via the strain-transformation matrix ${\pmb {\mathcal {T}}}_{{\varvec{\varepsilon }}}$ shown in Eq. (63).

$$\begin{aligned} {\pmb {\mathcal {T}}}_{{\varvec{\varepsilon }}}=\begin{bmatrix} {\pmb {\mathcal {T}}}_{11} &{}\quad {\pmb {\mathcal {T}}}_{12} \\ {\pmb {\mathcal {T}}}_{21} &{}\quad {\pmb {\mathcal {T}}}_{22} \end{bmatrix} \end{aligned}$$

(63)

with,

$$\begin{aligned} {\pmb {\mathcal {T}}}_{11}= & {} \begin{bmatrix} l^2_1 &{}\quad m^2_1 &{}\quad n^2_1 \\ l^2_2 &{}\quad m^2_2 &{}\quad n^2_2 \\ l^2_3 &{}\quad m^2_3 &{}\quad n^2_3 \end{bmatrix} \end{aligned}$$

(64)

$$\begin{aligned} {\pmb {\mathcal {T}}}_{12}= & {} \begin{bmatrix} l_1 m_1 &{}\quad m_1 n_1 &{}\quad n_1 l_1 \\ l_2 m_2 &{}\quad m_2 n_2 &{}\quad n_2 l_2 \\ l_3 m_3 &{}\quad m_3 n_3 &{}\quad n_3 l_3 \\ \end{bmatrix} \end{aligned}$$

(65)

$$\begin{aligned} {\pmb {\mathcal {T}}}_{21}= & {} \begin{bmatrix} 2 l_1 l_2 &{}\quad 2 m_1 m_2 &{}\quad 2 n_1 n_2 \\ 2 l_2 l_3 &{}\quad 2 m_2 m_3 &{}\quad 2 n_2 n_3 \\ 2 l_3 l_1 &{}\quad 2 m_3 m_1 &{}\quad 2 n_3 n_1 \end{bmatrix} \end{aligned}$$

(66)

$$\begin{aligned} {\pmb {\mathcal {T}}}_{22}= & {} \begin{bmatrix} l_1 m_2 + l_2 m_1 &{}\quad m_1 n_2 + m_2 n_1 &{} \quad n_1 l_2 + n_2 l_1 \\ l_2 m_3 + l_3 m_2 &{}\quad m_2 n_3 + m_3 n_2 &{}\quad n_2 l_3 + n_3 l_2 \\ l_3 m_1 + l_1 m_3 &{}\quad m_3 n_1 + m_1 n_3 &{}\quad n_3 l_1 + n_1 l_3 \end{bmatrix} \end{aligned}$$

(67)

where the terms $[l_1, m_1, n_1]$, $[l_2, m_2, n_2]$ and $[l_3, m_3, n_3]$ correspond to the direction cosines of the shell nodal-vectors $\mathbf{V}_{1i}$, $\mathbf{V}_{2i}$ and $\mathbf{V}_{3i}$ respectively, defined according to Eq. (68) [11].

$$\begin{aligned} \begin{aligned} l_1&= \text{ cos }[\bar{\mathbf{e}}_{x},{\hat{\mathbf{e}}}_{x}] \, ; \, m_1 = \text{ cos }[\bar{\mathbf{e}}_{y},{\hat{\mathbf{e}}}_{x}] \, ; \, n_1 = \text{ cos }[\bar{\mathbf{e}}_{z},{\hat{\mathbf{e}}}_{x}] \\ l_2&= \text{ cos }[\bar{\mathbf{e}}_{x},{\hat{\mathbf{e}}}_{y}] \, ; \, m_2 = \text{ cos }[\bar{\mathbf{e}}_{y},{\hat{\mathbf{e}}}_{y}] \, ; \, n_2 = \text{ cos }[\bar{\mathbf{e}}_{z},{\hat{\mathbf{e}}}_{y}] \\ l_3&= \text{ cos }[\bar{\mathbf{e}}_{x},{\hat{\mathbf{e}}}_{z}] \, ; \, m_3 = \text{ cos }[\bar{\mathbf{e}}_{y},{\hat{\mathbf{e}}}_{z}] \, ; \, n_3 = \text{ cos }[\bar{\mathbf{e}}_{z},{\hat{\mathbf{e}}}_{z}] \end{aligned} \end{aligned}$$

(68)

The resulting ${\pmb {\mathcal {T}}}_{{\varvec{\varepsilon }}}$ is a $(6\times 6)$ matrix which can be multiplied to $(6\times 1)$ strain vector (expressed in Voigt notation) to transform it from coordinate system $C_1$ to coordinate system $C_2$.