A phase field approach for damage propagation in periodic microstructured materials

TL;DR

This paper introduces a multiscale computational method combining homogenization and phase field fracture modeling to analyze damage evolution in periodic composite materials, accounting for microstructural effects.

Contribution

It presents a novel finite element-based multiscale approach integrating homogenization with phase field fracture modeling for damage analysis in microstructured materials.

Findings

Damage evolution depends on microstructural topology and inclusion content.

The method predicts tensile strength and post-peak behavior influenced by microstructure.

Full microscopic field reconstruction is possible through down-scaling relations.

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 methodology is developed by means of the original combination of homogenization methods with the phase field approach of fracture. 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 allows considering 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 topology. 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 the inclusion volumetric content and its shape. Down-scaling relations allow the full reconstruction of the microscopic fields at any point of the macroscopic model, as a simple post-processing operation.

Equations

$$\frac{D}{D{x}_j}\left(C_{ijkl}^{(m,\epsilon )}\frac{D{u}_k}{D{x}_l}\right)+b_i=0,$$

$$C_{ijkl}^{(m,\epsilon )}\left(\mathbf{x}+{\mathbf{v}}_{\alpha }\right)=C_{ijkl}^{(m,\epsilon )}\left(\mathbf{x}\right),\alpha =1,2,\forall \mathbf{x}\in \mathcal{A}$$

$$C_{ijkl}^{(m,\epsilon )}\left(\mathbf{x}\right)=C_{ijkl}^m\left(\mathbf{\xi }=\frac{\mathbf{x}}{\epsilon }\right),\forall \mathbf{\xi }\in \mathcal{Q}.$$

$$u_k\left(\mathbf{x},\frac{\mathbf{x}}{\epsilon }\right)=\sum _{l=0}^{+\infty }{\epsilon }^l{u}_k^{(l)}\left(\mathbf{x},\frac{\mathbf{x}}{\epsilon }\right)=u_k^{(0)}\left(\mathbf{x},\frac{\mathbf{x}}{\epsilon }\right)+\epsilon u_k^{(1)}\left(\mathbf{x},\frac{\mathbf{x}}{\epsilon }\right)+{\epsilon }^2{u}_k^{(2)}\left(\mathbf{x},\frac{\mathbf{x}}{\epsilon }\right)+\mathcal{O}({\epsilon }^3)$$

$$\frac{D}{D{x}_j}f\left(\mathbf{x},\mathbf{\xi }=\frac{\mathbf{x}}{\epsilon }\right)={\left(\frac{\partial f}{\partial x_j}+\frac{1}{\epsilon }\frac{\partial f}{\partial {\xi }_j}\right)|}_{\mathbf{\xi }=\frac{\mathbf{x}}{\epsilon }}={\left(\frac{\partial f}{\partial x_j}+\frac{1}{\epsilon }{f}_{,j}\right)|}_{\mathbf{\xi }=\frac{\mathbf{x}}{\epsilon }}$$

$$\{{\epsilon }^{-2}{\left(C_{ijkl}^m{u}_{k,l}^{(0)}\right)}_{,j}+{\epsilon }^{-1}\left\{{\left[C_{ijkl}^m\left(\frac{\partial u_k^{(0)}}{\partial x_l}+u_{k,l}^{(1)}\right)\right]}_{,j}+\frac{\partial }{\partial x_j}\left(C_{ijkl}^m{u}_{k,l}^{(0)}\right)\right\}+$$

$$+{\left[C_{ijkl}^m\left(\frac{\partial u_k^{(1)}}{\partial x_l}+u_{k,l}^{(2)}\right)\right]}_{,j}+\frac{\partial }{\partial x_j}\left[C_{ijkl}^m\left(\frac{\partial u_k^{(0)}}{\partial x_l}+u_{k,l}^{(1)}\right)\right]+\epsilon \{{\left[C_{ijkl}^m\left(\frac{\partial u_k^{(2)}}{\partial x_l}+u_{k,l}^{(3)}\right)\right]}_{,j}+$$

$${+\frac{\partial }{\partial x_j}\left[C_{ijkl}^m\left(\frac{\partial u_k^{(1)}}{\partial x_l}+u_{k,l}^{(2)}\right)\right]\}+O({\epsilon }^2)\}|}_{\mathbf{\xi }=\frac{\mathbf{x}}{\epsilon }}+b_i(\mathbf{x})=0.$$

$$u_k^{(0)}\left(\mathbf{x},\mathbf{\xi }\right)=U_k\left(\mathbf{x}\right)$$

$$u_k^{(1)}\left(\mathbf{x},\mathbf{\xi }\right)=N_{kp{q}_1}^{(1)}(\mathbf{\xi })\frac{\partial U_p\left(\mathbf{x}\right)}{\partial x_{q_1}}$$

$$u_k^{(2)}\left(\mathbf{x},\mathbf{\xi }\right)=N_{kp{q}_1{q}_2}^{(2)}(\mathbf{\xi })\frac{{\partial }^2{U}_p\left(\mathbf{x}\right)}{\partial x_{q_1}\partial x_{q_2}}$$

$$⟨\left(\cdot \right)⟩=\frac{1}{|\mathcal{Q}|}{\int }_{\mathcal{Q}}\left(\cdot \right)d\mathbf{\xi }=0$$

$${\left(C_{ijkl}^m{N}_{kp{q}_1,l}^{(1)}\right)}_{,j}+C_{ijp{q}_1,j}^m=0$$

$${\left(C_{ijkl}^m{N}_{kp{q}_1{q}_2,l}^{(2)}\right)}_{,j}+\frac{1}{2}[{\left(C_{ijk{q}_2}^m{N}_{kp{q}_1}^{(1)}\right)}_{,j}+C_{i{q}_{1}p{q}_2}^m$$

$$+C_{i{q}_{2}kl}^m{N}_{kp{q}_1,l}^{(1)}+{\left(C_{ijk{q}_1}^m{N}_{kp{q}_2}^{(1)}\right)}_{,j}+C_{i{q}_{2}p{q}_1}^m+C_{i{q}_{1}kl}^m{N}_{kp{q}_2,l}^{(1)}]=$$

$$+\frac{1}{2}⟨C_{i{q}_{1}p{q}_2}^m+C_{i{q}_{2}kl}^m{N}_{kp{q}_1,l}^{(1)}+C_{i{q}_{2}p{q}_1}^m+C_{i{q}_{1}kl}^m{N}_{kp{q}_2,l}^{(1)}⟩$$

$$u_k(\mathbf{x},\mathbf{\xi })=U_k(\mathbf{x})+\epsilon N_{kp{q}_1}^{(1)}(\mathbf{\xi })\frac{\partial U_p(\mathbf{x})}{\partial x_{q_1}}+{\epsilon }^2{N}_{kp{q}_1{q}_2}^{(2)}(\mathbf{\xi })\frac{{\partial }^2{U}_p(\mathbf{x})}{\partial x_{q_1}\partial x_{q_2}}+\mathcal{O}({\epsilon }^3)$$

