A phase field approach for damage propagation in periodic microstructured materials

Abstract

In the present work, the evolution of damage in periodic composite materials is investigated through a novel finite element-based multiscale computational approach. The proposed methodology is developed by means of the original combination of asymptotic homogenization with the phase field approach to nonlocal damage. This last is applied at the macroscale level on the equivalent homogeneous continuum, whose constitutive properties are obtained in closed form via a two-scale asymptotic homogenization scheme. The formulation considers different assumptions on the evolution of damage at the microscale (e.g., damage in the matrix and not in the inclusion/fiber), as well as the role played by the microstructural reinforcement, i.e. its volumetric content and shape. Numerical results show that the proposed formulation leads to an apparent tensile strength and a post-peak branch of unnotched and notched specimens dependent not only on the internal length scale of the phase field approach, as for homogeneous materials, but also on microstructural features. Down-scaling relations provide the full reconstruction of the microscopic fields at any point of the macroscopic model, as a simple post-processing operation.

Appendices

Appendix A. Homogenization of a bi-phase elastic material: benchmark test

Fig. 15

Heterogeneous microstructured model and the equivalent homogenized one loaded by $\mathcal {L}$-periodic harmonic body forces $\mathbf {b}(x_1)$. The periodic cell $\mathcal {A}$ has microstructural characteristic size equal to $\varepsilon$ and volume fraction $f=1/4$

The first order asymptotic homogenization technique described in Sect. 3 is here validated for all the admissible values of the phase field variable $\mathfrak {d}$, with $0\le \mathfrak {d}\le 1$, in order to assess its capabilities to accurately describe the global behavior of composite elastic materials subjected to damage. The periodic cell $\mathcal {A}$ of the considered two-phase elastic material is the one depicted in Fig. 4, with matrix made by an Alluminum-like material and circular or square inclusion with Silicum carbide constitutive parameters. The elastic tensors of the two phases in plane strain conditions are reported in Eq. (29) in the absence of damage ($\mathfrak {d}=0$). Only components of the elastic tensor of the matrix are supposed to be affected by damage evolution inside the material through multiplication of the undamaged Young modulus $E_{m,0}=E_{Al}$ by degradation function $g(\mathfrak {d})=(1-\mathfrak {d})^2+\mathcal {K}$, with residual stiffness $\mathcal {K}=0.005$.

Fig. 16

Dimensionless micro displacement components $\tilde{u}_1$ and $\tilde{u}_2$ evaluated along the mean horizontal line of the heterogeneous model (dashed lines) and dimensionless macro displacement components $\tilde{U}_1$ and $\tilde{U}_2$ as solutions of the homogenized model (continuous lines) and as a result of the upscaling of the numerical solution of the heterogeneous model (squares). a$\tilde{U}_1$ and $\tilde{u}_1$ vs $\tilde{x}_1$ for $n_b=1$, b$\tilde{U}_1$ and $\tilde{u}_1$ vs $\tilde{x}_1$ for $n_b=2$, c$\tilde{U}_2$ and $\tilde{u}_2$ vs $\tilde{x}_1$ for $n_b=1$, and d$\tilde{U}_2$ and $\tilde{u}_2$ vs $\tilde{x}_1$ for $n_b=2$. Phase field variable $\mathfrak {d}=0$ (blue), $\mathfrak {d}=0.25$ (magenta), and $\mathfrak {d}=0.5$ (red)

Considering for example a volume fraction $f=1/4$, the periodic medium is loaded by means of $\mathcal {L}$-periodic body forces $\mathbf {b}(\mathbf {x})$ of the form

$$\begin{aligned} b_j(x_1)=B_j\, e^{(i2\pi n_bx_1/L)}, \end{aligned}$$

(31)

with $j=1,2$, wave number $n_b$ and L the macrostructural characteristic size, see Fig. 15. In Eq. (31) i represents the imaginary unit, such that $i^2=-1$. In view of the $\mathcal {L}$-periodicity of volume forces, only a representative portion of the entire heterogeneous material has been analyzed. In particular, because of the invariance of body forces $\mathbf {b}$ with respect to $x_2$, the model problem is composed by 11 cells along the $\mathbf {e}_1$ direction having a total dimension equal to L and one cell along the $\mathbf {e}_2$ direction.

