Progressive delamination analysis of composite materials using XFEM and a discrete damage zone model

Abstract

The modeling of progressive delamination by means of a discrete damage zone model within the extended finite element method is investigated. This framework allows for both bulk and interface damages to be conveniently traced, regardless of the underlying mesh alignment. For discrete interfaces, a new mixed-mode force–separation relation, which accounts for the coupled interaction between opening and sliding modes, is proposed. The model is based on the concept of Continuum Damage Mechanics and is shown to be thermodynamically consistent. An integral-type nonlocal damage is adopted in the bulk to regularize the softening material response. The resulting nonlinear equations are solved using a Newton scheme with a dissipation-based arc-length constraint, for which an analytical Jacobian is derived. Several benchmark delamination studies, as well as failure analyses of a fiber/epoxy unit cell, are presented and discussed in detail. The proposed model is validated against available analytical/experimental data and is found to be robust and mesh insensitive.

Appendices

Appendix 1

For completeness, we provide the proof in this appendix that the power law criterion given in (33) can be recovered by using the proposed mixed-mode force–separation relation (31), if the discrete interface is subjected to a deformation history with a constant ratio $\beta = \delta _n /\delta _t$.

Without loss of generality, we only consider a monotonic loading with $\delta _n \ge 0$. Then the history ratio (29) can be rewritten as

$$\begin{aligned} \zeta&= \left[ \left( \dfrac{\delta _{n}}{\delta ^{\mathrm {cr}}_n} \right) ^{\alpha } + \left( \dfrac{| \delta _{t} |}{\delta ^{\mathrm {cr}}_t} \right) ^{\alpha } \right] ^{1/\alpha }\nonumber \\&= \delta _{n} \left[ \left( \dfrac{1}{\delta ^{\mathrm {cr}}_n} \right) ^{\alpha } + \left( \dfrac{1}{|\beta |\delta ^{\mathrm {cr}}_t} \right) ^{\alpha } \right] ^{1/\alpha } \end{aligned}$$

(80)

Accoding to the mixed-mode damage law (30), the critical separation $\bar{\delta }_i$ corresponding to the damage initiation can be calculated as

$$\begin{aligned} \bar{\delta }_n=\frac{1}{\eta },\,\,\,\,\bar{\delta }_t=\frac{1}{\beta \eta } \end{aligned}$$

(81)

where

$$\begin{aligned} \eta =\left[ \left( \dfrac{1}{\delta ^{\mathrm {cr}}_n} \right) ^{\alpha } + \left( \dfrac{1}{|\beta |\delta ^{\mathrm {cr}}_t} \right) ^{\alpha } \right] ^{1/\alpha } \end{aligned}$$

(82)

Given the proposed mixed-mode force–separation relation (31), the energy release rate components of mode I and mode II due to the complete decohesion are obtained by splitting the integrals (44) into two parts:

$$\begin{aligned} G_{\mathrm {I}}&= \int _{0}^{\infty } F_{n}(\delta _n, \delta _t) \, d\delta _{n} /l_{s} \nonumber \\&= \left[ \int _{0}^{\frac{1}{\eta }} K_{n}^{0} \delta _n \, d\delta _n + \int ^{\infty }_{\frac{1}{\eta }} \dfrac{K_{n}^{0} \delta _n}{\exp (\eta \delta _n -1)} \, d\delta _n \right] \Bigg /l_s \nonumber \\&= \left( \dfrac{K_n^0}{2 \eta ^2} + \dfrac{2K_n^0}{\eta ^2} \right) /l_s = \dfrac{2.5K_n^0}{\eta ^2 l_s}\end{aligned}$$

(83)