$$U_k(\mathbf{x})\doteq {⟨u_k\left(\mathbf{x},\frac{\mathbf{x}}{\epsilon }+\mathbf{\zeta }\right)⟩}_{\mathbf{\zeta }}$$

$${⟨g(\mathbf{\xi }+\mathbf{\zeta })⟩}_{\mathbf{\zeta }}=\frac{1}{|\mathcal{Q}|}{\int }_{\mathcal{Q}}g(\mathbf{\xi }+\mathbf{\zeta })d\mathbf{\zeta }=\frac{1}{|\mathcal{Q}|}{\int }_{\mathcal{Q}}g(\mathbf{\xi }+\mathbf{\zeta })d\mathbf{\xi }=⟨g(\mathbf{\xi })⟩$$

$${\bar{E}}_m=\frac{1}{|\mathcal{Q}|}{\int }_{\mathcal{L}}{\int }_{\mathcal{Q}}{\phi }_m\left(\mathbf{x},\mathbf{\xi }\right)d\mathbf{\xi }d\mathbf{x}={\int }_{\mathcal{L}}⟨{\phi }_m⟩d\mathbf{x}$$

$${\phi }_m=\frac{1}{2}{C}_{ijhk}^m\frac{D{u}_i}{D{x}_j}\frac{D{u}_h}{D{x}_k}=\frac{1}{2}{C}_{ijhk}^m\left({\delta }_{ip}{\delta }_{j{q}_1}+N_{ip{q}_1,j}^{(1)}\right)\left({\delta }_{hs}{\delta }_{k{r}_1}+N_{hs{r}_1,k}^{(1)}\right)\frac{\partial U_p}{\partial x_{q_1}}\frac{\partial U_s}{\partial x_{r_1}}+\mathcal{O}(\epsilon )$$

$$E_M={\int }_{\mathcal{L}}{\phi }_{M}d\mathbf{x}$$

$${\phi }_M=\frac{1}{2}{C}_{p{q}_{1}i{q}_2}\frac{\partial U_p}{\partial x_{q_1}}\frac{\partial U_i}{\partial x_{q_2}}.$$

$${\bar{E}}_m^0\doteq E_M$$

$$C_{p{q}_{1}i{q}_2}=⟨C_{rjkl}^m\left(N_{ri{q}_2,j}^{(1)}+{\delta }_{ir}{\delta }_{j{q}_2}\right)\left(N_{kp{q}_1,l}^{(1)}+{\delta }_{pk}{\delta }_{l{q}_1}\right)⟩.$$

$$\Pi \left(U_i,\mathfrak{d}\right)={\int }_{\Omega }\phi \left(H_{ij},\mathfrak{d}\right)\text{d}\Omega +{\int }_{\Omega }{\mathcal{G}}_c\gamma \left(\mathfrak{d},\frac{\partial \mathfrak{d}}{\partial x_j}\right)\text{d}\Omega -{\int }_{{\Gamma }_t}{t}_i{U}_i\text{d}\Gamma -{\int }_{\Omega }{b}_i{U}_i\text{d}\Omega .$$

$$H_{ij}=\frac{\partial U_i}{\partial x_j}.$$

$$\gamma \left(\mathfrak{d},\frac{\partial \mathfrak{d}}{\partial x_j}\right)=\frac{1}{2\ell }{\mathfrak{d}}^2+\frac{\ell }{2}{\left(\frac{\partial \mathfrak{d}}{\partial x_j}\right)}^2$$

$$\delta \Pi (U_i,\delta U_i,\mathfrak{d},\delta \mathfrak{d})$$

$$\delta \Pi (U_i,\delta U_i,\mathfrak{d},\delta \mathfrak{d})={\int }_{{\Gamma }_t}\frac{\partial U_i}{\partial x_j}{C}_{ijhk}{n}_k\delta U_h\text{d}\Gamma -{\int }_{\Omega }\frac{\partial }{\partial x_k}\left(\frac{\partial U_i}{\partial x_j}{C}_{ijhk}\right)\delta U_h\text{d}\Omega +$$

$${\int }_{\Omega }\frac{1}{2}{H}_{ij}\frac{\partial C_{ijhk}(\mathfrak{d})}{\partial \mathfrak{d}}{H}_{hk}\delta \mathfrak{d}\text{d}\Omega +{\int }_{\Omega }\frac{G_C}{\ell }\mathfrak{d}\delta \mathfrak{d}\text{d}\Omega +{\int }_{\Gamma }{G}_C\ell \frac{\partial \mathfrak{d}}{\partial x_j}{n}_j\delta \mathfrak{d}\text{d}\Gamma +$$

$$-{\int }_{\Omega }{G}_C\ell \frac{{\partial }^2\mathfrak{d}}{\partial x_j^2}\delta \mathfrak{d}\text{d}\Omega -{\int }_{{\Gamma }_t}{t}_i\delta U_{i}d\Gamma -{\int }_{\Omega }{b}_i\delta U_i\text{d}\Omega$$

$$\frac{\partial }{\partial x_j}\left(C_{ijhk}(\mathfrak{d})\frac{\partial U_h}{\partial x_k}\right)+b_i=0,$$

Full text

1 Introduction

The increasing demand in achieving lightweight structures with superior performance in terms of damage tolerance and load-bearing capacities has motivated the incorporation of composite materials in many different applications in the last decades. As a result, at present, composite-like structures can be found in biomechanics, renewable energy systems, aerospace and automotive sectors, to quote some practical fields with strong societal impact. The irruption of new manufacturing techniques has encouraged this trend with renovated interest, allowing the production of tailor-made composite materials with unprecedent capabilities by means of controlling the microstructural arrangement to attain the specific practical needs.

In this context, it is well established that the versatility and applicability of composite materials strongly depend upon of their inherent failure mechanisms that take place at different scales. These failure events can be generally categorized as meso-and-macro-scale mechanisms (Maimí et al., 2007; Reinoso et al., 2017a) and micro-scale failure modes (Herráez et al., 2015; Arteiro et al., 2015), with a cumbersome effective connection across the scales. From a mechanical point of view, the heterogenous character of composite materials induces their inherent anisotropic response at the structural level, being this behavior hardly predictable and thus representing one of the major drawbacks for their practical use. In order to overcome current limitations, the potential routes for achievement of superior responses in next generation of composite materials would require a thorough understanding of the influence of the micro-scale arrangements on macroscopic responses, in terms of constituents properties and their spatial distribution within the material, and the reliable prediction of failure events at such materials.