Field equations for the first order homogenized material in terms of the overall elastic tensor $\mathbb {C}(\mathfrak {d})$ take the form

$$\begin{aligned} \frac{\partial }{\partial x_1}\left( C_{1j1j}(\mathfrak {d})\frac{\partial U_{j}}{\partial x_1} \right) +b_j(x_1)=0 \end{aligned}$$

(32)

whose analytical solution in terms of macro displacement $\mathbf {U}$ depends only upon $x_1$ in view of formula (31) and reads

$$\begin{aligned} U_1(x_1)= & {} \frac{B_1}{C_{1111}(\mathfrak {d})}\left( \frac{L}{2\pi }\right) ^2 \frac{1}{n_b^2}\,e^{(i2\pi n_bx_1/L)}\nonumber \\ U_2(x_1)= & {} \frac{B_2}{C_{1212}(\mathfrak {d})}\left( \frac{L}{2\pi }\right) ^2 \frac{1}{n_b^2}\,e^{(i 2\pi n_b x_1/L)} \end{aligned}$$

(33)

If one considers for example only the imaginary part of expressions (33), components of macro displacement field results

$$\begin{aligned} U_1(x_1)= & {} \frac{B_1}{C_{1111}(\mathfrak {d})}\left( \frac{L}{2\pi }\right) ^2 \frac{1}{n_b^2}\, \sin ({2\pi n_bx_1/L})\nonumber \\ U_2(x_1)= & {} \frac{B_2}{C_{1212}(\mathfrak {d})}\left( \frac{L}{2\pi }\right) ^2 \frac{1}{n_b^2}\, \sin ({ 2\pi n_b x_1/L}) \end{aligned}$$

(34)

Comparison between the homogenized model solution, as expressed in Eq. (34), and the solution obtained from a finite element analysis of the heterogeneous model subjected to periodic boundary conditions, is depicted in Fig. 16 for square inclusion and three different values of the phase field variable, namely $\mathfrak {d}=0,0.25,0.5$. Analogous results have been obtained for the case of circular inclusion. The heterogeneous macro solution has been computed from the corresponding microscopic one by means of up-scaling relation (14), by performing a mean of the micro field over each cell. Furthermore, for each value of $\mathfrak {d}$, components of the micro displacement field of the heterogeneous model are represented in Fig. 16, where the micro solution is evaluated at nodes having $x_2=0$, see Fig. 15. In particular, dimensionless components of the macro and micro displacement field, defined as

$$\begin{aligned}&\tilde{U}_1=\frac{U_1(x_1)}{L}, \; \tilde{U}_2=\frac{U_2(x_1)}{L}, \nonumber \\&\tilde{u}_1=\frac{u_1(x_1,x_2=0)}{L}, \tilde{u}_2=\frac{u_2(x_1,x_2=0)}{L}, \end{aligned}$$

(35)

are represented in Fig. 16 as functions of the dimensionless length $\tilde{x}_1=x_1/L.$ Unit dimensionless amplitudes

$$\begin{aligned} \tilde{B}_1=\frac{B_1 L}{C_{1111}(\mathfrak {d}=0)}=1, \; \tilde{B}_2=\frac{B_2 L}{C_{1111}(\mathfrak {d}=0)}=1, \end{aligned}$$

(36)

are considered in the analyses, where $C_{1111}(\mathfrak {d}=0)$ is a component of the fourth order overall elastic tensor evaluated for the undamaged material $(\mathfrak {d}=0)$. Capabilities of the first-order asymptotic homogenization technique applied for this class of periodic microstructured elastic materials are assessed by the good agreement obtained in all the analyzed cases between the numerical solution of the heterogeneous model and the solution of the homogenized one. Therefore, at the macroscale, the elastic behavior of a periodic heterogeneous medium can be accurately described by the present homogenization procedure through derivation of the overall constitutive properties of the equivalent Cauchy medium depending upon the value of the phase field parameter $\mathfrak {d}$ and the chosen degradation function $g(\mathfrak {d})$.

Appendix B. Finite element formulation of the coupled model