$$\begin{aligned} G_{\mathrm {II}}&= \int _{0}^{\infty } F_{t}(\delta _n, \delta _t) \, d\delta _{t} /l_{s} \nonumber \\&= \left[ \int _{0}^{\frac{1}{\beta \eta }} K_{t}^{0} \delta _t \, d\delta _t + \int ^{\infty }_{\frac{1}{\beta \eta }} \dfrac{K_{t}^{0} \delta _t}{\exp (\beta \eta \delta _t -1)} \, d\delta _t \right] \Bigg /l_s \nonumber \\&= \left( \dfrac{K_t^0}{2 \beta ^2 \eta ^2} + \dfrac{2K_t^0}{\beta ^2 \eta ^2} \right) /l_s = \dfrac{2.5K_t^0}{\beta ^2 \eta ^2 l_s} \end{aligned}$$

(84)

Considering the expressions for the initial stiffnesses (28) and damage coefficient (23) results:

$$\begin{aligned} \dfrac{G_{\mathrm {I}}}{G_{\mathrm {IC}}} = \dfrac{1}{\eta ^2 {\delta ^{\mathrm {cr}}_n}^2},\,\,\,\, \dfrac{G_{\mathrm {II}}}{G_{\mathrm {IIC}}} = \dfrac{1}{|\beta |^2 \eta ^2 {\delta ^{\mathrm {cr}}_t}^2} \end{aligned}$$

(85)

Finally, the proof is completed:

$$\begin{aligned} \left( \dfrac{G_{\mathrm {I}}}{G_{\mathrm {IC}}} \right) ^{\alpha /2} +\left( \dfrac{G_{\mathrm {II}}}{G_{\mathrm {IIC}}} \right) ^{\alpha /2}&= \dfrac{1}{\eta ^{\alpha }} \left[ \left( \dfrac{1}{\delta ^{\mathrm {cr}}_n} \right) ^{\alpha } + \left( \dfrac{1}{|\beta |\delta ^{\mathrm {cr}}_t} \right) ^{\alpha } \right] \nonumber \\&= \dfrac{\eta ^{\alpha }}{\eta ^{\alpha }} = 1 \end{aligned}$$

(86)

Appendix 2

Considering the bulk damage law (12) and the mixed-mode force–separation relation (31) for discrete interfaces, the partial derivatives appear in Eqs. (70)–(73) are computed explicitly as follows:

$$\begin{aligned} d^{\prime }&=\dfrac{\partial {D^{\Omega }}}{\partial {\kappa }} =\left\{ \begin{array}{lcc} 0 &{} \text{ if } &{} \kappa \le \epsilon ^{\mathrm {cr}} \\ B \exp (-B(\kappa -\epsilon ^{\mathrm {cr}})) &{} \text{ if } &{} \kappa > \epsilon ^{\mathrm {cr}} \end{array} \right. \end{aligned}$$

(87)

$$\begin{aligned} \kappa ^{\prime }&=\dfrac{\partial {\kappa }}{\partial {\tilde{\epsilon }_{eq}}} =\left\{ \begin{array}{ccc} 0 &{} \text{ if } &{} \text{ unloading } \\ 1 &{} \text{ if } &{} \text{ loading } \end{array} \right. \end{aligned}$$

(88)

The partial derivative $\mathbf {g}$ is equal to

$$\begin{aligned} \mathbf {g}=\dfrac{\partial {\tilde{\epsilon }_{eq}}}{\partial {\varvec{\epsilon }}}= \dfrac{\partial {\tilde{\epsilon }_{eq}}}{\partial {\varvec{e}}} \dfrac{\partial {\varvec{e}}}{\partial {\varvec{I}}} \dfrac{\partial {\varvec{I}}}{\partial {\varvec{\epsilon }}} \end{aligned}$$

(89)

where $\varvec{e}$ and $\varvec{I}$ are arranged into column vectors as follows:

$$\begin{aligned} \varvec{I}&=\{I_1\,I_2 \}^{\mathrm {T}}\,\,\,\,\, \text{ with }\, I_1=\epsilon _{xx}+\epsilon _{yy},\, I_2=\epsilon _{xx} \epsilon _{yy}-\gamma ^2_{xy}/4 \end{aligned}$$

(90)

$$\begin{aligned} \varvec{e}&=\{e_1\,e_2\,e_3 \}^{\mathrm {T}}\,\,\text{ with }\, e_1=\dfrac{I_1}{2}+\dfrac{\sqrt{{I_1}^2-4I_2}}{2},\nonumber \\ e_2&=\dfrac{I_1}{2}-\dfrac{\sqrt{{I_1}^2-4I_2}}{2},\, e_3=h\dfrac{\nu }{\nu -1}I_1 \end{aligned}$$

(91)

In the above equations, $\nu$ is the Poisson ratio, $h=1$ for plane-stress cases, and $h=0$ for plane-strain cases. Then the partial derivative of the Eq. (89) are given by

$$\begin{aligned} \dfrac{\partial {\tilde{\epsilon }_{eq}}}{\partial {\varvec{e}}}&= \left\{ \dfrac{\langle e_1 \rangle }{\sqrt{\langle e_1 \rangle ^2 +\langle e_2 \rangle ^2+\langle e_3 \rangle ^2}}\, \dfrac{\langle e_2 \rangle }{\sqrt{\langle e_1 \rangle ^2 +\langle e_2 \rangle ^2+\langle e_3 \rangle ^2}}\,\right. \nonumber \\&\qquad \left. \dfrac{\langle e_3 \rangle }{\sqrt{\langle e_1 \rangle ^2 +\langle e_2 \rangle ^2+\langle e_3 \rangle ^2}} \right\} \end{aligned}$$

(92)

$$\begin{aligned} \dfrac{\partial {\varvec{e}}}{\partial {\varvec{I}}}&= \begin{bmatrix} \dfrac{1}{2}+\dfrac{I_1}{2\sqrt{{I_1}^2-4I_2}}&- \dfrac{1}{\sqrt{{I_1}^2-4I_2}} \\ \dfrac{1}{2}-\dfrac{I_1}{2\sqrt{{I_1}^2-4I_2}}&\dfrac{1}{\sqrt{{I_1}^2-4I_2}} \\ h\dfrac{\nu }{\nu -1}&0 \end{bmatrix} \end{aligned}$$

(93)

$$\begin{aligned} \dfrac{\partial {\varvec{I}}}{\partial {\varvec{\epsilon }}}&= \begin{bmatrix} 1&1&0 \\ \epsilon _{yy}&\epsilon _{xx}&-\dfrac{\gamma _{xy}}{2} \end{bmatrix} \end{aligned}$$

(94)

Finally, the partial derivative $\mathbf {M}$ for $\delta _n \ge 0$ is evaluated by differentiating the Eq. (31) as

$$\begin{aligned} \mathbf {M}=\dfrac{\partial {\mathbf {F}^L}}{\partial {\varvec{\delta }^L}}= \left\{ \begin{array}{lcc} \mathbf {K}^0 &{} \text{ if } &{} \zeta \le 1 \\ \dfrac{\mathbf {K}^0}{\exp (\zeta -1)} &{} \text{ if } &{} \zeta > 1\, \text{(unloading) } \\ \dfrac{\mathbf {K}^0}{\exp (\zeta -1)} \left\{ \mathbf {I}-\zeta ^{1-\alpha } \begin{bmatrix} \dfrac{|\delta _t|^{\alpha }}{(\delta ^{\mathrm {cr}}_t)^{\alpha }} &{} \dfrac{\delta _t \delta _n^{\alpha -1}}{(\delta ^{\mathrm {cr}}_n)^{\alpha }}\\ \dfrac{\mathrm {sign}(\delta _t)\delta _n|\delta _t|^{\alpha -1}}{(\delta ^{\mathrm {cr}}_t)^{\alpha }} &{} \dfrac{\delta _n^{\alpha }}{(\delta ^{\mathrm {cr}}_n)^{\alpha }} \end{bmatrix} \right\} &{} \text{ if } &{} \zeta > 1\, \text{(loading) } \end{array} \right. \end{aligned}$$

(95)

with $\mathbf {I}$ the identity matrix.