Investigation of composite materials having periodic or quasi-periodic arrangements resorts micromechanical approaches, which describe the material behavior in detail, but in general result in labor-intensive analyses. Multiscale techniques, based on homogenization approaches, could conveniently overcome these drawbacks, allowing to gather a synthetic but accurate description of the heterogeneous material behavior both for static and dynamic problems. Homogenization methods take into account the effects of the microscopic phases on the overall constitutive properties of such composites, also in the presence of multi-field phenomena. Specifically, homogenization methods allow replacing an heterogeneous material with a homogeneous equivalent one that can be modeled through either a Cauchy or a nonlocal continuum. Homogenization techniques have been a matter of intensive research within the last decades. In general sense, the local and/or nonlocal homogenization methods can be classified with respect to the underlying fundamental hypotheses as: the asymptotic techniques (Bensoussan et al., 1978; Bakhvalov and Panasenko, 1984; Gambin and Kröner, 1989; Hubert and Palencia, 1992; Allaire, 1992; Meguid and Kalamkarov, 1994; Boutin, 1996; Andrianov et al., 2008; Panasenko, 2009; Tran et al., 2012; Bacigalupo, 2014; Fantoni et al., 2017, 2018), the variational-asymptotic techniques (Willis, 1981; Smyshlyaev and Cherednichenko, 2000; Smyshlyaev, 2009; Bacigalupo and Gambarotta, 2014a, b; Bacigalupo et al., 2014; Del Toro et al., 2019) and many identification approaches, involving the analytical (Sevostianov et al., 2005; Bigoni and Drugan, 2007; Bacca et al., 2013a, b, c; Bacigalupo and Gambarotta, 2013; Sevostianov and Giraud, 2013; Rizzi et al., 2019a, b), and the computational techniques (Forest and Sab, 1998; Ostoja-Starzewski et al., 1999; Kouznetsova et al., 2002; Forest, 2002; Feyel, 2003; Kouznetsova et al., 2004; Kaczmarczyk et al., 2008; Yuan et al., 2008; Bacigalupo and Gambarotta, 2010; De Bellis and Addessi, 2011; Forest and Trinh, 2011; Addessi et al., 2013; Trovalusci et al., 2015; Otero et al., 2018; Dirrenberger et al., 2019).

Regarding the prediction of failure events in engineering materials and structures, the advent of new computational capabilities has promoted the generation of different numerical tools including diffusive crack methods (Bazant and Jirasek, 2002; Comi, 1999; Peerlings et al., 2001; Dimitrijevic and Hackl, 2011), strong discontinuity procedures (Moes et al., 1999; Linder and Armero, 2007; Oliver et al., 2006) and cohesive-like crack approaches (Camacho and Ortiz, 1996; Ortiz and Pandolfi, 1999; Paggi and Wriggers, 2012; Turon et al., 2018), among many others, where most of them rely on the exploitation of finite element(FE)-based procedures. Recent variational formulations and crack tracking algorithms based on the analogy between linear elastic fracture mechanics and standard dissipative systems can be found in (Salvadori and Fantoni, 2016; Salvadori et al., 2019). Derived from its versatility for the estimation of failure mechanisms due to crack initiation and growth, the seminal variational approach of fracture developed by Francfort and Marigo (1998), being denominated as the phase field approach of fracture, endows a smeared crack idealization that permits overcoming most of the limitations of alternative numerical methods. This methodology attains a regularized modeling of Griffith-like fracture (Griffith, 1921) and can be conceived as a nonlocal damage method, which in its original format is especially suitable for triggering fracture in brittle materials. In this concern, Bourdin et al. (2008, 2000) comprehensively developed the corresponding numerical treatment and extended the applicability of the phase field method for different applications (see also Pham and Marigo (2013); Pham et al. (2011)), providing a robust implementation. Alternatively, Miehe et al. (2010a) came up with a rigorous analysis with regard to the thermodynamic considerations within a engineering point of view. One of the main ingredients for the outstanding progresses on phase field methods stems from the fact that this variational approach does offer very appealing aspects and can easily be implemented into multi-field finite element frameworks. In view of the strong potential of the phase field methods, recent developments encompassed its application to cohesive-like fracture (Verhoosel and de Borst, 2013), coupled damage-plasticity (Ambati et al., 2015; Miehe et al., 2015a, 2016), shells (Miehe et al., 2014; Areias et al., 2016; Reinoso et al., 2017b), thermo-elastic (Miehe et al., 2015b) and hydrogen embrittlement (Martínez-Pañeda et al., 2018) applications, defining alternative degradation functions (Wu, 2017; Sargado et al., 2017), among many others. Owing to its modular formulation, the phase-field approach to fracture has proven to be a powerful tool for fracture characterization of many different materials such as arterial walls (Gültekin et al., 2018), and anisotropic media (Teichtmeister et al., 2017; Bleyer and Alessi, 2018; Quintanas-Corominas et al., 2019), to quote a few of them. Additional contributions in the field have addressed the accuracy of phase field methods evaluating the reliability of different implementation schemes with special attention on the local irreversibility constraint (Linse et al., 2017). Recently, the conjecture of crack nucleation in brittle materials discussed in (Tanné et al., 2018) has provided a further insight into the physical interpretation to the length scale in phase field formulations $\ell$, leading to a direct link of this variable with respect the material strength $\sigma_{c}$, see (Nguyen et al., 2016a).

With regard to heterogeneous media, being of special interest for composites, new developments of phase field methods generally accounted for the combination of bulk fracture with cracking along the existing interfaces. Specifically, Nguyen et al. (2016c) employed a level set function to describe interface displacement jump using an unique damage variable for both fracturing events, i.e. bulk and interface cracking. Conversely, Khisamitov and Meschke (2018) proposed a split of the dissipative term in the variational formalism distinguishing between bulk and interface fracture mechanisms, whereas an alternative formulation was developed in (Hansen-Dörr et al., 2019). The distinction between bulk and interface energy dissipation was seminally exploited by the authors through the coupling of the phase field and cohesive-like methods for triggering bulk and interface cracking, respectively, being denominated as the PF-CZM (Paggi and Reinoso, 2017). The robustness of this methodology has been assessed by its successful application to poly-crystalline materials (Paggi et al., 2018), layered ceramics (Carollo et al., 2018) and micro-mechanics of fiber-reinforced composites (Guillén-Hernández et al., 2019) in comparison with the so-called coupled criterion (Mantič and García, 2012).

In this respect, though there exist different numerical methods that allow reliable characterization of complex fracture phenomena at macro- and micro-scales, to the best authors’ knowledge, the coupling between both scales of observations for the efficient material tailoring has not been addressed so far with the exploitation of the phase field approach of fracture. Consequently, the primary goal of the current study is to establish a novel numerical methodology which links a variational phase field approach of fracture at the macro-scale with a homogenization-based techniques in order to account for the micro-structural information and the spread of damage depending on the material microstructure geometry and properties. This framework represents an interesting development towards the efficient multi-scale modeling of damage events in reinforced composite through the combination of the phase field method of fracture and a careful microstructure treatment by means of asymptotic homogenization techniques.

The manuscript is arranged as follows. Section 2 describes the developed methodology. The treatment of the corresponding material modeling at both scales is outlined in Section 3. The description of the damage evolution at the macro-scale is presented in Section 4. The proposed methodology is examined in Section 5 for different benchmark applications, whereas the main conclusions are given in Section 6.

2 Synopsis of the proposed approach

The present section highlights the key steps of the methodology herein proposed for multi-scale simulation of damage within heterogeneous materials by originally combining asymptotic homogenization and the phase field (nonlocal) approach to damage, as sketched in Fig. 1. In this respect, let consider a microstructured periodic medium, which is made by different phases and it is characterized by the macroscopic length scale $L$. At the microscale, each periodic cell is characterized by a microscopic length scale $\varepsilon$ which, for the validity of scale separation, has to be much smaller than $L$, i.e. $L\gg\varepsilon$.