This appendix details the finite element formulation that has been exploited and implemented in the finite element software FEAP in order to solve the coupled system (27) in terms of the macro displacements $U_i(\mathbf {x})$ and the phase field variable $\mathfrak {d}(\mathbf {x})$.

The weak form of balance equations (27) of the coupled field problem detailed in Sect. 4, taking into account boundary conditions (28), reads

$$\begin{aligned}&- \int _{\varOmega }\frac{\partial \psi _{U_h}}{\partial x_k}\,C_{ijhk}\left( \mathfrak {d}\right) \frac{\partial U_i}{\partial x_j}\,\text {d}\varOmega + \int _{\Gamma _t} \psi _{U_h}\,t_h \,\text {d}\Gamma \nonumber \\&+ \int _{\varOmega }\psi _{U_h}\, b_h \, \text {d}\varOmega =0 \;\;\forall \psi _{U_h}\;s.t.\;\psi _{U_h}=0\;\mathrm {on}\;\Gamma _u \nonumber \\&\int _{\varOmega }\ell ^2\frac{\partial \psi _{\mathfrak {d}}}{\partial x_j}\,\frac{\partial \mathfrak {d}}{\partial x_j}\,\text {d}\varOmega + \int _{\varOmega }\psi _{\mathfrak {d}}\,\mathfrak {d}\,\text {d}\varOmega \nonumber \\&+ \int _{\varOmega }\frac{\ell }{2 G_C}\psi _\mathfrak {d} H_{ij}\,\frac{\partial C_{ijhk}(\mathfrak {d})}{\partial \mathfrak {d}}H_{hk}\,\text {d}\varOmega =0\;\; \forall \psi _{\mathfrak {d}}\nonumber \\ \end{aligned}$$

(37)

where $\psi _{U_h}$ and $\psi _{\mathfrak {d}}$ are taken as test functions. Considering the finite dimensional space $V_h$, for which $\{N_j|j=1,2,...,N_h\}$ is a basis, in the finite element discretization the macro displacement field $\mathbf {U}(\mathbf {x})$ and the phase field $\mathfrak {d}(\mathbf {x})$ are approximated as linear combinations of shape functions $N_j(\mathbf {x})$ and nodal unknowns $U_{ij}$ and $\mathfrak {d}_j$

$$\begin{aligned} U_i(\mathbf {x})=\sum _{j=1}^{N_h}N_j(\mathbf {x})U_{ij}, \quad \mathfrak {d}(\mathbf {x})=\sum _{j=1}^{N_h}N_j(\mathbf {x})\mathfrak {d}_{j} \end{aligned}$$

(38)

Analogous approximations are considered for test functions $\psi _{U_i}$ and $\psi _{\mathfrak {d}}$, whose nodal unknowns are indicated respectively as $\delta U_{ij}$ and $\delta \mathfrak {d}_j$

$$\begin{aligned} \psi _{U_i} (\mathbf {x})=\sum _{j=1}^{N_h}N_j(\mathbf {x})\delta U_{ij}, \quad \psi _\mathfrak {d}(\mathbf {x})=\sum _{j=1}^{N_h}N_j(\mathbf {x})\delta \mathfrak {d}_{j} \end{aligned}$$

(39)

In every single finite element e, one can define matrices $\mathbf {B}_U$ and $\mathbf {B}_{\mathfrak {d}}$ as

$$\begin{aligned} \mathbf {B}_U=\mathbf {D}_U\mathbf {N}_U, \quad \mathbf {B}_\mathfrak {d}=\mathbf {D}_\mathfrak {d}\mathbf {N}_\mathfrak {d} \end{aligned}$$

(40)

where, in a two dimensional setting, matrices $\mathbf {D}_U$ and $\mathbf {D}_{\mathfrak {d}}$ contain the derivatives with respect to coordinates $x_1$ and $x_2$

$$\begin{aligned} \mathbf {D}_U= & {} \left[ \begin{array}{c c} \frac{\partial }{\partial x_1} &{} 0 \\ 0 &{} \frac{\partial }{\partial x_2} \\ \frac{\partial }{\partial x_2} &{} \frac{\partial }{\partial x_1} \end{array} \right] ,\nonumber \\ \mathbf {D}_{\mathfrak {d}}= & {} \left[ \begin{array}{c} \frac{\partial }{\partial x_1}\\ \frac{\partial }{\partial x_2} \end{array} \right] , \end{aligned}$$