Such a medium can be properly described trough a multiscale homogenization technique, which allows overcoming the computation cost of explicitly discretizing the heterogeneous microstructure using the finite element method. Thus, through the advocation of such a homogenization method, a concise but accurate description of material behavior and its mechanical damage can be gained. In this study, although the formulation could also be extended to other degradation mechanisms, we restrict our analysis to a situation where damage takes place only in the matrix, while the inclusions/fibers are assumed to behave obeying a linear elastic response with no degradation. This is a reasonable assumption in the majority of engineering applications concerning linear elastic stiff particle or fiber-reinforced composite materials. Hence, damage at the microscale is herein modeled as the degradation of the elastic properties of the matrix.

A two-scale asymptotic homogenization procedure is therefore advocated in order to homogenize the microstructured composite medium for any damage level, whose phases are modeled as first-order elastic continua, to an equivalent first-order elastic continuum.

At the macroscale, specimens of any geometry and subjected to any loading conditions can be considered by employing the finite element method and the phase field approach to (nonlocal) damage propagation. Usually, the damaged constitutive tensor entering the stress-strain relation is degraded in the phase field approach according to a prescribed degradation function $g(\mathfrak{d})$ which does not take into account the damage mechanism occurring at the microscale. Following a standard procedure with the context of Damage Mechanics, depending on the phase field damage variable $\mathfrak{d}$, the elastic parameters of the matrix will be degraded by the degradation function $g(\mathfrak{d})=\left(1-\mathfrak{d}\right)^{2}$, while that of the particle/fiber will be left undamaged. This will establish the coupling with asymptotic homogenization, which is used to compute the degraded constitutive tensor based on the above assumption on the damage mechanism. As a result, the overall homogenized damaged constitutive tensor will not be simply equal to the undamaged one re-scaled by $g(\mathfrak{d})$ as in the standard phase field approach, but it will depend on the damage mechanisms at the microscale and eventually on the microscale topology (inclusion shape and volumetric content) and linear elastic parameters of the material constituents.

In this computational framework, from the practical standpoint, the asymptotic homogenization scheme is be called at each Gauss point of the finite element discretization at the macroscale, in order to compute the homogenized constitutive tensor of the damaged material, based on the damage-like phase field variable $\mathfrak{d}$. This resembles the FE2 method, where two nested finite element models are solved, one at the macroscale, and another at the microscale. To reduce the computation cost associated with the present formulation, we herein explore a look-up table concept. Basically, the current asymptotic homogenization is applied off-line to a series of unit periodic cell problems for a series of values of $\mathfrak{d}$ ranging from zero to unity, i.e intact and fully degraded states, respectively. This information is subsequently employed to degrade the elastic properties of the matrix. The obtained homogenized constitutive tensor components, for different values of $\mathfrak{d}$, are then interpolated to determine their closed-form (analytical) dependency upon $\mathfrak{d}$, that to be used by the phase field approach for the computation of the damaged elastic strain energy density of the damaged material. Similarly, the first derivative of the components of the homogenized damaged constitutive tensor with respect to $\mathfrak{d}$ are also analytically determined, to be inserted into the computation of the tangent constitutive operator at the macroscale, which is required for fully consistent implicit iterative-incremental solution schemes.

Finally, as a post-process of the simulations, the homogenized stress and strain fields can be visualized. Moreover, by exploiting the down-scaling relations established by the asymptotic homogenization approach, the microscale fields can also be reconstructed and inspected, for any material point. Algorithm 1 outlines the complete numerical-analytical approach herein developed for a given pseudo time increment throughout the simulation.

In summary, the main goal of the proposed formulation is to study the evolution of damage inside the bulk depending on the microstructural features of the material itself. In particular, by stating that damage takes place only in the matrix, the proposed methodology allows investigating the role of the material microstructure (inclusion shape and volumetric content) onto the macroscopic response, both in terms of load-displacement relation and in terms of damage propagation. Moreover, it is worth to state that the proposed methodology precludes the use of computational demanding FE2 schemes, and whose use can be widely exploited in different applications.

A range of values for the characteristic phase field length scale $\ell$ are also examined, under the assumption that diffuse damage is spread enough in the space to avoid strain localization and therefore guarantee the applicability of homogenization theory.

As shown in the results obtained from the numerical simulations, the proposed coupled model is able to predict an apparent strength of the system that is not only affected by $\ell$ (which is usually linked to the apparent material strength), as in the standard phase field approaches for homogenous materials that do not account for different degradation scenarios at the microscale, but also on the microstructural topology and the inclusion volumetric content.

3 Periodic elastic material modeled at two scales

To establish a general scheme, one considers a linear elastic material characterized by a periodic microstructure, whose phases can be described as first-order continua (Fig. 2). In a two-dimensional perspective, vector $\mathbf{x}=x_{1}\mathbf{e}_{1}+x_{2}\mathbf{e}_{2}$ identifies the position of each material point in the orthogonal coordinate system with origin at point $O$ and base $\mathbf{e}_{1},\mathbf{e}_{2}$ as shown in Fig. 2-(b).

Advocating the periodicity of the medium microstructure, a periodic cell $\mathcal{A}=[0,\varepsilon]\times[0,\delta\varepsilon]$ is identified by two orthogonal periodicity vectors $\mathbf{v}_{1}=\varepsilon\mathbf{e}_{1}$ and $\mathbf{v}_{2}=\delta\varepsilon\,\mathbf{e}_{2}$, where $\varepsilon$ characterizes the size of cell $\mathcal{A}$ as depicted in Fig. 2. As is usually done in the context of asymptotic homogenization, the unit cell $\mathcal{Q}=[0,1]\times[0,\delta]$ is through obtained rescaling the periodic cell $\mathcal{A}$ by the characteristic length $\varepsilon$. The microscopic displacement $\mathbf{u}(\mathbf{x})=u_{i}(\mathbf{x})\mathbf{e}_{i}$ is induced by body forces $\mathbf{b}(\mathbf{x})=b_{i}(\mathbf{x})\mathbf{e}_{i}$ whose periodicity is much greater than the microstructural one for the validity of the scale separation condition. Body forces in fact, are here assumed to be $\mathcal{L}$-periodic, and to have vanishing mean values on $\mathcal{L}$, being $\mathcal{L}=[0,L]\times[0,\delta L]$ a true representative portion of the whole body. Consistently with what above said, the structural length $L$ has to be much greater than the microstructural one (${L}\gg\varepsilon$). Moreover, diffuse damage without strain localization phenomena is also required to allow the rigorous applicability of homogenization theory.

Under the hypothesis of quasi-static loading processes, the partial differential equation governing the elastic problem is given componentwise as

[TABLE]

being $C_{ijkl}^{(m,\varepsilon)}$ the components of the fourth order micro elasticity tensor, where the superscript $m$ refers to the microscale and $\varepsilon$ to the characteristic size of $\mathcal{A}$. In Eq. (1), the symbol $D/Dx_{j}$ denotes the generalized derivative with respect to the variable $x_{j}$. The micro-constitutive elastic tensor is $\mathcal{Q}$-periodic, meaning that