(41)

while $\mathbf {N}_U$ and $\mathbf {N}_{\mathfrak {d}}$ collect the shape functions

$$\begin{aligned} \mathbf {N}_U= & {} \left[ \begin{array}{c c c c c c c} N_1 &{} 0 &{} N_2 &{} 0 &{} ... &{} N_{Nnod} &{} 0\\ 0 &{} N_1 &{} 0 &{} N_2 &{} ... &{} 0 &{} N_{Nnod} \end{array} \right] , \nonumber \\ \mathbf {N}_{\mathfrak {d}}= & {} \left[ \begin{array}{c c c c} N_1&N_2&...&N_{Nnod} \end{array} \right] , \end{aligned}$$

(42)

being $N_{nod}$ the number of single element nodes. Thus, the weak form (37) can be written over each element domain $\varOmega _e$ as

$$\begin{aligned}&- {{\varvec{\delta }}} \mathbf {U}^T \int _{\varOmega _e} \mathbf {B}_U^T\, \mathbf {C}(\mathfrak {d}) \, \mathbf {B}_U \,\text {d}\varOmega \,\mathbf {U}+{{\varvec{\delta }}} \mathbf {U}^T \int _{\Gamma _{e_t}} \mathbf {N}_U^T\, \mathbf {t} \,\text {d}\Gamma \,\nonumber \\&+ {{\varvec{\delta }}} \mathbf {U}^T \int _{\varOmega _e} \mathbf {N}_U^T\, \mathbf {b} \,\text {d}\varOmega =0, \nonumber \\&{{\varvec{\delta }}} \mathfrak {d}^T G_C \ell \int _{\varOmega _e}\mathbf {B}_{\mathfrak {d}}^T\, \mathbf {B}_{\mathfrak {d}}\,\text {d}\varOmega \,\mathfrak {d} + {{\varvec{\delta }}} \mathfrak {d}^T \frac{G_C}{\ell }\int _{\varOmega _e}\mathbf {N}_\mathfrak {d}^T \mathbf {N}_{\mathfrak {d}}\,\text {d}\varOmega \,\mathfrak {d} \nonumber \\&+ \frac{1}{2}\mathbf {\delta }\mathfrak {d}^T \int _{\varOmega _e}\mathbf {N}_{\mathfrak {d}}^T \mathbf {U}^T\mathbf {B}_U^T\,\frac{\partial \mathbf {C}(\mathfrak {d})}{\partial \mathfrak {d}}\mathbf {B}_U\,\text {d}\varOmega \, \mathbf {U} = 0 \end{aligned}$$

(43)

The numerical solution of the coupled problem (43) is obtained by means of an iterative Newton-Raphson procedure, for which residual vectors of the displacement field $\mathbf {R_U}$ and of the phase field $\mathbf {R}_{\mathfrak {d}}$ are defined as

$$\begin{aligned}&\mathbf {R}_U =-\int _{\varOmega _e} \mathbf {B}_U^T\, \mathbf {C}(\mathfrak {d}) \, \mathbf {B}_U \,\text {d}\varOmega \,\mathbf {U} + \int _{\Gamma _{e_t}} \mathbf {N}_U^T\, \mathbf {t} \,\text {d}\Gamma \,\nonumber \\&+ \int _{\varOmega _e} \mathbf {N}_U^T\, \mathbf {b} \,\text {d}\varOmega , \end{aligned}$$

(44a)

$$\begin{aligned}&\mathbf {R}_\mathfrak {d} = -G_C\ell \int _{\varOmega _e}\mathbf {B}_{\mathfrak {d}}^T \mathbf {B}_{\mathfrak {d}}\,\text {d}\varOmega \,\mathfrak {d} - \frac{G_C}{\ell } \int _{\varOmega _e}\mathbf {N}_\mathfrak {d}^T \mathbf {N}_{\mathfrak {d}}\,\text {d}\varOmega \,\mathfrak {d}\nonumber \\&- \frac{1}{2} \int _{\varOmega _e}\mathbf {N}_{\mathfrak {d}}^T \mathbf {U}^T\mathbf {B}_U^T\,\frac{\partial \mathbf {C}(\mathfrak {d})}{\partial \mathfrak {d}}\mathbf {B}_U\,\text {d}\varOmega \, \mathbf {U} \end{aligned}$$