[TABLE]

and its components depend only on the fast variable ${\bm{\xi}}=\mathbf{x}/\varepsilon$:

[TABLE]

Two variables in fact can be identified for the separation of the macro and the micro scales, namely the macroscopic (slow) one $\mathbf{x}\in\mathcal{A}$ and the microscopic (fast) one ${\bm{\xi}}\in\mathcal{Q}$. In this regard, the micro displacement results to be a function of both the slow and the fast variables, i.e. $u_{i}=u_{i}(\mathbf{x},{\bm{\xi}}=\frac{\mathbf{x}}{\varepsilon})$. In such a context, homogenization techniques reveal to be efficient in accurately describing the overall behavior of a microstructured composite material, replacing the periodic continuum with an homogeneous equivalent one requiring a much lower computation cost for its analysis and simulation. Analytical and numerical solutions of Eq. (1) could be particularly complex to obtain because of the rapidly oscillating $\mathcal{Q}$-periodic coefficients. In what follows, the formulation of an equivalent first-order elastic continuum is described, together with the exact closed-form expression of the overall constitutive elastic tensor components.

In the obtained equivalent homogenized medium, the macro displacement field in each material point $\mathbf{x}$ is denoted with $\mathbf{U}(\mathbf{x})=U_{i}(\mathbf{x})\mathbf{e}_{i}$. It is worth mentioning that variability of source terms acting on the medium has to be much greater than the microstructural length scale $\varepsilon$, in order to preserve the scales separability principle. If volume forces are $\mathcal{L}$-periodic, then macroscopic displacement $\mathbf{U}(\mathbf{x})$ will be $\mathcal{L}$-periodic too, otherwise suitable boundary conditions should be taken into account in order to determine the macroscopic field (Fantoni et al., 2017, 2018).

3.1 Asymptotic expansion of the microscopic field equation

Splitting the contributions of the slow and fast variables in accordance with the scale separation condition, the micro displacement field can be asymptotically expanded in terms of the micro characteristic length scale $\varepsilon$ (Bakhvalov and Panasenko, 1984) as

[TABLE]

Taking into account the differentiation rule for a generic function $f(\mathbf{x},{\bm{\xi}}=\frac{\mathbf{x}}{\varepsilon})$:

[TABLE]

and substituting the asymptotic expansion (4) into the governing field equation (1), one obtains

[TABLE]

A set of recursive partial differential problems can now be obtained from Eq.(3.1) by equating terms at the different orders of $\varepsilon$. Details about the solutions of such recursive differential problems can be found in (Smyshlyaev and Cherednichenko, 2000; Bacigalupo, 2014; Fantoni et al., 2017). In particular, the solvability condition of the partial differential problem at the order $\varepsilon^{-2}$ in the class of $\mathcal{Q}$-periodic functions (Bakhvalov and Panasenko, 1984), implies that the solution at the order $\varepsilon^{-2}$ does not depend on the fast variable ${\bm{\xi}}$, taking the form

[TABLE]

In force of solution (7) of the elastic differential problem at the order $\varepsilon^{-2}$, the solvability condition of resulting differential problem at the order $\varepsilon^{-1}$ and the $\mathcal{Q}$-periodicity of $\mathbb{C}^{m}$, lead to a solution at the first order of the form

[TABLE]

where $N_{kpq_{1}}^{(1)}$ is the first order perturbation function, which accounts for the influence of microstructural inhomogeneities.

Accounting for the solutions (7) and (8) of differential problems at the order $\varepsilon^{-2}$ and $\varepsilon^{-1}$, respectively, and of the $\mathcal{Q}$-periodicity of the micro constitutive tensor components $C_{ijkl}^{m}$ and of the perturbation function $N_{kpq_{1}}^{(1)}$, at the order $\varepsilon^{2}$ the micro displacement solution takes the form

[TABLE]

where $N_{kpq_{1}q_{2}}^{(2)}$ is the second order perturbation function which, again, depends only on the geometrical and physico-mechanical properties of the considered material.

Determination of perturbation functions derives from the solution of non-homogeneous recursive differential problems at the different order $\varepsilon$, being those problems referred to as cell problems. Such differential problems are characterized by $\mathcal{Q}$-periodic source terms having a vanishing mean value over unit cell $\mathcal{Q}$. Consequently, their solutions result to be $\mathcal{Q}$-periodic and they are enforced to have a zero mean value over $\mathcal{Q}$ in order to guarantee their uniqueness.

The following normalization condition is therefore imposed for all perturbation functions:

[TABLE]

where $|\mathcal{Q}|$ represents the unit cell area, namely $|\mathcal{Q}|=\delta$. The cell problem at the order $\varepsilon^{-1}$ in terms of the first order perturbation function $N_{kpq_{1},l}^{(1)}$ reads

[TABLE]

Once $N_{kpq_{1}}^{(1)}$ is determined, from Eq. (3.1) one obtains the cell problem at the order $\varepsilon^{0}$ which, symmetrized with respect to indices $q_{1}$ and $q_{2}$, takes the form

[TABLE]

Determination of perturbation functions as solutions of cell problems at the different orders of $\varepsilon$ allows obtaining the down-scaling relation, expressing the micro displacement field $\mathbf{u}(\mathbf{x},{\bm{\xi}})$ as an asymptotic expansion of powers of $\varepsilon$ in terms of the $\mathcal{L}$-periodic macro field $\mathbf{U}(\mathbf{x})$ and its gradients. The down-scaling relation is expressed as

[TABLE]

Denoting with ${\bm{\zeta}}\in\mathcal{Q}$ a variable such that the vector $\varepsilon{\bm{\zeta}}\in\mathcal{A}$ defines the translation of the medium with respect to $\mathcal{L}$-periodic body forces $\mathbf{b}(\mathbf{x})$ (Smyshlyaev and Cherednichenko, 2000; Bacigalupo, 2014), the macro displacement field can in turn be expressed in terms of the micro field trough the up-scaling relation, as the mean value of $\mathbf{u}$ over unit cell $\mathcal{Q}$

[TABLE]

In particular, the variable ${\bm{\zeta}}$ removes rapid fluctuations of coefficients and is such that invariance property

[TABLE]

is proved to hold for $\mathcal{Q}$-periodic functions. Eq. (15), together with normalization condition (10) enforced for all perturbation functions, leads to the up-scaling relation (14).

3.2 Closed form of the homogenized elastic tensor

The overall elastic tensor is determined by means of a generalized macro-homogeneity condition, which establishes an energetic equivalence between the macro and the micro scales (Smyshlyaev and Cherednichenko, 2000; Bacigalupo, 2014).

The microscopic mean strain energy $\bar{E}_{m}$ is defined in terms of the micro strain energy density $\varphi_{m}$ as

[TABLE]

where $\varphi_{m}$ reads

[TABLE]

in view of the down-scaling relation (13). The strain energy at the macroscale is defined as

[TABLE]

in terms of the macro strain energy density $\varphi_{M}$ expressed in the form

[TABLE]

In Eq. (19), $C_{pq_{1}iq_{2}}$ denotes a generic component of the overall elastic tensor. Denoting with $\bar{E}_{m}^{0}$ the mean micro strain energy truncated at the zeroth order, the generalized macro-homogeneity condition implies that

[TABLE]

from which the components of the first order macroscopic tensor $C_{pq_{1}iq_{2}}$ can be defined in terms of the components of the microscopic constitutive elastic tensor and in terms of the first order perturbation function. The closed form of the overall elastic tensor components is

[TABLE]

Better approximation of the solution of the heterogeneous problem could be obtained by means of higher order homogenization techniques, modifyng the generalized macro-homogeneity condition (20) in a proper way. Higher order approaches consider both higher order terms in the asymptotic expansion performed at the microscale and the non local character of the macroscopic constitutive tensor. Alternatively, the solution of the average field equations of infinite order by means of perturbation methods allows to obtain a higher order approximation by suitably truncating the asymptotic expansion of the macro field $\mathbf{U}$ in powers of $\varepsilon$ (Bacigalupo, 2014). Higher order homogenization schemes could be explored in follow-up studies, motivating the development of higher order phase field formulations at the macroscale, for a consistent modelling of higher order effects across the material length scales.

4 Damage evolution in a microstructured material trough a phase field model approach

Let consider an arbitrary, homogeneous body $\Omega\in\mathbb{R}^{2}$. Note that for the sake of simplicity, we restrict our analysis to bidimensional geometries. As described in Sec. 3, vector $\mathbf{x}=x_{1}\mathbf{e}_{1}+x_{2}\mathbf{e}_{2}$ identifies the position of each material point inside the bulk at the macroscale. The boundary of the body $\Gamma\in\mathbb{R}$ is the union of Dirichlet $\Gamma_{U}$ and a Neumann $\Gamma_{t}$ parts, with $\Gamma_{U}\cup\Gamma_{t}=\Gamma$ and $\Gamma_{U}\cap\Gamma_{t}=\emptyset$.

Under the validity of the small strains assumption, the material response to the following quasi-static external actions is sought: body forces $\mathbf{b}(\mathbf{x})$ in $\Omega$, tractions $t$ on $\Gamma_{t}$, and displacements $\bar{U}$ on $\Gamma_{U}$. Accordingly to the variational approach of fracture described in (Bourdin et al., 2008; Miehe et al., 2010b), the total potential energy of a system $\Pi$ can be written as the balance between the internal contribution, which is the sum of a global energy storage functional and a dissipation functional due to damage, and the contribution due to external loading:

[TABLE]

Note that the first term on the right-hand side of Eq. (22) is the global energy storage functional in which the strain energy density per unit volume $\varphi(H_{ij},\mathfrak{d})$ depends on both the macro displacement field trough $H_{ij}$ and the phase field $\mathfrak{d}$. In the small-strain context, $H_{ij}$ is defined as the displacement gradient:

[TABLE]

According to (Miehe et al., 2010b, c), the phase field $\mathfrak{d}\in[0,1]$ is an internal auxiliary variable characterizing, for $\mathfrak{d}=0$, the undamaged condition and, for $\mathfrak{d}=1$, the fully damaged condition of the material.

The second term in Eq. (22) defines the work needed to create a diffusive damage and its rate provides the dissipation power due to damage in the bulk. In this context, $\mathcal{G}_{C}$ would correspond to a Griffith-type fracture energy, and $\gamma(\mathfrak{d},\partial\mathfrak{d}/\partial x_{j})$ represents the crack surface density function depending on the phase field and on its gradient, which leads to the nonlocality of the formulation and therefore its mesh independency. It can be expressed as

[TABLE]

where $\ell$ is an internal length scale parameter governing the regularization/diffusion of damage (Kuhn and Müller, 2014; Nguyen et al., 2016b). Here, we request to $\ell$ to be big enough to avoid crack localization and model situations where diffuse damage takes place. In order to introduce damage only in tension, the energy density in the bulk is left undamaged in compression. On the other hand, it is degraded in tension by assuming that only the Young’s modulus of the matrix is affected by damage, i.e., by setting $E_{m}(\mathfrak{d})[(1-\mathfrak{d})^{2}+\mathcal{K}]E_{m,0}$, where $E_{m,0}$ is the undamaged value of the Young’s modulus of the matrix and $\mathcal{K}$ is a residual material stiffness introduced to avoid numerical instabilities for $\mathfrak{d}=1$. The overall effective constitutive tensor components $C_{ijhk}(\mathfrak{d})$ is finally computed at the microscale by asymptotic homogenization and it is provided to the finite element computation scheme at the macroscale. In this regard, it is important to highlight that the elastic strain energy density is not degraded as in the standard phase field approach used in homogenous materials. Therefore, since the damaged constitutive tensor is the outcome of asymptotic homogenization, it will also depend upon the volumetric content and the shape of the inclusion/fiber.

The first order variation of functional (22) reads:

[TABLE]

Integration by parts of the first and fourth terms of the right-hand side of Eq. (25) allows to rephrase the variation $\delta\Pi$ as

[TABLE]

where symmetry properties of the elastic tensor $\mathbb{C}$ and the condition $\delta U_{i}=0$ on $\Gamma_{U}$ have been exploited.

The Euler-Lagrange equations associated with the displacement and phase field problems yield the coupled field equations

[TABLE]

in the domain $\Omega$, along with Neumann boundary conditions:

[TABLE]

The solution of the nonlinear system (4), together with boundary conditions (4), allows to describe the evolution of macroscopic mechanical quantities and of the phase field $\mathfrak{d}(\mathbf{x})$ inside the homogenized material.

5 Numerical examples: influence of the phase field internal length scale and microstructure upon damage evolution

The main objective of the present section is to show the dependence of the evolution of damage inside the material to the shape and volumetric content of the inclusion, as well as to the phase field internal length $\ell$. To this aim, let consider the specimens sketched in Fig. 3, having width $L=1$ mm and height $2L$, both subjected to uniform tensile loading through the application of an imposed far-field displacement $\Delta$ to all the finite element nodes belonging to the lower and upper sides. The case depicted in Fig. 3-(b) distinguishes from the one of Fig. 3-(a) because of the presence of an initial defect, which is a horizontal straight edge crack with mouth at $x_{2}=L$ and crack tip at $x_{1}=L/2$. In both cases, the specimen is supposed to be constituted by a microstructured material composed by a matrix embedding a circular or a square inclusion. Damage is assumed to develop only within the matrix.

In order to accelerate coupling between the macroscale and the microscale computations, asymptotic homogenization of the microstructured material as described in Section 3 is performed off-line, for a set of admissible values of the phase field variable $\mathfrak{d}$ ranging from zero to unity. The numerical results of the coupled problem (4) for the test problem without initial defect and with an initial defect are discussed in Sections 5.2 and 5.3, respectively.

5.1 Homogenization of the material with microstructure