(44b)

The specific form of elemental stiffness matrices reads

$$\begin{aligned}&\mathbf {K}_{UU}^e=-\frac{\partial \mathbf {R}_U}{\partial \mathbf {U}}=\int _{\varOmega _e} \mathbf {B}_U^T\, \mathbf {C}(\mathfrak {d}) \, \mathbf {B}_U \,\text {d}\varOmega , \end{aligned}$$

(45a)

$$\begin{aligned}&\mathbf {K}_{U\mathfrak {d}}^e=-\frac{\partial \mathbf {R}_U}{\partial \mathfrak {d}}= \int _{\varOmega _e} \mathbf {B}_U^T\, \frac{\partial \mathbf {C}(\mathfrak {d})}{\partial \mathfrak {d}} \, \mathbf {B}_U\,\mathbf {U}\,\mathbf {N}_{\mathfrak {d}} \,\text {d}\varOmega , \end{aligned}$$

(45b)

$$\begin{aligned}&\mathbf {K}_{\mathfrak {d}U}^e=-\frac{\partial \mathbf {R}_{\mathfrak {d}}}{\partial U}= \int _{\varOmega _e} \mathbf {N}_{\mathfrak {d}}^T \,\mathbf {U}^T\,\mathbf {B}_U^T\,\frac{\partial \mathbf {C}(\partial \mathfrak {d})}{\partial \mathfrak {d}}\mathbf {B}_U\,\text {d}\varOmega , \end{aligned}$$

(45c)

$$\begin{aligned}&\mathbf {K}_{\mathfrak {d}\mathfrak {d}}^e=- \frac{\partial \mathbf {R}_{\mathfrak {d}}}{\partial \mathfrak {d}}=G_C\ell \int _{\varOmega _e} \mathbf {B}_{\mathfrak {d}}^T\,\mathbf {B}_{\mathfrak {d}}\,\text {d}\varOmega \nonumber \\&+\frac{G_C}{\ell }\, \int _{\varOmega _e} \mathbf {N}_{\mathfrak {d}}^T\, \mathbf {N}_{\mathfrak {d}}\,\text {d}\varOmega \nonumber \\&+ \frac{1}{2} \int _{\varOmega _e} \mathbf {N}_{\mathfrak {d}}^T \, \mathbf {U}^T\mathbf {B}_U^T \frac{\partial ^2 \mathbf {C}(\mathfrak {d})}{\partial \mathfrak {d}^2} \mathbf {B}_U \mathbf {U}\,\mathbf {N_{\mathfrak {d}}}\,\text {d}\varOmega \end{aligned}$$

(45d)

Consistently with the linearization of the resulting nonlinear system of Eq. (43), at each iteration of the Newton-Raphson loop, the following linear system has to be solved

$$\begin{aligned} \left[ \begin{array}{c c} \mathbf {K}_{UU} &{} \mathbf {K}_{U\mathfrak {d}} \\ \mathbf {K}_{\mathfrak {d} U} &{} \mathbf {K}_{\mathfrak {d}\mathfrak {d}} \end{array} \right] \left[ \begin{array}{c} \varDelta \mathbf {U}\\ \varDelta \mathfrak {d} \end{array} \right] = \left[ \begin{array}{c} \mathbf {R}_U\\ \mathbf {R}_{\mathfrak {d}} \end{array} \right] \end{aligned}$$

(46)

where $\varDelta \mathbf {U}$ and $\varDelta \mathfrak {d}$ are discretized according to (38) and the elemental stiffness matrices (45) and the residual vectors (44) have been assembled in the corresponding global ones.¹

Coupling with asymptotic homogenization is established by taking the closed-form expression for $\mathbf {C}(\mathfrak {d})$ and $\partial \mathbf {C}(\mathfrak {d})/\partial \mathfrak {d}$ provided by an off-line computation based on asymptotic homogenization, for different values of $\mathfrak {d}$.