Specimens depicted in Fig. 3 are supposed to be made by a material with microstructure whose periodic cell is shown in Fig. 4. The periodic cell is composed of an Aluminum matrix undergoing damage, with an undamaged linear elastic constitutive behavior characterized by a Young’s modulus $E_{m,0}=E_{Al}=60$ GPa and Poisson ratio $\nu_{Al}=0.3$, and a linear elastic inclusion corresponding to Silicon carbide with Young’s modulus $E_{SiC}=340$ GPa and $\nu_{SiC}=0.18$.

Under plane strain conditions, the elastic tensor of the undamaged matrix $\mathbb{C}_{Al}$ and of the inclusion $\mathbb{C}_{SiC}$ are:

[TABLE]

Two different inclusion shapes have been taken into account, namely circular (Fig. 4-(a)) with five different values of the volumetric content $f=1/4,1/8,1/16,1/32,1/100$, and square (Fig. 4-(b)) with a volumetric content $f=1/4$. The volumetric content represents the ratio between the area of the inclusion and that of the periodic cell $\mathcal{A}$ in the cross-section of the specimen, as customary.

After computing the perturbation functions $N_{kpq_{1}}^{(1)}$ as the solutions of the cell problem (11), the overall elastic constitutive tensor has been computed by means of the closed form (21). Considering, for example, a volumetric content $f=1/4$, we get:

[TABLE]

where subscripts circ and sq stand for circular and square topology, respectively.

Adopting a residual stiffness $\mathcal{K}=0.005$ in correspondence of $\mathfrak{d}=1$, the elastic constitutive tensor of the damaged matrix $\mathbb{C}_{Al}$ have been multiplied by degradation function $g(\mathfrak{d})=(1-\mathfrak{d})^{2}+\mathcal{K}$, with the phase field variable in the range $0\leq\mathfrak{d}\leq 1$. The components of the overall elastic tensor $\mathbb{C}$ have then been computed for each value of $\mathfrak{d}$. Fig. 5 shows, for example, the components of $\mathbb{C}_{circ}$ (Fig. 5-(a)) and of $\mathbb{C}_{sq}$ (Fig. 5-(b)) as a function of the phase field $\mathfrak{d}$, for a microstructure configuration characterized by a volumetric content $f=1/4$. In correspondence of $\mathfrak{d}=0$, the overall elastic constitutive tensors $\mathbb{C}_{circ}$ and $\mathbb{C}_{sq}$ correspond to those of the undamaged material, given by Eq. (30). Their values monotonically decrease as $\mathfrak{d}$ approaches unity, where the value $\mathcal{K}$ is retrieved.

Considering a constant value for $G_{C}=6$ N/mm for the homogenized material, the coupled system of equations (4) has been numerically solved for each value of the imposed displacement $\Delta$, for the unnotched and the notched specimen configurations. Details regarding the finite element framework adopted to obtain the numerical solution, and implemented in the finite element software FEAP (Zienkiewicz and Taylor, 1977), are described in the Appendix A. Five different values of the internal length scale $\ell$ have been considered, in order to investigate the role of this parameter in triggering damage inside the microstructured material. To avoid strain localization which would violate the applicability of homogenization, values of $\ell$ have been set intermediate between the value of the macroscopic length scale $L$ and the value of the microscopic one $\varepsilon$, namely $\ell=0.05,0.1,0.2,0.4$ mm.

5.2 Case 1: tensile test on a plain unnotched specimen

One considers the plain unnotched specimen depicted in Fig. 3-(a), where a uniaxial tensile test is performed. Fig. 6 depicts the variation of the homogenized mean stress $\bar{T}_{22}$ with respect to the homogenized mean strain $\bar{H}_{22}$, by varying the phase field length scale $\ell$, for a circular inclusion with a different volumetric content $f$ in the subfigures from (a) to (e), and for a the square inclusion with $f=1/4$ in the subfigure (f). In particular, the average stress $\bar{T}_{22}$ is computed as the sum of the reaction forces on the upper (or lower) boundary of the specimen, divided by the width $L$ of the specimen and its unit out-of-plane thickness. The average strain $\bar{H}_{22}$ is computed by dividing the imposed far-field displacement $\Delta$ by $L$.

According to previous results reported in the literature (Kuhn and Müller, 2014; Nguyen et al., 2016b; Paggi and Reinoso, 2017), the apparent strength $\bar{T}_{22_{max}}$, representing the maximum of the average stress-strain curves, increases as the length scale $\ell$ decreases. This trend can be intuitively explained by the fact that larger values of $\ell$ indicate a more widespread damage in the bulk. However, by comparing the different subfigures in Fig. 6, one notices that the present model provides an apparent strength also dependent upon the volumetric content of the reinforcement and on the shape of the inclusion, which is a major novelty as compared to the standard phase field approach for homogeneous materials. Examining each stress-strain curve more in detail, the system response is almost linear till $\bar{T}_{22_{max}}$ is reached. The stiffness of such an elastic part is the same by varying $\ell$, which is consistent with the fact that all the curves in a subfigure present the same microstructure with the same volumetric content $f$. For an axial strain larger than the value of the strain corresponding to the apparent strength, the post-peak response presents softening. This is a major difference from the prediction of the phase field in the case of a homogenous material, and it is caused by the progressive stress transfer from the damaged matrix to the undamaged inclusion, by increasing damage. This redistribution of stresses cannot take place if the material is homogenous. The value of the residual stress are increasing by reducing the length scale $\ell$, as an effect of the different evolution and spread of damage. By comparing Figs. 6-(a) and 6-(f), which refer to the same value of $f$, but correspond to different inclusion shapes, one can notice how the apparent strength and the corresponding value of $\bar{H}_{22_{max}}$ are slightly influenced by the change of topology, while the post peak behavior is greatly affected by the inclusion shape, showing steeper branches and lower residual stresses in the case of a square inclusion. Furthermore, Fig. 6-(f) shows that the residual stresses are slightly influenced by $\ell$, for a square topology of the cell inclusion.

To highlight the effect of the volumetric content, Fig. 7 depicts, for a circular inclusion, the average stress $\bar{T}_{22}$ as a function of $\bar{H}_{22}$ by keeping $\ell$ constant and varying $f$. The apparent strength $\bar{T}_{22_{max}}$ increases as volume fraction $f$ increases, and this is coherent with the choice to weaken the constitutive properties only for the matrix and not for the inclusion. For the same reason, residual stresses are greater as volume fraction $f$ increases.

To summarize the above results, Fig. 8 shows how the apparent strength of the system, $\bar{T}_{22_{max}}$, and the correspondent values of the strain, $\bar{H}_{22_{max}}$, depend on $\ell$, for the different inclusion volumetric contents and also for different shape of the inclusion. Curves referring to the square topology and $f=1/4$ show lower values than those for a circular inclusion and the same value of $f$, both for $\bar{T}_{22_{max}}$ and $\bar{H}_{22_{max}}$.

5.3 Case 2: tensile test for a specimen with an edge notch

Let consider a uniaxial tensile test for the notched specimen depicted in Fig. 3-(b). Fig. 9 shows the average stress $\bar{T}_{22}$ as a function of average strain $\bar{H}_{22}$ for this test, where $\bar{T}_{22}$ and $\bar{H}_{22}$ are computed as for the specimen without an initial notch. In particular, Fig. 9 shows the homogenized stress-strain curves for a material microstructure with a circular inclusion and different values of $f$, for two distinct values of the internal length scale $\ell$, namely $\ell=0.1$ mm (Fig. 9-(a)) and $\ell=0.4$ mm (Fig. 9-(b)). Once again, the two-scale coupled formulation predicts an apparent strength $\bar{T}_{22_{max}}$ which is an increasing function of $f$, while the correspondent average deformation $\bar{H}_{22_{max}}$ is a decreasing function of $f$. The mechanical response is quite affected by $\ell$ as can be assessed by comparing Figs. 9-(a) and 9-(b). Softening branches occur for $\ell=0.4$ mm, while for $\ell=0.1$ mm the post-peak response is quite brittle. Fig. 10 shows $\bar{T}_{22_{max}}$ and the corresponding strain $\bar{H}_{22_{max}}$ as functions of $\ell$, for all the different values of $f$ herein considered. The apparent strength of the system and the corresponding homogenized strain are monotonically decreasing functions of the internal length scale $\ell$. Moreover, the comparison between Figs. 8 and 10 shows that, as expected, $\bar{T}_{22_{max}}$ and $\bar{H}_{22_{max}}$ are lower than the corresponding values for the unnotched specimen.

Figures 11 and 12 show two contour plots of the phase field variable $\mathfrak{d}$ for $\ell=0.1$ mm and $\ell=0.4$ mm, respectively, at two consecutive pseudo-time steps of the quasi-static computation in the case of a circular inclusion and a volumetric content $f=1/4$. In particular Fig. 11-(a) refers to $\bar{H}_{22}=0.0062$, Fig. 11-(b) to $\bar{H}_{22}=0.0072$ , Fig. 12-(a) to $\bar{H}_{22}=0.014$, and Fig. 12-(b) to $\bar{H}_{22}=0.0172$, respectively. As expected, the spread of the damaged zone increases by increasing the internal length $\ell$.

Without any loss of generality, the cases shown in Fig. 12 have been considered as possible examples to show the applicability of the down-scaling relations established by asymptotic homogenization to reconstruct the local microscopic fields. In fact, once the macro displacement field $\mathbf{U}(\mathbf{x})$ is known in each point $\mathbf{x}$ of the macroscale model, it is possible to reconstruct the micro displacement field $\mathbf{u}\left(\mathbf{x},{\bm{\xi}}\right)$ on the periodic cell by means of the down-scaling relation (13), herein truncated to the first order of $\varepsilon$. For example, considering a point of coordinates $\mathbf{x}=\{0.75,0.12\}$ mm inside the specimen sketched in Fig. 3-(b), the computed dimensionless micro displacement $\tilde{u}_{1}(\mathbf{x},{\bm{\xi}})=u_{1}(\mathbf{x},{\bm{\xi}})/U_{1}(\mathbf{x})$ is shown in Fig. 13 for $\ell=0.4$ mm and a circular inclusion with $f=1/4$, for two different values of $\bar{H}_{22}$. In particular, Fig. 13-(a) refers to an imposed average strain equal to $\bar{H}_{22}=0.014$, and a macro displacement $U_{1}(\mathbf{x})=-2.0113\,10^{-3}$ mm, while Fig. 13-(b) refers to $\bar{H}_{22}=0.0172$ and $U_{1}(\mathbf{x})=-2.2932\,10^{-3}$ mm.

The dimensionless micro displacement $\tilde{u}_{2}(\mathbf{x},{\bm{\xi}})=u_{2}(\mathbf{x},{\bm{\xi}})/U_{2}(\mathbf{x})$ in the $\mathbf{e}_{2}$ direction, at the same point $\mathbf{x}=\{0.75,0.12\}$ mm, is shown in Fig. 14. Fig. 14-(a) refers to an imposed strain equal to $\bar{H}_{22}=0.014$ and $U_{2}(\mathbf{x})=-1.1741\,10^{-2}$ mm, while Fig.14-(b) refers to $\bar{H}_{22}=0.0172$ and $U_{2}(\mathbf{x})=-1.4217\,10^{-2}$ mm.

6 Conclusions

The evolution of damage in composite materials with periodic or quasi-periodic microstructure has been investigated in the present work trough a multiscale finite element-based computational approach, integrating the phase field method with an homogenization scheme. Specifically, the overall constitutive properties of the heterogeneous material have been determined in closed form by means of a two-scale asymptotic homogenization technique and the damage propagation inside the equivalent homogeneous medium has been described at the macroscale via a phase field approach.

The current technique allowed triggering macroscopic damage-like by means of tracking the matrix degradation at the microscale. In this regard, differing from alternative FE2-techniques, this method provided a fast and efficient computational tool for the link between the two-scales using the concept of a look-up table scheme. Therefore, it was precluded the instantaneous computation of both scales throughout the simulation.

The performance of the present technique has been demonstrated by means of several representative examples. These results clearly evidenced that, for cases under analysis, the overall response of the specimen was affected by the material length scale, which is usually related to the apparent material strength of the material, and by the micromechanical information regarding the inclusion shape and volume content as well as the damage extent into the matrix. Moreover, through these numerical applications, it has been shown that these aspects played a significant role in the the peak response of the specimen and the posterior softening evolution upon failure, featuring a progressive softening instead of an abrupt drop, as in most of the phase field methods of fracture.

Finally, it is worth to recall that the influence of the microscale over the macroscale results is one of the major novelties of the present approach in comparison to standard phase field methods for homogeneous materials with no specific scale separation.

Acknowledgements

AB would like to acknowledge the financial support by National Group of Mathematical Physics (GNFM-INdAM). MP would like to acknowledge the financial support of the Italian Ministry of Education, University and Research to the Research Project of National Interest (PRIN 2017). The authors would like to thank the IMT School for Advanced Studies Lucca for its support to the stays of FF and JR in the IMT Campus as visiting researchers in 2019, making possible the realization of this joint work.

Appendix A 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 (4) 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 (4) of the coupled field problem detailed in Section 4, taking into account boundary conditions (4), reads

[TABLE]

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}$

[TABLE]

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}$

[TABLE]

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

[TABLE]

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}$

[TABLE]

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

[TABLE]

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

[TABLE]

The numerical solution of the coupled problem (A) 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

[TABLE]

The specific form of elemental stiffness matrices reads

[TABLE]

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

[TABLE]

where $\Delta\mathbf{U}$ and $\Delta\mathfrak{d}$ are discretized according to (32) and the elemental stiffness matrices (39) and the residual vectors (38) have been assembled in the corresponding global ones222Finally note that following (Miehe et al., 2010a), the current formulation is equipped with a viscous crack resistance parameter, leading to the modification of the operators associated with the phase field variable. An alternative solution scheme would encompass a staggered Jacobi-type method which can be easily recalled by eliminating the coupling stiffness matrices and adopting an alternate minimization procedure.

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}$.