# A computational homogenization framework with enhanced localization criterion for macroscopic cohesive failure in heterogeneous materials

Lu Ke, Frans P van Der Meer

<figure><img src="image_1.png" /></figure>

## To cite this version:

Lu Ke, Frans P van Der Meer. A computational homogenization framework with enhanced localization criterion for macroscopic cohesive failure in heterogeneous materials. Journal of Theoretical, Computational and Applied Mechanics, 2022, Volume 2, pp.1-28. 10.46298/jtcam.7707. hal-03298979v3

HAL Id: hal-03298979

https://hal.science/hal-03298979v3

Submitted on 4 Mar 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.7707

OAI hal-03298979v3

**History**

Received July 27, 2021

Accepted Feb 9, 2022

Published Feb 28, 2022

**Associate Editor**

Laurence Brassart

**Reviewers**

Thierry Massart

Phu Nguyen

Alfredo Huespe

**Open Review**

OAI hal-03589655

**Supplementary Material**

Data permalink
DOI 10.5281/zenodo.6044934

**Licence**

CC BY 4.0

©The Authors

# A computational homogenization framework with enhanced localization criterion for macroscopic cohesive failure in heterogeneous materials

Lu KE and Frans P. VAN DER MEER

Faculty of Civil Engineering and Geosciences, Delft University of Technology, The Netherlands

Computational homogenization allows to let the macroscopic constitutive behavior of materials emerge from microscale simulations without loss of generality with respect to microstructure and microscale constitutive response. Although computationally demanding, computational homogenization works very well for the hardening response of materials where the macroscopic stress and strain fields are smooth. However, in case of softening materials, when localization of deformation takes place, special care is needed to ensure objectivity of the method. In this paper, a generic multiscale computational homogenization approach for modeling onset and propagation of cracks in heterogeneous materials that is capable of considering various microscale mechanisms is presented. The common acoustic tensor bifurcation criterion is reinforced by an additional condition to help detect the localization mode more robustly. After the onset of macroscale localization, a key scale transition parameter is needed to translate the macroscopic displacement jump to an averaged strain over the micromodel domain. Then the macroscale crack is governed by a homogenized traction-separation relation evaluated from the underlying micromodel in which micro-failure accumulates. The scale transition parameter is studied for a range of different scenarios and endowed with a geometrical interpretation. Various numerical tests have been performed to confirm the objectivity and validity of the framework. The framework is generic in the sense that no assumptions on the microscale constitutive or kinematic representation of material failure are made in the scale transition. The framework is also highly compatible with the first order computational homogenization, which minimizes the additional complexity of adding macroscopic crack growth to the computational implementation.

**Keywords:** multiscale modeling, computational homogenization, localization, phantom node method, cohesive zone model

## 1 Introduction

Synthetic composite materials such as concrete and fibre-reinforce polymers are widely used in a variety of engineering sectors due to their superior mechanical performance and/or excellent durability compared to traditional single-component materials. From design perspective, this type of heterogeneous materials provides much flexibility to meet requirements for different applications. However, exploiting this design freedom requires ability to predict the homogenized response of heterogeneous materials accurately. Complicated phenomenological constitutive models are often adopted to describe the overall nonlinear material behavior with apparent properties on the macroscale. Apart from the troublesome calibration process with a number of internal variables, the heterogeneities and direct interactions among constituents on the fine scale are difficult to be taken into account by such kind models, which limits the accuracy as well as the ability to explain underlying physical mechanisms. As an alternative approach for analysis of heterogeneous materials, the multiscale method (Budarapu et al. 2019) explicitly connects the mechanical behavior and failure processes across different scales of observation.

Among a wide range of multiscale approaches (Matouš et al. 2017; Fish 2006; Miehe and Bayreuther 2007; Tadmor et al. 1996; Hill 1965; van der Sluis et al. 2000), the current contribution is particularly confined to a two-scale homogenization-based approach which is computational
homogenization (CH), also referred to as FE² (Feyel and Chaboche 2000), that provides a nested solution of two boundary value problems (BVPs) defined on each of scales. The point-wise overall response (usually stress) on the macroscale is directly evaluated, on the fly, from the lower scale problem which accounts for heterogeneous description of the material morphology and properties. Two fundamental assumptions, namely *Separation of scales* and *Periodicity of microstructures*, together with the *Hill-Mandel macrohomogeneity principle* (Hill 1965) serve as the basis of the classical first-order computational homogenization (FOCH) approach (Kouznetsova et al. 2001) for macroscopic problems within linear and hardening regimes. However, after entering the softening regime when localization happens, the two assumptions become questionable, which leads to two severe numerical consequences: Firstly, the macroscopic BVP is ill-posed (Nguyen et al. 2011), which requires the regularization of the macroscopic model. Secondly, results are not objective with respect to the size of the representative volume element (RVE) (Gitman et al. 2007). The second of these is related to the fact that the existence of a softening RVE is challenged as material loses statistical homogeneity upon strain localization.

To overcome above deficiencies of FOCH for strain localization cases, various models has been developed in recent years. For cases with weak discontinuities (Rudnicki and Rice 1975) (discontinuous strain field) such as shear bands, a second-order computational homogenization (SOCH) scheme has been proposed by Kouznetsova et al. (2002). The main improvement resides in the scale transition which not only includes the macroscopic deformation gradient but also its Lagrangian gradient. This approach is very useful for homogenization of shells and beams (Helfen and Diebels 2014; Coenen et al. 2010). However, SOCH cannot properly deal with softening materials exhibiting deformation beyond a quadratic nature in the displacements. Towards strong discontinuities (Simo et al. 1993) (discontinuous displacement field), several special homogenization schemes based on specifically selective RVE-area averaging rules have been developed to obtain homogenized cohesive laws for macroscale material points in which the localized deformation is represented. The first approach in this category was developed in the context of masonry cracking (Massart et al. 2007). It essentially relies on the proper incorporation of a localization band at the macroscale, along with a two-fold FOCH to resolve the response in both damaging and unloading material parts. On the macroscale, cracks are discretized using the embedded discontinuity approach (Linder and Armero 2007) while on the microscale, localization bands appear as a consequence of softening material laws. Gitman et al. (2008) proposed an alternative solution, the coupled-volume approach, establishing a rigorous coupling between a macroscopic finite element size and the damaging RVE. However, it more resembles a domain decomposition approach with an embedded scale refinement. By adopting extended finite element method (XFEM) on both macro and microscale, Belytschko et al. (2008) proposed the multiscale aggregating discontinuities method of which the essential feature is the injection of an equivalent crack at the macroscopic model. Based on the introduced failure-zone averaging technique which guarantees the objectivity of softening RVEs, Nguyen et al. (2011) and Nguyen et al. (2012a) proposed a continuous-discontinuous CH framework to model the transition of microscopic diffusive damage to macroscopic cohesive failure for tensile cracking problems. Regarding the same type of two-scale problems, Sánchez et al. (2013) proposed a failure-oriented multiscale variational formulation in which a special macro-micro mechanism of kinematical coupling is defined for the softening RVEs. The subsequent additional boundary conditions (BC) applied over the localized domain for those RVEs lead to the objectivity of the formulation. This work has been further extended by Toro et al. (2016a) and Toro et al. (2016b) to cases where microscopic failure mode is a set of cohesive cracks forming a dominant failure path. By establishing a direct scale transition relation between the macroscopic and equivalent microscopic displacement jump terms, Coenen et al. (2012a) and Bosco et al. (2015) proposed a framework in which the macroscopic traction is included in the macroscopic problem as Lagrange multipliers. Similarly, based on the direct assumption of definition on macroscopic effective traction or displacement jump, the work by Turteltaub and Suárez-Millán (2020) established an energetically-consistent CH framework that upscales the microscopic cohesive traction to macroscale discontinuity. Instead of resorting to cohesive models, Khoei and Saadat (2019) proposed a CH framework which adopts non-local damage model on both scales for the purpose of regularization through establishing an additional transition relation for non-local damage terms. These CH models have been successful in dealing
with multiscale strain localization problems. However, they all have a certain assumption on the
particularly targeted failure pattern on microsale and most of them implement a homogenization
scheme based on the integration over a dominant localization zone or cracking surface in the
micro model. Subsequently, a dedicated identification technique of such averaging regions has to
be put forward, which is, for complex micro structure exhibiting various dissipation mechanisms,
not trivial. For the case where distributed localization bands exist in the softening RVE, it is even
harder to identify a single major percolation path (Turteltaub et al. 2018). Besides, the high
complexity of each model impedes the smooth incorporation of newly-developed parts into
existing FEM codes. Thus, an approach for strain localization problems with more generality that
imposes minimal restrictions on the micromodel formulation is worth further exploration.

Regarding the above pursuit, Hirschberger et al. (2009) has given some enlightenment
even though their paper focused on pre-existent interfaces and material layers rather than
arbitrary propagating cracks on both scales. The proposed macroscopic displacement-jump strain
tensor term, $h^{-1}(\|\mathbf{u}^M\| \otimes^{\text{sym}} \mathbf{n}_d^M)$ in which $\|\mathbf{u}^M\|$ is the displacement jump, $\mathbf{n}_d^M$ is the normal
vector to the discontinuity and $h$ is a scaling parameter, makes it possible that FOCH could
be extended to localization cases. The same strain tensor term has been applied to cases with
propagating cracks to define the softening micromodel boundary conditions in (Mercatoris and
Massart 2009a) in which the parameter $h$ is interpreted as the width of an embedded localization
band following (Massart et al. 2007). Oliver et al. (2015) used the same relation in a multiscale
framework based on the Continuum Strong Discontinuity approach on both scales. The most
distinctive feature of this work is the extension of classical FOCH towards softening materials
without substantial changes in the formulation. For material points within discontinuities, the
above strain term is added to the macroscopic bulk strain to formulate boundary conditions for a
microscopic BVP. Then a similar procedure as FOCH could be followed to obtain homogenized
cohesive laws. Basically, the same FEM code blocks for bulk material points could be reused
for localized points with minor changes, which makes the approach less intrusive than other
frameworks. Adopting the similar homogenization methodology, Svenning et al. (2017a) and
Svenning et al. (2019) made improvements by introducing the localization aligned weakly periodic
BCs to achieve better accuracy and convergence rate. Following the same path, the current
contribution presents a generic CH framework which is tailored to dealing with macroscopic
crack onset and propagation for heterogeneous materials without restrictive assumptions on the
microscopic failure processes and remains as less intrusive to the existing code as possible. The
major features are as follows:

* The key scale transition parameter, *h*, to transform the macroscopic displacement jump into a
strain term is studied for different micro-failure mechanisms, which is crucial for maintaining the
energetic equivalence, or extended *Hill-Mandel macrohomogeneity principle*, between scales.

* Regarding the localization criterion on the macroscale, a robust scheme based on the singularity
condition of the acoustic tensor is proposed to accurately detect the localization direction, which
is essential for capturing arbitrary crack growth on the macroscale.

* For describing failure, discontinuous approaches are implemented on both macroscale and
microscale. To be specific, phantom-node method (van der Meer et al. 2012) with cohesive cracks
is adopted on the macroscale while a model where interelement cohesive elements are inserted
on the fly (Camacho and Ortiz 1996) is used on the microscale for efficiently simulating scattered
micro-cracks. Generality with respect to micromodel formulation is illustrated by also including
plasticity in the micromodel.

## 2 A two-scale failure homogenization framework

In the present framework, we consider a two-scale failure problem. On the macroscale, regardless
of the heterogeneous nature, the bulk material in the intact solid domain is modelled as
homogeneous continuum with effective constitutive laws defined through the classic FOCH. As
the microscopic failure process gradually develops and damage accumulates, strain localization is
triggered on the macroscale. A discontinuity represented by a cohesive crack is inserted into the
macroscopic strain-localized region. The point-wise cohesive tractions are obtained through a
modified FOCH procedure with new scale-transition relations that relate macro-discontinuity
kinematics to the response of softening microscopic BVPs.

Note that due to the nature of the scale transition relations, any kind of microscopic dissipative processes can be incorporated into the framework.

In the following sections, two superscripts, namely □M and □m, are adopted to denote items that are defined on the macroscale and microscale, respectively.

## 2.1 Macroscopic cohesive failure problem

As schematically shown in Figure 1, the macroscopic problem can be described by a solid domain ΩM containing discontinuities represented by cohesive cracks ΓdM (lines in two-dimensions or surfaces in three-dimensions) to which the outward unit normal vector is denoted by n$_{d}$M, with the convention that n$_{d}$M = −n$_{d}$M is pointing to the positive side +ΓdM while +n$_{d}$M = −n$_{d}$M is pointing to the opposite side. Cohesive tractions t$_{d}$M = −t$_{d}$M are present across the crack and related to the

<figure><img src="image_3.png" /><figcaption>Figure 1 Macroscopic failure problem with cohesive cracks.</figcaption></figure>

displacement jump which is defined as [[u$^{M}$]] = +u$^{M}$ − −u$^{M}$ where u$^{M}$ is the displacement field. Prescribed tractions t̂ are imposed on Neumann boundary Γt̂M while prescribed displacements û are applied on Dirichlet boundary Γu$^{M}$.

The governing equations for quasi-static analysis in the absence of body forces are

$$ \nabla \cdot \boldsymbol{\sigma}^{\mathrm{M}} = \mathbf{0} \quad \text{in} \quad \Omega \setminus \Gamma_d^{\mathrm{M}} \qquad (1) $$

$$ \boldsymbol{\sigma}^{\mathrm{M}} \cdot \mathbf{n}^{\mathrm{M}} = \hat{\mathbf{t}} \quad \text{on} \quad \Gamma_t^{\mathrm{M}} \qquad (2) $$

$$ \mathbf{u}^{\mathrm{M}} = \hat{\mathbf{u}} \quad \text{on} \quad \Gamma_u^{\mathrm{M}} \qquad (3) $$

and

$$ +\mathbf{t}_d^M + -\mathbf{t}_d^M = \mathbf{0}; \quad \mathbf{n}_d^M \cdot \boldsymbol{\sigma}^M = -\mathbf{t}_d^M \quad \text{on} \quad \Gamma_d^M, \qquad (4) $$

where σM is the Cauchy stress tensor, ∇ is the gradient operator and nM is the normal vector to the boundary ∂ΩM. Across the crack, Cauchy theorem is exploited to link the stress σM to the cohesive traction t$_{d}$M. For simplicity, the dependence on position xM is omitted.

Under the small displacement-gradient assumption, the infinitesimal strain field at the bulk point reads

$$ \boldsymbol{\varepsilon}^{\mathrm{M}} = \nabla^{\mathrm{sym}} \mathbf{u}^{\mathrm{M}} = \frac{1}{2} (\nabla \mathbf{u}^{\mathrm{M}} + (\nabla \mathbf{u}^{\mathrm{M}})^{\mathrm{T}}) \quad \text{in} \quad \Omega^{\mathrm{M}} \setminus \Gamma_d^{\mathrm{M}}, \qquad (5) $$

where the superscript T denotes transposition and □sym refers to the symmetric part of the term. The weak formulation of the equilibrium can be interpreted as the balance of virtual work

$$ \int_{\Omega^{\mathrm{M}} \setminus \Gamma_d^{\mathrm{M}}} \nabla^{\mathrm{sym}} \delta \mathbf{u}^{\mathrm{M}} : \boldsymbol{\sigma}^{\mathrm{M}} d\Omega + \int_{\Gamma_d^{\mathrm{M}}} \delta [[\mathbf{u}^{\mathrm{M}}]] \cdot \mathbf{t}_d^{\mathrm{M}} d\Gamma = \int_{\Gamma_t^{\mathrm{M}}} \delta \mathbf{u}^{\mathrm{M}} \cdot \hat{\mathbf{t}} d\Gamma, \qquad (6) $$

where the left-hand-side is the internal contribution from the bulk material and cohesive cracks while the right-hand-side is the external virtual work. The problem is completed by constitutive laws for bulk and cohesive points which can be denoted as

$$ \boldsymbol{\sigma}^{\mathrm{M}} = \boldsymbol{\sigma}^{\mathrm{M}}(\boldsymbol{\varepsilon}(\mathbf{u}^{\mathrm{M}}), \kappa); \quad \mathbf{t}_d^{\mathrm{M}} = \mathbf{t}_d^{\mathrm{M}}([[\mathbf{u}^{\mathrm{M}}]], \boldsymbol{\varepsilon}(\mathbf{u}^{\mathrm{M}}), \omega), \qquad (7) $$
where $κ$ and $ω$ are internal variables. Equation (7) also includes possible dependency of the cohesive traction on the bulk strain. Although most cohesive formulations relate the traction only to the displacement jump, it can be beneficial to include dependency on the stress in the surrounding material (see e.g. van der Meer and Sluys 2009), and in the proposed multiscale formulation this coupling is also present. In the computational homogenization approach, material history is tracked on the microscale if history-dependent laws such as plasticity and damage are adopted. Thus, no internal variable is explicitly defined on the macroscale, which is exactly the advantage of CH approach since constitutive relations for the macroscopic problem that represent the complex evolution of the microscopic state are not needed.

## 2.2 Kinematics scale transition

In the current deformation-driven CH framework, the boundary conditions of the microscopic BVP follow from macroscopic deformations. Thus, kinematic relations are needed for the macro-to-micro scale transition.

### 2.2.1 Macroscopic deformation measures

Deformation measures are defined for two parts of the macroscale domain, namely the continuous part $\Omega_c^M \stackrel{\text{def}}{=} \Omega^M \setminus \Gamma_d^M$ and the discontinuous part $\Gamma_d^M$. For bulk points in the domain $\Omega_c^M$, the macrostrain field $\epsilon_c^M$ is described according to Equation (5) as

$$ \epsilon_c^M = \nabla^{\text{sym}} \mathbf{u}^M. \qquad (8) $$

For cohesive points on $\Gamma_d^M$, in order to include all discontinuity information in one single macroscopic quantity and allow for a scale transition that does not require any information on what is happening inside the micromodel, we define the cohesive macrostrain field following (Hirschberger et al. 2009) and (Oliver et al. 2015) as

$$ \epsilon_d^M = \epsilon_b^M + \epsilon_f^M \quad \text{with} \quad \epsilon_b^M \stackrel{\text{def}}{=} \frac{1}{2}(\epsilon_{\Gamma_d^M}^+ + \epsilon_{\Gamma_d^M}^-) \qquad (9) $$

and

$$ \epsilon_f^M \stackrel{\text{def}}{=} \frac{1}{h} (\|\mathbf{u}^M\| \otimes^{\text{sym}} \mathbf{n}_d^M), \qquad (10) $$

where $\epsilon_b^M$ is the bulk strain average from $+\Gamma_d^M$ and $-\Gamma_d^M$; $\epsilon_f^M$ is the second-order discontinuity tensor. The latter quantity is formulated by introducing a scaling parameter $h$ which plays an important role in the scale transition.

### 2.2.2 Microstrain average

On the microscale, the linear kinematics description from FOCH is adopted. Consequently, the micro displacement field is defined as

$$ \mathbf{u}^M = \mathbf{u}^M + \epsilon^M \cdot \mathbf{x}^M + \tilde{\mathbf{u}}^M, \qquad (11) $$

where $\epsilon^M$ is the macrostrain ($\epsilon_c^M$ or $\epsilon_d^M$) of the material point to which the RVE relates; $\tilde{\mathbf{u}}^M$ is the displacement fluctuation field inside the RVE domain. Due to the assumption of separation of scales, macroscopic quantities $\mathbf{u}^M$ and $\epsilon^M$ are constant inside the micromodel domain, i.e. independent of $\mathbf{x}^M$. Note that if fracture is present in the micromodel and a discrete cracking model (e.g. XFEM, cohesive elements) is used, $\tilde{\mathbf{u}}^M$ will be discontinuous.

The domain average of microstrain is defined as

$$ \langle \epsilon^M \rangle_{\Omega^M} \stackrel{\text{def}}{=} \frac{1}{|\Omega^M|} \int_{\Omega^M} \epsilon^M d\Omega = \frac{1}{|\Omega^M|} \int_{\Omega^M} \nabla^{\text{sym}} \mathbf{u}^M d\Omega, \qquad (12) $$

where $\langle \square \rangle_{\Omega^M}$ denotes the domain averaging operation; $|\Omega^M|$ is the measure of the domain $\Omega^M$ which is the area for two-dimensions and volume for three-dimensions.

Replacing Equation (11) into Equation (12) with application of divergence theorem leads to

$$ \langle \epsilon^M \rangle_{\Omega^M} = \epsilon^M + \frac{1}{|\Omega^M|} \int_{\partial \Omega^M} \tilde{\mathbf{u}}^M \otimes^{\text{sym}} \mathbf{n}^M d\Gamma, \qquad (13) $$
where $\mathbf{n}^m$ is the outward normal to the domain boundaries $\partial\Omega^m$. With appropriate boundary conditions, the second term on the RHS vanishes. In this work, strong periodic BCs are used. However, the present framework is independent micro-BCs as long as they fulfill the requirement of energetic equivalence (Hill-Mandel macrohomogeneity principle).

### 2.2.3 Definition of scale transition parameter $h$

From Equation (10), $h$ can be regarded as a numerical characteristic length in the present multiscale material model which is required to obtain objective and realistic results as indicated by Bažant (2010). $h$ serves to translate the displacement jump to an averaged strain over the RVE with finite size; a smearing operation where $[\mathbf{u}^M]$ is smeared over a band with width $h$ to obtain a strain. Note that in the actual micromodel the deformation is not actually smeared, but for defining BCs a smeared deformation measure is needed that includes the strain as well as the displacement jump. The value of $h$ is acquired from the underlying microscopic problem accounting for the developed failure mechanism and morphology of the microstructure. However, its definition varies in literature.

The simplest definition is given by Hirschberger et al. (2009) adopting the RVE height as $h$ in case of modeling a horizontal material layer. In that case, a straight material interface crosses the RVE from the left boundary to the right boundary. In the later work of Souza and Allen (2011), $h$ is defined as the ratio between the volume of the RVE and vectorial area of the external boundary projected onto macrocrack normal vector, which accounts for a slanted macrocrack. Note that both definitions do not consider any failure detail inside the RVE.

Assuming discrete micro-crack failure, Oliver et al. (2015) proposed the definition

$$h = |\Omega^m|/|\Gamma_d^m| \qquad (14)$$

where $|\Gamma_d^m|$ is the measure of open micro-cracks (length in two-dimension and area in three-dimension). Similarly, in case of a continuum model with smeared micro-damage, we could replace $|\Gamma_d^m|$ with the length of the dominant damage band.

The basic assumption in the previous definition is that the entire bulk material elastically unloads. However, if the plasticity is present on the microscale as for the case of composites delamination failure (Jalalvand et al. 2016; van der Meer 2016), then the microscopic dissipative process is a competition between plasticity and cohesive cracking. This could leads to a scenario where both mechanisms co-exist even after onset of localization. In this case, it is not obvious how to define $\Gamma_d^m$.

Another factor that also influences $\Gamma_d^m$ is the choice of the macroscopic localization criterion. As an essential element of the framework, such a criterion is adopted to detect the moment of introducing the macroscale discontinuity. However different criteria exist (Verhoosel et al. 2010; Neilsen and Schreyer 1993; Nguyen et al. 2012a). These criteria are satisfied along the simulation at different moments, which yields different times of localization. For instance, Massart et al. (2007) has shown that the loss of ellipticity might be fulfilled in the hardening region, i.e. before the limit point is reached in the load-deflection response in case of masonry material failure. Consequently, the microscale failure mechanism detected at the localization moment depends on the choice of the particular criterion. One scenario is that micro-cracks/damage zones are still dispersed rather than localizing into one single crack/localization band, as schematically indicated in Figure 2, which therefore leads to an estimate for the fracture length $|\Gamma_d^m|$ that is not representative for the post peak response.

Svenning et al. (2017b) defined $h$ based on the assumption of an effective discontinuity surface $\Gamma_{eff}^m$ as the plane through the center of RVE parallel to the macroscopic crack plane

$$h = |\Omega^m|/|\Gamma_{eff}^m| \qquad (15)$$

with, for a square RVE in 2D

$$|\Gamma_{eff}^m| = \begin{cases} \sqrt{|\Omega^m|}/\cos \alpha, & 0^\circ \le \alpha \le 45^\circ \\ \sqrt{|\Omega^m|}/\sin \alpha, & 45^\circ \le \alpha \le 90^\circ \end{cases} \qquad (16)$$
<figure><img src="image_4.png" /><figcaption>**Figure 2** Microscale failure mechanism at different load steps. Black lines in the micromodel show cohesive cracks that are damaging, while gray lines represent unloading cracks.</figcaption></figure>

where $\alpha$ is the angle between the effective discontinuity and the global x-axis. Note that for a circular type of RVE (Hofman et al. 2021), Equation (16) is not valid while Equation (15) still holds.

An incorrect choice for $h$ eventually affects the trans-scale energetic equivalence as discussed in Section 2.3.2. Therefore a study into the proper choice for $h$ for a range of different cases will be presented in the numerical examples section of this paper.

## 2.3 Micro-to-macro homogenization

After solving the microscale BVP, microscale quantities need to be upscaled to the macroscale in a homogenization step. The main requirement of this process is the energetic consistency condition which will be addressed later on.

### 2.3.1 Stress and traction homogenization

Regarding macroscale bulk points in the $\Omega_c^M$, the homogenized stress is acquired through the domain averaging as in classical FOCH as

$$ \sigma_c^M = \langle \sigma^m \rangle_{\Omega^m} = \frac{1}{|\Omega^m|} \int_{\Omega^m} \sigma^m d\Omega. \qquad (17) $$

With equilibrium at the microscale and divergence theorem, the volume integral on the RHS can be rewritten as a surface integral (Kouznetsova 2002)

$$ \frac{1}{|\Omega^m|} \int_{\Omega^m} \sigma^m d\Omega = \frac{1}{|\Omega^m|} \int_{\partial\Omega^m} t^m \otimes x^m d\Gamma, \qquad (18) $$

where $t^m = \sigma^m \cdot n^m$ is the traction vector over the RVE boundary.

For extracting cohesive tractions for points on the macroscale discontinuity $\Gamma_d^M$, combining Equation (17) with the Cauchy theorem gives

$$ t_d^M = \sigma_d^M \cdot n_d^M = \langle \sigma^m \rangle_{\Omega^m} \cdot n_d^M, \qquad (19) $$

where $\sigma_d^M$ is the homogenized stress from applying the macrostrain $\epsilon_d^M$ to the damaged RVE. Note that tractions are computed from integration over the whole RVE domain rather than just the failure zone.

Within the nonlinear finite element framework, the consistent tangent of the constitutive relation is crucial to maintain fast convergence. With the above homogenized stress and traction at hand, the tangents for bulk and cohesive integration points are given as

$$ \delta \sigma_c^M = D^M : \delta \epsilon_c^M \qquad (20) $$

$$ \delta t_d^M = T_{\epsilon}^M : \delta \epsilon_b^M + T_d^M : \delta [\mathbf{u}]^M, \qquad (21) $$

where $D^M$ is the bulk tangent stiffness while $T_{\epsilon}^M$ and $T_d^M$ are, respectively, cohesive tangent stiffness term related to the macroscopic bulk strain and displacement jump.

From Equation (10), Equation (19) and Equation (20), it follows

$$ \delta t_d^M = \frac{1}{h} \cdot n_d^M \cdot D^M : \delta \epsilon_b^M + \frac{1}{h} \cdot n_d^M \cdot D^M \cdot n_d^M \cdot \delta [\mathbf{u}]^M. \qquad (22) $$
Comparing Equation (21) and Equation (22), it is clear that

$$ T_{\varepsilon}^{M} = \frac{1}{h} \cdot n_{d}^{M} \cdot D^{M}; \qquad T_{d}^{M} = \frac{1}{h} \cdot n_{d}^{M} \cdot D^{M} \cdot n_{d}^{M}, \qquad (23) $$

where both cohesive tangent terms are explicitly related to bulk tangent stiffness $D^M$, scaling parameter $h$ and macro-crack normal vector $n_d^M$.

### 2.3.2 Energy consistency

As mentioned previously, the fundamental requirement for a valid two-scale CH is that the energy remains consistent across scales. Being fulfilled in FOCH, the Hill-Mandel macrohomogeneity principle states that the local rate of work on the macroscale is equal to the volume average of the rate of work performed on the underlying RVE, which can be expressed as

$$ \sigma^M : \dot{\varepsilon}^M = \langle \sigma^m : \dot{\varepsilon}^m \rangle_{\Omega^m}. \qquad (24) $$

For the bulk material part, with the homogenized stress, see Equation (17), the kinematics scale transition relation (13) and the chosen boundary conditions, Equation (24) can be proved valid, see for instance (Kouznetsova 2002). However, for the softening RVE, the standard Hill-Mandel condition is not directly applicable. When discrete cracking is the microscale failure mechanism, following (Turteltaub et al. 2018), the extended version of such energetic equivalence principle for localization points reads

$$ \sigma_d^M : \dot{\varepsilon}_d^M = \langle \sigma^m : \dot{\varepsilon}^m \rangle_{\Omega^m} + \frac{1}{|\Omega^m|} \int_{\Gamma_d^m} t_d^m \cdot [[\dot{\mathbf{u}}]]^m d\Gamma. \qquad (25) $$

Replacing Equation (9) into Equation (25) and using the Cauchy theorem yields

$$ \sigma_d^M : \dot{\varepsilon}_b^M + \frac{1}{h} \cdot t_d^M \cdot [[\dot{\mathbf{u}}]]^M = \langle \sigma^m : \dot{\varepsilon}^m \rangle_{\Omega^m} + \frac{1}{|\Omega^m|} \int_{\Gamma_d^m} t_d^m \cdot [[\dot{\mathbf{u}}]]^m d\Gamma, \qquad (26) $$

where $[[\dot{\mathbf{u}}]]^m$ is the rate of micro-crack jump. When continuum damage is the microscale failure mechanism, the second term on the RHS of Equation (26) drops out.

In both cases the scaling parameter $h$ plays a role in regulating the energetic equivalence between the scales in the fracture homogenization process. Thus, the correctness of definition of $h$ has an influence on the validity of the homogenization.

## 2.4 Macroscopic localization analysis

As discussed in Section 2.2, detection of the initiation and propagation of the macroscale cohesive crack is crucial for obtaining physically consistent results. For this purpose, one could adopt criteria based on the macroscopic stress or strain field (Verhoosel et al. 2010; Wells and Sluys 2001; Nguyen et al. 2012b; Svenning et al. 2017b). In such criteria, a stress or strain quantity computed from the homogenized macroscopic stress or strain at every integration point is compared to a failure threshold. Despite the fact that microstructural effects have influence on the macroscopic stress or strain field, a failure stress or strain threshold still has to be defined empirically on the macroscale, which violates the principle that the heterogeneous material response is described on the microscale.

Another type of approach is the discontinuous bifurcation analysis (Neilsen and Schreyer 1993) based on the macroscopic bulk tangent stiffness tensor or the corresponding acoustic tensor. One option is to adopt the limit point criterion, as used in (Massart et al. 2007; Svenning et al. 2019; Nguyen et al. 2012a), that is

$$ \det(D^M) = 0 \qquad (27) $$

which, once being fulfilled, will lead to the moment of bifurcation. However, this criterion could not directly yield the localization direction. Besides, use of the limit point criterion might be inappropriate for complex loadings that involve shear and/or compression stresses (Nguyen et al. 2012a) which leads to the non-symmetry of the tangent stiffness tensor.
A stricter criterion, i.e. loss of ellipticity (Hill 1962) which indicates the existence of a discontinuity in the velocity gradient in solids, serves as an alternative and is widely adopted in literature (Oliver et al. 2015; Sánchez et al. 2013; Souza and Allen 2011; Coenen et al. 2012a; Mercatoris et al. 2009b). This criterion can be phrased as:

$$ \text{Find } \mathbf{n}_d^M \Rightarrow \mathbf{Q} \cdot \mathbf{m} = \mathbf{0} \quad \text{at the first time } t_B \text{ for some } \mathbf{m} \\ \text{with } \mathbf{Q} = \mathbf{n}_d^M \cdot \mathbf{D}^M \cdot \mathbf{n}_d^M \quad \text{and} \quad ||\mathbf{n}_d^M|| = ||\mathbf{m}|| = 1 \qquad (28) $$

where $\mathbf{Q}$ is the acoustic tensor; $\mathbf{m}$ is the polarization vector indicating the initial displacement jump velocity. Since Condition (28) detects the singularity of the acoustic tensor, the equivalent form

$$ f_Q(\mathbf{n}_d^M) \stackrel{\text{def}}{=} \det(\mathbf{Q}) = 0 \qquad (29) $$

is commonly adopted as an alternative.

In practice, given the high-nonlinearity of $f_Q(\mathbf{n}_d^M)$, a sweeping procedure is usually used with the parameterization of the discontinuity normal vector $\mathbf{n}_d^M$ in terms of spherical coordinates to obtain the solution. After each deformation increment, the global minimum of $f_Q(\mathbf{n}_d^M)$ is found by sweeping all spherical angles of $\mathbf{n}_d^M$ (or by some iterative algorithm (Oliver et al. 2010) such as Newton's method) and then it is checked against Condition (29). Once the global minimum turns zero, it is identified that bifurcation occurs and $\mathbf{m}$ is also computed as the zero-eigenvector of the acoustic tensor. Macroscopic localization is subsequently introduced.

Due to the finite load increments in the finite element scheme, it will generally be the case that a perfect match of Condition (29) never happens. Instead, a negative global minimum of $f_Q(\mathbf{n}_d^M)$ is detected. Due to the fact that the shape of function $f_Q(\mathbf{n}_d^M)$ is evolving along the course of loading as the damage gradually accumulates, one scenario where the global minimum changes from one spherical angle of $\mathbf{n}_d^M$ to another angle as soon as it turns negative could happen, leading to a wrong prediction of the macroscopic localization direction. The same fact has been reported by Massart et al. (2007) and Mercatoris et al. (2009b) who tracked the evolution of the eigenspectrum of the acoustic tensor for a microscopic mansonry-structure model. In these papers, the authors proposed to take the last local minimum that makes eigenspectrum turn negative when the limit point is reached, instead of the global minimum, to detect the bifurcation orientation which matches the microstructural damage patterns. However, when multiple negative local minima exist at that moment this procedure may be unsuccessful in identifying the critical direction.

In order to systematically obtain a correct prediction regardless of the loading scheme, an extra condition is proposed here to reinforce Conditions (28) and (29) with the capability to discern the correct localization angle more robustly. At the load step when $f_Q(\mathbf{n}_d^M)$ turns negative for the first time, if multiple negative local minima are detected, we select that non-positive stationary point of which the failure mode is mostly aligned with the bulk strain field $\varepsilon^M$:

$$ \text{Find } (\mathbf{m}, \mathbf{n}_d^M) \Rightarrow \max(\varepsilon^M : (\mathbf{m} \otimes \text{sym} \mathbf{n}_d^M)). \qquad (30) $$

The performance of the proposed additional condition is further illustrated in Section 4.1.

## 3 Numerical implementations

Because the proposed framework deals with a local transition from bulk stress-strain behavior to cohesive traction-separation behavior, a discrete scheme that introduces a discontinuity on the fly is a natural match, such as can be achieved with XFEM (Moës and Belytschko 2002), the phantom node method (Song et al. 2006) or CutFEM (Burman et al. 2015). In this work, the phantom node method is chosen on the macroscale to model discontinuities. On the microscale, considering the potential high number of micro-cracks present in the RVE, a model that inserts cohesive elements into the mesh during the simulation (Camacho and Ortiz 1996; Liu et al. 2021) is adopted. Strong periodic BCs are applied on all RVEs with the acknowledgment of other BCs in the literature (e.g. Coenen et al. 2012b; Larsson et al. 2011; Svenning et al. 2017a; Goldmann et al. 2018) which aim to alleviate the shortcomings of strong periodic BCs when localization takes place. The
current implementation is for 2D macro-problems allowing for use of 3D micromodel in case of macroscopic plane strain conditions. Below, Voigt notation is used for definition of stress and strain and matrix notation is used for matrix-vector operations.

### 3.1 Macroscale model

The phantom node method simulates a cohesive crack by duplicating elements along the crack path on top of existing elements. Consider a three node triangular element with original nodes $[n_1, n_2, n_3]$ and one bulk Gauss point (GP) occupying a domain $\Omega^M$. In the pre-peak regime, an RVE is associated to the original bulk GP, being denoted as bulk RVE with classical FOCH. After microscopic damage initiates and accumulates at some point, the macroscopic localization is detected. At this moment, the element is crossed by a newly inserted crack segment $\Gamma_d^M$ which divides the element into two subdomains, $\Omega_A^M$ and $\Omega_B^M$, where $\Omega_A^M \cup \Omega_B^M = \Omega^M$. The original element is now replaced by two new elements A and B. Nodes are duplicated, except for the nodes on the element edge that contains the crack tip. Eventually, three new nodes are added to the current mesh, i.e. $\tilde{n}_1, \tilde{n}_2, \tilde{n}_3$. The connectivities are, respectively, nodes$_A = [n_1, \tilde{n}_2, \tilde{n}_3]$ and nodes$_B = [\tilde{n}_1, n_2, n_3]$. Both new elements are partially active within the subdomains $\Omega_A^M$ and $\Omega_B^M$, respectively. This leads to the discontinuous displacement description in the element domain $\Omega^M$:

$$ \mathbf{u}^M(\mathbf{x}^M) = \begin{cases} \mathbf{N}(\mathbf{x}^M)\mathbf{u}_A, & \mathbf{x}^M \in \Omega_A^M \\ \mathbf{N}(\mathbf{x}^M)\mathbf{u}_B, & \mathbf{x}^M \in \Omega_B^M \end{cases} \qquad (31) $$

where $\mathbf{N}(\mathbf{x}^M)$ is the standard shape function matrix. The difference between displacement fields of both elements along the crack is defined as the displacement jump:

$$ \llbracket \mathbf{u}(\mathbf{x}^M) \rrbracket^M = \mathbf{N}(\mathbf{x}^M) (\mathbf{u}_A - \mathbf{u}_B), \quad \mathbf{x}^M \in \Gamma_d^M. \qquad (32) $$

The strain field is derived by Equation (8),

$$ \varepsilon^M(\mathbf{x}^M) = \begin{cases} \varepsilon_A(\mathbf{x}^M) = \mathbf{B}(\mathbf{x}^M)\mathbf{u}_A, & \mathbf{x}^M \in \Omega_A^M \\ \varepsilon_B(\mathbf{x}^M) = \mathbf{B}(\mathbf{x}^M)\mathbf{u}_B, & \mathbf{x}^M \in \Omega_B^M \end{cases} \qquad (33) $$

where $\mathbf{B}(\mathbf{x}^M)$ is the standard strain-displacement matrix. The bulk macrostrain in at the location of the discontinuity is further defined as the average of the bulk strain at the opposite sides of the crack

$$ \varepsilon_b^M = \frac{1}{2}\mathbf{B}(\mathbf{u}_A + \mathbf{u}_B) \qquad (34) $$

and the smeared strain from the macroscopic displacement jump is evaluated as

$$ \varepsilon_f^M = \frac{1}{h} \mathbf{A} \mathbf{N}(\mathbf{u}_A - \mathbf{u}_B), \quad \text{where} \quad \mathbf{A}^\top = \begin{bmatrix} \mathbf{n}_{d1}^M & 0 & \mathbf{n}_{d2}^M \\ 0 & \mathbf{n}_{d2}^M & \mathbf{n}_{d1}^M \end{bmatrix}, \qquad (35) $$

and h is evaluated from the underlying RVE by using Equation (15).

Separate integration schemes are constructed for both partially-active elements A and B as well as for the cohesive crack segment, see Figure 3. Triangulation with three GPs are used for the bulk elements to calculate integration weights. Note that in the case that the constant strain triangular element is used only one bulk RVE is present in the uncracked element, which can be cloned to all newly allocated bulk and cohesive GPs, with the deformation state at the moment of bifurcation. Besides, for both element A and B all GPs are associated with the same cloned RVE. In the case where more than one GP is used for uncracked macro-elements, transfer of state variables between old and new GPs is required. However, since in multiscale modeling the state variable is an RVE it is not possible to do the mapping algebraically. Thus, a cloning operation is performed based on the distance between old and new GPs. The new RVE is copied from the closest old GP.

If bulk material around the crack is assumed to unload elastically, then one can even just use the secant stiffness for the new bulk GPs, as discussed in (Nguyen et al. 2012b). In some
<figure><img src="image_5.png" /><figcaption>Figure 3 Overlapping elements in phantom node model with cohesive segments on the macroscale.</figcaption></figure>

situation using the secant stiffness is also beneficial for robustness since sudden considerable unloading of bulk RVEs may cause convergence problems in the microscopic BVPs. However for general loading scenarios, the possibility of additional nonlinearity in bulk GPs of cracked elements cannot be excluded. In the current implementation, cloning is first attempted. If difficult convergence is observed in a newly cloned micromodel, the secant stiffness is used instead.

It follows from Equation (21) that the stiffness matrix contribution from the cohesive traction consists of two parts, one related to the macroscopic strain and one related to the macroscopic displacement jump. Together with the bulk contribution from the active part of each of the two overlapping elements, the consistent tangent matrix for the two overlapping elements is

$$ K = \begin{bmatrix} \mathbf{K}_A^{\text{bulk}} & 0 \\ 0 & \mathbf{K}_B^{\text{bulk}} \end{bmatrix} + \begin{bmatrix} \mathbf{K}_{\varepsilon}^{\text{coh}} & \mathbf{K}_{\varepsilon}^{\text{coh}} \\ -\mathbf{K}_{\varepsilon}^{\text{coh}} & -\mathbf{K}_{\varepsilon}^{\text{coh}} \end{bmatrix} + \begin{bmatrix} \mathbf{K}_{\llbracket u \rrbracket}^{\text{coh}} & -\mathbf{K}_{\llbracket u \rrbracket}^{\text{coh}} \\ -\mathbf{K}_{\llbracket u \rrbracket}^{\text{coh}} & \mathbf{K}_{\llbracket u \rrbracket}^{\text{coh}} \end{bmatrix}, \quad (36) $$

with the matrices

$$ K_i^{\text{bulk}} = \int_{\Omega_i^M} \mathbf{B}^{\text{T}} \mathbf{D}^M \mathbf{B} \mathrm{d}\Omega, \quad i = A, B, \quad K_{\varepsilon}^{\text{coh}} = \frac{1}{2} \int_{\Gamma_d^M} \mathbf{N}^{\text{T}} \mathbf{A}^{\text{T}} \mathbf{D}^M \mathbf{B} \mathrm{d}\Gamma, \quad K_{\llbracket u \rrbracket}^{\text{coh}} = \frac{1}{h} \int_{\Gamma_d^M} \mathbf{N}^{\text{T}} \mathbf{A}^{\text{T}} \mathbf{D}^M \mathbf{A} \mathrm{N} \mathrm{d}\Gamma. $$

The homogenized bulk tangent matrix $\mathbf{D}^M$ can be obtained through static condensation which excludes the dependent degree of freedom (dofs) on the global stiffness matrix of the microscale problem at the converged state, that is (Kouznetsova et al. 2001)

$$ \mathbf{D}^M = \frac{1}{|\Omega^m|} \mathbf{H} \mathbf{K}_{pp}^{\star} \mathbf{H}^{\text{T}} \quad \text{with} \quad \mathbf{H} = [\mathbf{H}_1 \ \mathbf{H}_2 \ \mathbf{H}_4], \quad \mathbf{H}_q^{\text{T}} = \frac{1}{2} \begin{bmatrix} 2x_{q1}^m & 0 & x_{q2}^m \\ 0 & 2x_{q2}^m & x_{q1}^m \end{bmatrix}, \quad q = 1, 2, 4, \quad (37) $$

where $x_q^m$ denotes the coordinates of controlling nodes $q$ in the RVE as shown in Figure 4. $K_{pp}^{\star}$ is

<figure><img src="image_6.png" /><figcaption>Figure 4 Strong periodic BCs on RVE, which reads $\mathbf{u}_{\Gamma_{01}} = \mathbf{u}_{\Gamma_{00}} + \mathbf{u}_2 - \mathbf{u}_1$ and $\mathbf{u}_{\Gamma_{11}} = \mathbf{u}_{\Gamma_{10}} + \mathbf{u}_4 - \mathbf{u}_1$.</figcaption></figure>

the microscopic reduced stiffness matrix with $p$ representing dofs associated to the controlling nodes.
The macroscale solution algorithm is Newton-Raphson iterative procedure with the extra check for crack insertion or propagation at the end of each converged load step. Adaptive load stepping is also used to help convergence. For more details on the solution algorithm, see (van der Meer 2012).

## 3.2 Microscale model

On the microscale, cohesive elements are used to model cracks. In order to reduce the computational cost and improve robustness in presence of extensive microcracking, the cohesive elements are placed between neighbouring finite elements during the simulation, wherever the stress measure exceeds a critical value following the work by Camacho and Ortiz (1996).

Six-node triangular elements are used to discretize the RVE domain in 2D problems. Cracks may initiate at the mid-node along the element edge and always propagate along element edges, see Figure 5.

<figure><img src="image_7.png" /><figcaption>**Figure 5** Insertion of cohesive elements during simulation.</figcaption></figure>

An effective stress criterion is adopted to test the stress state for local failure (Camacho and Ortiz 1996)

$$ \bar{\sigma} \ge f_t \quad \text{with} \quad \bar{\sigma} = \begin{cases} \sqrt{(t_{dn}^m)^2 + \beta(|t_{ds}^m|^2)}, & \|u\|_n^m \ge 0 \\ \sqrt{\beta(|t_{ds}^m| - \mu|t_{dn}^m|)}, & \|u\|_n^m < 0, \end{cases} \qquad (38) $$

where $f_t$ is the tensile strength, $t_{dn}^m$ and $t_{ds}^m$ is the normal and shear traction component, $\beta$ is a shear stress factor and $\mu$ is the friction coefficient. This criterion is checked at every mid-node when equilibrium is met at each load step. After crack insertion the Newton-Raphson procedure is re-entered to find equilibrium again for the same load step.

Because the cohesive element is inserted at a non-zero stress state, it is necessary to adopt an initially rigid cohesive law. To achieve this without introducing a singularity at zero-opening, we shift the origin of a non-rigid cohesive law such that the traction at zero displacement jump and zero damage matches the surrounding bulk stress on the mid-node from the bulk element, as in (van der Meer et al. 2012).

It is worth noting that for using periodic BCs, special care must be given to the splitting nodes on the boundary of the RVE. As shown in Figure 6, because of the periodicity, a crack running through a boundary actually re-enters the domain on the opposite boundary. So splitting a node

<figure><img src="image_8.png" /><figcaption>**Figure 6** Periodic cracking pattern in pairs of RVE boundary nodes.</figcaption></figure>

at the edge of the domain must be consistent with the splitting of the node on the opposite edge in order to maintain the ability to apply periodic constraints. It means such a pair of opposing
nodes should be split within one load step. When the crack reaches the edge from one side only, it is kept closed just like the node that is at the tip of a crack. Subsequently, when a crack reaches the corresponding node on the opposite boundary, both nodes are split to achieve a crack that effectively crosses the boundary.

The model has also been extended to dealing with 3D problems where a slice of material is modelled to consider longitudinal shear case in fiber reinforced composites. A 2D mesh as input is extruded to a 3D mesh with only one element in the thickness direction. Thus six-node triangular elements are extended to twelve-node wedge elements. The stress evaluation in Equation (38) now considers the contribution from the area along element edges but is still evaluated at the mid-node. The cracking operation, however, still remains on the original 2D mesh which means the two opposite six-node faces of the wedge element always have a mirrored structure.

# 4 Numerical tests

## 4.1 Bifurcation mode analysis

In order to demonstrate the importance and the capability of the localization condition proposed in Section 2.4, microstructures with randomly distributed fiber inclusions inside matrix are analyzed for transverse tension (TT) and longitudinal shear (LS) cases associated to one macroscopic material point. Two micromodels are adopted in simulations, as shown in Figure 7. A 2D

<figure><img src="image_9.png" /><figcaption>**Figure 7** Micromodels for bifurcation mode analysis: (left) TT case; (right) LS case.</figcaption></figure>

plane-strain mesh is used for the TT case while a 3D slice mesh is used for the LS case. Both meshes are periodic and the periodic BCs can be applied through linear constraints between dofs. Material parameters of the two cases are summarized in Table 1. Fibers and matrix are

<table><caption>**Table 1** Material parameters for matrix and fibers.</caption><thead><tr><td>Component</td><td>Young's modulus E (MPa)</td><td>Poisson's ratio v</td></tr></thead><tbody><tr><td>Fiber</td><td>12970</td><td>0.46</td></tr><tr><td>Matrix</td><td>5070</td><td>0.35</td></tr></tbody></table>

characterized as isotropic linear elastic. Matrix cracking and fiber/matrix debonding are considered with the cohesive zone model with the shifted cohesive law described in Section 3.2. Interfacial parameters are given in Table 2.

<table><caption>**Table 2** Material parameters of cohesive law in Section 4.1.</caption><thead><tr><td>Crack location</td><td>Penalty stiffness K<sub>p</sub> (N/mm<sup>3</sup>)</td><td>Fracture strength f<sub>t</sub> (MPa)</td><td>Fracture energy G<sub>c</sub> (N/mm)</td><td>Interaction coefficient β</td><td>Friction coefficient μ</td></tr></thead><tbody><tr><td>Matrix</td><td>1.0 × 10<sup>7</sup></td><td>121</td><td>0.09</td><td>0.4</td><td>0.1</td></tr><tr><td>Fiber/Matrix</td><td>1.0 × 10<sup>7</sup></td><td>42</td><td>0.04</td><td>0.4</td><td>0.1</td></tr></tbody></table>

The macrostrain is applied on controlling nodes through deformation increments. The determinant of the homogenized acoustic tensor is being monitored throughout the simulation until a few load steps into the softening regime. For the 2D macroscale problem, the normal vector $\mathbf{n}_d^M$ is expressed as $\mathbf{n}_d^M = [\cos \theta \sin \theta]^T$ with $\theta \in [-90^\circ, 90^\circ]$, which subsequently leads to
<figure><img src="image_10.png" /><figcaption>**Figure 8** Bifurcation mode analysis: [left column] $f_Q(\theta)$ normalized by its maximum value; [right column] deformation of micromodels at time of bifurcation. [top row] TT case; [bottom row] LS case; [*] global minimum, [•] correct localization angle.</figcaption></figure>

$f_Q(\mathbf{n}_d^M) = f_Q(\theta)$. The evolution of $f_Q(\theta)$ for both cases is shown in left column of Figure 8. The determinant is normalized by the maximum value. For each curve, the global minimum is marked by a red star while the actual orientation of the microscopic localization band is indicated with a black dot. The corresponding deformations ($u_{x1}$ and $u_{x3}$, respectively) of micromodels at time of bifurcation are shown in the right column of Figure 8.

By choosing the global minimum of $f_Q(\theta)$ as the desired solution, the predicted normal vector angle for the TT case is 2.537°, which is accurate enough. However, for the LS case, the obtained angle 0° does not correspond to the orientation of the actual localization band in the micromodel. If we consider the minimum of $f_Q(\theta)$ for later steps which is representative for what would need to be used in case of larger time steps the global minimum deviates from the correct for both cases. However, by adopting the proposed Condition (30) with the extra information from the bulk strain field, the correct localization angle is always obtained.

## 4.2 Micromodel size independence

### 4.2.1 Regular micro-crack pattern

To demonstrate the objectivity of homogenized results with respect to the micromodel size, a strip under uniaxial tensile loading with pre-defined regular micro-crack pattern is solved by the proposed homogenization framework. A horizontal displacement $\delta$ is applied at the right edge. As shown in Figure 9(a), the tested strip with a dimension of $3d^M \times d^M$ ($d^M = 0.0114$ mm) is discretized with constant strain triangular elements. The two-scale analysis is only restricted to the central zone of two elements while the left and right zones are modelled by an isotropic linear elastic law with Young's modulus $E = 3000$ MPa and Poisson's ratio $\nu = 0.18$. The central vertical dashed line indicates the pre-defined potential crack location. Note that in such a simple case, the macro-crack is fixed rather than detected from the micromodel. Besides, considering the pre-defined vertical crack pattern, the parameter $h$ is set equal to the RVE length.

The adopted micromodel is shown in Figure 9(b) with a dimension of $d^m \times d^m$. Four different sizes are investigated, i.e. $d^m = \frac{1}{5}d^M$, $\frac{1}{10}d^M$, $\frac{1}{30}d^M$ and $\frac{1}{50}d^M$. The material is assumed to be
<figure><img src="image_11.png" /><figcaption>**Figure 9** Uniaxial tensile test with pre-defined micro-crack patterns: (a) macroscale problem description (shown in FE discretization); (b) micromodel which is allowed to open at the central vertical line and discretized by six-node triangular elements.</figcaption></figure>

isotropic linear elastic with Young's modulus $E = 2000$ MPa and Poisson's ratio $\nu = 0.18$. Failure can take place as a cohesive crack along the central vertical mesh line with cohesive properties as given in Table 3. A plane strain condition is assumed for both macroscale and microscale.

<table><caption>**Table 3** Material parameters of cohesive law for example in Section 4.2.1. See Table 2 for parameter definition.</caption><thead><tr><th>K<sub>p</sub> (N/mm<sup>3</sup>)</th><th>f<sub>t</sub> (MPa)</th><th>G<sub>c</sub> (N/mm)</th><th>coefficient β</th><th>coefficient μ</th></tr></thead><tbody><tr><td>1.0 × 10<sup>7</sup></td><td>3.0</td><td>0.0005</td><td>0.4</td><td>0.1</td></tr></tbody></table>

Figure 10 shows the macroscale force-displacement curves and homogenized traction-separation relations in normal direction for all micromodels. It is evident that all results

<figure><img src="image_12.png" /><figcaption>**Figure 10** Objective results with respect to the size of the micromodel with regular crack pattern. (a) force-displacement curves; (b) homogenized cohesive laws in normal direction.</figcaption></figure>

coincide, which confirms that the homogenized traction-separation relation is independent of the dimensions of the micromodel.

### 4.2.2 Complex micro-crack pattern

Micromodels that consist of a matrix with a periodic array of fiber inclusions and an unstructured mesh are now considered in the same uniaxial tensile loading strip test to show the micromodel size independence in the case of complex micro-crack patterns. The tested strip of $3d^M \times d^M$ ($d^M = 0.1$ mm), as shown in Figure 11(a), is discretized with quadrilateral elements with four integration points. The two-scale analysis is again restricted to the central element while other elements are modelled with the homogenized linear elastic properties from the micromodel. Due to the random nature of the micro-crack patterns, the macroscopic crack orientation is now predicted at the onset of localization rather than pre-defined.

The three micromodels are shown in Figure 11(b) with $d^m = 0.01$ mm. Based on the arrangement of fibers, they are coined as micromodel $1 \times 1, 2 \times 2$ and $3 \times 3$, respectively. Fibers and matrix are characterized as isotropic linear elastic with properties given in Table 1. Matrix cracking and fiber/matrix debonding are handled by the cohesive zone model with the interfacial parameters given in Table 2 except for the fracture strength of matrix being 69 MPa and fracture
<figure><img src="image_13.png" /><figcaption>**Figure 11** Uniaxial tensile test with complex micro-crack pattern: (a) macroscale problem description (shown in FE discretization); (b) micromodels discretized by six-node triangular elements; (c) final micro-failure patterns (black lines indicate completely failed cracks while grey lines indicate unloading cracks).</figcaption></figure>

energy for matrix and interface being 0.15 N/mm and 0.1 N/mm. A plane strain condition is assumed for both macroscale and microscale.

The final failure patterns of micromodels are shown in Figure 11(c). As expected, randomly distributed cohesive cracks emerge and gradually localize into one opening crack plane which separates the microscopic domain into two parts. Note that the parameter *h* in this test is set equal to the RVE length as well. However, it is very close to the value that would be obtained with Equation (16) since in this mode-I case the predicted macroscopic crack orientation (−90.67°, −90.52° and −88.54°) is almost vertical. The strategy to determine *h* in complex micro-crack patterns is further discussed in Section 4.4. It should be noted that the calculation of *h* does not influence the examination of micromodel size objectivity as long as it follows the same definition in all test cases.

Figure 12 shows the macroscale force-displacement curves and homogenized traction-separation relations in normal direction for all micromodels. Although small variations can be

<figure><img src="image_14.png" /><figcaption>**Figure 12** Objective results with respect to the size of the micromodel with complex crack pattern: (a) force-displacement curves; (b) homogenized cohesive laws in normal direction.</figcaption></figure>

observed due to the different micro-failure patterns, the results are very similar, which again confirms the objectivity with respect to the micromodel size.

### 4.3 Three point bending test

A three point bending test with an array of voids over the height at mid-span is considered to check the ability of the proposed approach to simulate a problem with a propagating crack.
DNS analysis is carried out to provide a reference solution for verifying the accuracy of the multiscale framework. Furthermore, the insensitivity of the macroscale response with respect to the macroscale mesh size is examined.

Figure 13 shows the set-up of the macroscale problem and the micromodel adopted for the two-scale simulation. A notch is located at the center bottom of the beam to trigger cracking. A

<figure><img src="image_15.png" /><figcaption>**Figure 13** Three point bending test: (a) macroscale problem description; (b) micromodel which is allowed to open at the central vertical line.</figcaption></figure>

vertical prescribed displacement increment $\delta$ is applied at the top mid-point. The microstructure consists of an isotropic linear elastic material with a central void. The Young's modulus is $E = 30000$ MPa and Poisson's ratio $\nu = 0.15$. Cohesive cracks are only allowed along the middle vertical line that is highlighted in red. The used shifted cohesive law properties are as same as the ones in Table 3 with the exception of the fracture energy which is set to $G_c = 1.0$ N/mm. A plane-strain condition is assumed for both scales.

Figure 14(a) shows the finite element mesh for the DNS simulation. Note that cohesive elements are also inserted on the fly at the central line with the same shifted cohesive law as used in the micromodel. Three different macroscale meshes are used for FE² analysis, as shown in

<figure><img src="image_16.png" /><figcaption>**Figure 14** Finite element mesh for three point bending test: (a) DNS analysis mesh; (b) macroscale mesh for two-scale CH analysis with different discretizations along the central band.</figcaption></figure>

Figure 14(b). In order to restrict the computational cost two-scale CH simulation only takes place in the central row of elements, which is in agreement with the DNS set-up. The number of elements used for the FE² analysis region is 80, 120 and 160 elements, respectively. For both DNS and FE² analysis, the elements away from the mid-plane are modelled by the homogenized elastic
stiffness matrix that obtained from a priori homogenization of the micromodel (RVE). As in previous examples, the *h* in FE² analysis is set equal to the RVE length.

Force-displacement curves are given in Figure 15. It is observed that the sensitivity of the

<figure><img src="image_17.png" /><figcaption>**Figure 15** Force-displacement curves of three point bending tests.</figcaption></figure>

macroscale response to the macro mesh size is negligible although the smoothness of the curve is affected by the size of the elements, which proves the objectivity of the proposed framework with respect to the macroscale discretization. Along the curve, small drops are present at each step when crack segment is inserted, which is due to the fact that the crack grows element-wise in the macromodel. Comparing the DNS and FE² analysis results, an exact agreement can be observed, which verifies the reliability of the proposed two-scale CH framework for crack growth analysis.

### 4.4 Energy equivalence: study on scale transition parameter *h*

As discussed in Section 2.2, the scale transition factor *h* directly influences the energy equivalence between two scales and subsequently the reliability of homogenized results. In the earlier examples in this paper, *h* has been set equal to the RVE length. In this section, three cases with different pre-defined micro-crack patterns and one case with random micro-crack pattern are studied to investigate the proper evaluation of the parameter *h*. The cases are selected to represent non-trivial crack growth scenario's in complex micromodels, where cracks will be tortuous or dispersed and the crack pattern may still evolve during softening. Note that although a discrete cracking model is adopted on the microscale, the test could also be performed with a continuum damage model. Figure 17, Figure 18, Figure 19 and Figure 20 show the test set-up and FE meshes for case A, case B, case C and case D, respectively. For case A, the material failure emerges as inclined straight cohesive cracks while for case B zigzag cracks separate the material. For both cases, the material is isotropic linear elastic with Young's modulus *E* = 2000 MPa and Poisson's ratio *ν* = 0.18. For case C, both cohesive cracks and matrix plasticity are present. A hardening plasticity model developed by Melro et al. (2013) is adopted with the hardening curve shown in Figure 16 and plastic Poisson's ratio being 0.22. Two non-matching cracks are

<figure><img src="image_18.png" /><figcaption>**Figure 16** Hardening curves for the plasticity model.</figcaption></figure>

pre-defined, such that also in the post-peak response, plastic deformations will continue to develop in part of the domain. Cohesive properties for case A to case C are given in Table 4. In
<table><caption>**Table 4** Material parameters of cohesive cracks in Section 4.4. See Table 2 for parameter definition.</caption><thead><tr><th>Case</th><th>K<sub>p</sub> (N/mm<sup>3</sup>)</th><th>f<sub>t</sub> (MPa)</th><th>G<sub>c</sub> (N/mm)</th><th>β</th><th>μ</th></tr></thead><tbody><tr><td>A & B</td><td>1.0 × 10<sup>8</sup></td><td>3.0</td><td>0.0005</td><td>1.0</td><td>0.1</td></tr><tr><td>C</td><td>1.0 × 10<sup>8</sup></td><td>3.0</td><td>0.0002</td><td>1.0</td><td>0.1</td></tr></tbody></table>

case D, a microstructure with periodic heterogeneous material is considered, while the material properties and cohesive laws are the same as in the example of Section 4.2.2.

<figure><img src="image_19.png" /><figcaption>**Figure 17** *h* study—case A: (a) macroscale problem shown in FE discretization; (b) FE mesh for DNS analysis with red lines dictating initial zero-thickness cohesive elements and the adopted micromodel.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18** *h* study—case B: (a) macroscale problem shown in FE discretization; (b) FE mesh for DNS analysis with red lines dictating initial zero-thickness cohesive elements and the adopted micromodel.</figcaption></figure>

On the macroscale, uniaxial tensile tests are considered with different specimen dimensions. As in Section 4.2, the FE² approach is only active in the central row of elements while the others are modelled as elastic with the same stiffness parameters as the material on the microscale. FE meshes for DNS analysis are also given together with the adopted micromodel mesh for each case. As shown, the DNS mesh is constructed by repeating the micromodel mesh in the central area. Zero-thickness cohesive elements are initially inserted along potential crack locations (shown in red) in the DNS mesh of case A to case C with a bilinear cohesive law possessing the same parameters in Table 4. However, for case D, the cohesive elements are inserted on the fly as in the micromodel. A plain strain condition is adopted on both scales. A prescribed displacement increment is applied horizontally on the right edge of the bar specimen.
<figure><img src="image_21.png" /><figcaption>**Figure 19** *h* study—case C: (a) macroscale problem shown in FE discretization; (b) FE mesh for DNS analysis with red lines dictating initial zero-thickness cohesive elements and the adopted micromodel.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20** *h* study—case D: (a) macroscale problem shown in FE discretization; (b) FE mesh for DNS analysis and the adopted micromodel.</figcaption></figure>

As mentioned in Section 2, *h* is introduced into the current framework as a characteristic length to avoid spurious RVE size sensitivity when entering the softening regime. From a physical point of view, this characteristic length represents the amount of unloading material perpendicular to the crack path. Thus, we propose to endow *h* with such geometrical interpretation and calculate it as the crack spacing between crack paths in the periodically deformed micro-structure, as indicated by *d*~*c*~ in Figure 21. This is in line with the idea of a macroscopic localization band width, from (Massart et al. 2007), into which the failure is localized. Note that this computing method is equally applicable to the case where continuum damage model is adopted in the microscale model for which the distance between damage bands should be considered. The definition by Equation (14) based on the length of loading cracks is also examined for all tests.

Horizontal displacement results are given for case B in Figure 22 to show that in the DNS test fracture eventually localizes into the opening of only one crack plane rather than multiple parallel cracks. In other words, only one line of the cohesive elements is open while the others just remain inactive.

Force-displacement curves from all tests are shown in Figure 23. It is evident that using $h = d_c$ always gives reliable results compared to DNS simulations while Equation (14) does not
<figure><img src="image_23.png" /><figcaption>**Figure 21** Periodic-crack spacing in micro-structures.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 22** Horizontal displacement (blue = 0 mm and yellow = 3.2 × 10⁻⁴ mm) at the end of simulation of DNS analysis in case B.</figcaption></figure>

guarantee energetically equivalent results for the zigzag or random cracking case and the case with plasticity. For case A, the same results are obtained from two *h* formulae which is not surprising because for this case the two formulae of *h* yield the same value. For case B, the crack tortuosity is not well captured by Equation (14) which could be improved by some correction based on the micro-crack information (Toro et al. 2016a). Regarding the case C, when using Equation (14) the non-smoothness of the force-displacement curve reflects the fact that the total length of the opening cohesive cracks is not constant during softening for this case. Updating the variable |Γ_d^m| and correcting for tortuosity becomes very cumbersome for this case. The tortuosity correction might become impossible except if the definition of the tortuosity starts from a known *d*_c in which case it is definitely more convenient to immediately set *h* = *d*_c. Similar results can be observed in case D, where the small value of *h* predicted by Equation (14) at the onset of localization makes the micromodel deform massively even for a small macroscopic displacement jump, leading to a sharply descending homogenized cohesive traction-separation relation. This comes from the fact that the most of microscopic cracks are actually loading at the bifurcation moment, as shown in Figure 24. Further evolution causes a sudden change of *h* to a very high value which leads to the spurious snapback on the macroscopic response. At the lowest point on the macroscopic force-displacement curve, it is even not possible to update *h* based on the loading crack length since all cracks are unloading at this moment, which ends the simulation. Compared to the DNS results, it is clear that with the *h* predicted by Equation (14) the results are very inaccurate. Thus, it is concluded that it is more sensible to adopt *h* as a constant value *d*_c at the onset of localization.

Furthermore, it can easily be deduced that setting *h* equal to the length of the micromodel (*h* = *d*^m) would work for cases B and C, where *d*_c = *d*^m, but not for case A, where *d*_c ≠ *d*^m. Finally, it is worth noting that for the 2D square micromodel case the proposed *h* = *d*_c is equivalent to

<figure><img src="image_25.png" /><figcaption>**Figure 23** Force-displacement curves of *h* study tests: (a) case A; (b) case B; (c) case C; (d) case D.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 24** Cohesive cracks in the micromodel at the moment of bifurcation in case D. Black lines indicate loading cracks while grey lines indicate unloading cracks.</figcaption></figure>

Equation (16). In case of predefined microscale failure mechanisms, $d_c$ can be determined in a natural way based on the geometry. In general multiscale analysis such as case D, where $h$ has to be detected based on the orientation of the localization band, Equation (16) can be used. For case A to C, it is confirmed that the acoustic tensor analysis does yield accurate macroscale localization directions, i.e. a vertical macroscopic crack for cases B and C and an inclined crack under 45° for case A.

## 4.5 Arbitrary macrocrack test

Finally, in order to show the ability of the proposed multiscale framework to simulate crack propagation for a case with non-trivial crack path, a mixed-mode failure test is performed. A similar test has been carried out numerically in (Wu 2018) with a mono-scale damage model and studied experimentally in (Gálvez et al. 1998). Although with different material parameters, the macrocrack is expected to follow an analogous path which emerges from the notch top and propagates in a curved trajectory toward the beam top. Macroscopic dimensions of the case are given in Figure 25. The microstructure consists of a matrix embedded with a circular inclusion.

<figure><img src="image_27.png" /><figcaption>**Figure 25** Mix-mode failure test: (a) macroscale problem description; (b) micromodel discretized by six-node triangular elements.</figcaption></figure>

Both materials are isotropic and linear elastic with Young's modulus $E_{matrix}$ = 68 000 MPa and $E_{inclusion}$ = 38 000 MPa and Poisson's ratio $\nu$ = 0.2 for both. Failure is assumed to occur inside the matrix and along the matrix/inclusion interface with shifted cohesive laws in Table 5.

<table><caption>**Table 5** Material parameters of cohesive law for example in Section 4.5. See Table 2 for parameter definition.</caption><thead><tr><th>Crack location</th><th>K<sub>p</sub> (N/mm<sup>3</sup>)</th><th>f<sub>t</sub> (MPa)</th><th>G<sub>c</sub> (N/mm)</th><th>β</th><th>μ</th></tr></thead><tbody><tr><td>Matrix</td><td>1.0 × 10<sup>4</sup></td><td>4.0</td><td>0.12</td><td>0.4</td><td>0.1</td></tr><tr><td>Matrix/Inclusion</td><td>1.0 × 10<sup>4</sup></td><td>2.8</td><td>0.075</td><td>0.4</td><td>0.1</td></tr></tbody></table>

A macroscale mesh with four-node quadrilateral elements adopting 2 × 2 Gauss integration scheme is shown in Figure 26. To limit the computational cost, the FE² approach is only active in elements inside the orange area while the others are modelled as elastic with the homogenized stiffness parameters from the micromodel. A plane-strain condition is assumed for both scales. The analysis is performed using the dissipation-based arc-length method (Verhoosel et al. 2009) to tackle the potential snap-back behavior of the macroscopic response with an adaptive-stepping scheme (van der Meer et al. 2012). For calculating the scale transition parameter $h$, the localization angle is obtained via acoustic tensor analysis in combination with Equation (15).
The final crack pattern is shown in Figure 26, which captures the expected crack propagation path rather well. Snapshots of two deformed micromodels at time-step indicated by point A

<figure><img src="image_28.png" /><figcaption>**Figure 26** Macroscale mesh for mix-mode failure test and curved macroscale crack with snapshots of micromodels at time-step indicated by point A in Figure 27.</figcaption></figure>

in Figure 27 are also shown. In Figure 27, the macroscopic force-displacement curve and the

<figure><img src="image_29.png" /><figcaption>**Figure 27** Macroscale response: (a) force-displacement curve; (b) homogenized traction-separation law of first cohesive point on macro-crack.</figcaption></figure>

homogenized traction-separation relation for the first cohesive GP are shown. Clearly, snap-back behavior is captured. However, some spurious snapback occurs due to the remeshing from phantom-node procedure when a crack inserts into a relatively coarse mesh (van der Meer 2012). From the obtained traction-separation relation, it can be confirmed that the problem indeed experienced a mixed-mode cracking due to the same order of magnitude of shear and normal tractions. Moreover, it is observed that at the end of analysis, traction tends to be constant rather than decreasing to zero, which prevents further failure at the corresponding material point. This artificial locking is caused by the imposing of strong periodic BCs as can be seen in the snapshot in Figure 26. Due to the non-conforming macroscopic crack direction with the periodicity, total failure of the RVE domain is not always guaranteed to occur. Remedies could be achieved by choosing alternative BCs which aim to resolve such problems. Note that this phenomenon does not violate the validity of current framework in terms of simulating a curved crack. However, it does highlight that it is not straightforward to choose an appropriate value of scale transition parameter *h* for the general case. The most damaged micromodel in Figure 26 has a crack pattern similar to case C from Section 4.4, but may still develop to a scenario with multiple cracks
through the RVE, similar to case A. The question how to find the most suitable *h* for such case is
related to the choices that need to be made to resolve the boundary condition issue.

5 Conclusion

In this paper, a two-scale computational homogenization framework is presented to address the strain softening issue in the multiscale material failure modeling problem. Discrete macroscopic crack onset and propagation for heterogeneous materials are derived from the micromodel in which the underlying microstructure undergoes accumulating micro-fracture failure. The framework keeps high compatibility with well-established first-order computational homogenization code in the sense that the down-scaling and homogenization procedures for both are the same for bulk and cohesive integration points. The key scale transition parameter which also serves as a numerical characteristic length is examined and endowed with a geometrical interpretation which leads to a proper computing method. The common acoustic tensor bifurcation criterion is reinforced by an additional criterion to help detect the localization mode more robustly. Various numerical tests have been performed to verify the consistency of the homogenized results from the presented computational homogenization framework through comparison with direct numerical simulation. Objectivity of the results in terms of macroscale finite element mesh size and micromodel size has been confirmed as well. Finally, the ability of the framework to track an arbitrary propagating macro-crack is demonstrated. The framework is generic in the sense that no assumptions are made in the scale transition on the microscale constitutive models. Cohesive cracks have been used on the micromodel, with and without plasticity in the surrounding bulk, and the framework is equally applicable to cases with continuum damage in the bulk material.

An open point concerns on-the-fly detection of the orientation of the microscopic localization band as is needed for the scale transition parameter *h*. It should be noted however, that strong periodic BCs do not support arbitrary localization orientations. It is recommended to approach these two problems as interconnected in future research: more appropriate formulations for the boundary conditions may already require detection of the orientation, after which determination of *h* is relatively straightforward.

References

Bažant, Z. (2010). Can multiscale-multiphysics methods predict softening damage and structural failure? *International Journal for Multiscale Computational Engineering* 8(1):61-67. [DOI],
[HAL].

Belytschko, T., S. Loehnert, and J. H. Song (2008). Multiscale aggregating discontinuities: A method for circumventing loss of material stability. *International Journal for Numerical Methods in Engineering* 73(6):869-894. [DOI].

Bosco, E., V. Kouznetsova, and M. Geers (2015). Multi-scale computational homogenization-localization for propagating discontinuities using X-FEM. *International Journal for Numerical Methods in Engineering* 102(3-4):496-527. [DOI].

Budarapu, P. R., X. Zhuang, T. Rabczuk, and S. Bordas (2019). Multiscale modeling of material failure: Theory and computational methods. *Advances in Crystals and Elastic Metamaterials, Part 2*. Vol. 52. Advances in Applied Mechanics. Elsevier. Chap. 1, pp 1–103. [DOI].

Burman, E., S. Claus, P. Hansbo, M. Larson, and A. Massing (2015). CutFEM: Discretizing geometry and partial differential equations. *International Journal for Numerical Methods in Engineering* 104(7):472-501. [DOI], [OA].

Camacho, G. and M. Ortiz (1996). Computational modelling of impact damage in brittle materials. *International Journal of Solids and Structures* 33(20-22):2899-2938. [DOI].

Coenen, E., V. Kouznetsova, and M. Geers (2010). Computational homogenization for heterogeneous thin sheets. *International Journal for Numerical Methods in Engineering*. [DOI].

Coenen, E., V. Kouznetsova, and M. Geers (2012a). Multi-scale continuous-discontinuous framework for computational homogenization localization. *Journal of the Mechanics and Physics of Solids* 60(8):1486-1507. [DOI].

Coenen, E., V. Kouznetsova, and M. Geers (2012b). Novel boundary conditions for strain localization analyses in microstructural volume elements. *International Journal for Numerical Methods in*
*Engineering* 90(1):1–21. [DOI].

Feyel, F. and J. L. Chaboche (2000). FE² multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials. *Computer Methods in Applied Mechanics and Engineering* 183(3-4):309–330. [DOI].

Fish, J. (2006). Bridging the scales in nano engineering and science. *Journal of Nanoparticle Research* 8:577–594. [DOI].

Gitman, I., H. Askes, and L. J. Sluys (2007). Representative volume: Existence and size determination. *Engineering Fracture Mechanics* 74(16):2518–2534. [DOI].

Gitman, I., H. Askes, and L. J. Sluys (2008). Coupled-volume multi-scale modelling of quasi-brittle material. *European Journal of Mechanics-A/Solids* 27(3):302–327. [DOI].

Goldmann, J., J. Brummund, and V. Ulbricht (2018). On boundary conditions for homogenization of volume elements undergoing localization. *International Journal for Numerical Methods in Engineering* 113(1):1–21. [DOI].

Gálvez, J., M. Elices, G. Guinea, and J. Planas (1998). Mixed mode fracture of concrete under proportional and nonproportional loading. *International Journal of Fracture* 94(3):267–284. [DOI].

Helfen, C. and S. Diebels (2014). Computational homogenisation of composite plates: Consideration of the thickness change with a modified projection strategy. *Computers and Mathematics with Applications*. [DOI], [OA].

Hill, R. (1962). Acceleration waves in solids. *Journal of the Mechanics and Physics of Solids* 10(1):1–16. [DOI].

Hill, R. (1965). A self-consistent mechanics of composite materials. *Journal of the Mechanics and Physics of Solids* 13(4):213–222. [DOI].

Hirschberger, C., S. Ricker, P. Steinmann, and N. Sukumar (2009). Computational multiscale modelling of heterogeneous material layers. *Engineering Fracture Mechanics*. [DOI], [HAL].

Hofman, P., L. Ke, and F. van der Meer (2021). Circular Microstructural Volume Elements With Periodic Boundary Conditions for Localization Problems. *VIII Conference on Mechanical Response of Composites* (Online, 2021). [DOI], [OA].

Jalalvand, M., G. Czél, J. Fuller, M. Wisnom, L. Canal, C. González, and J. LLoorca (2016). Energy dissipation during delamination in composite materials – An experimental assessment of the cohesive law and the stress-strain field ahead of a crack tip. *Composites Science and Technology* 134:115–124. [DOI], [OA].

Khoei, A. and M. Saadat (2019). A nonlocal computational homogenization of softening quasi-brittle materials. *International Journal for Numerical Methods in Engineering*. [DOI].

Kouznetsova, V. (2002). Computational homogenization for the multi-scale analysis of multi-phase materials. PhD thesis. Technische Universiteit Eindhoven. [DOI].

Kouznetsova, V., W. Brekelmans, and F. Baaijens (2001). An approach to micro-macro modeling of heterogeneous materials. *Computational Mechanics* 27(1):37–48. [DOI].

Kouznetsova, V., M. Geers, and W. Brekelmans (2002). Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme. *International Journal for Numerical Methods in Engineering*. [DOI].

Larsson, F., K. Runesson, S. Saroukhani, and R. Vafadari (2011). Computational homogenization based on a weak format of micro-periodicity for RVE-problems. *Computer Methods in Applied Mechanics and Engineering* 200(1-4):11–26. [DOI].

Linder, C. and F. Armero (2007). Finite elements with embedded strong discontinuities for the modeling of failure in solids. *International Journal for Numerical Methods in Engineering* 72(12):1391–1433. [DOI].

Liu, Y., F. van der Meer, L. J. Sluys, and L. Ke (2021). Modeling of dynamic mode I crack growth in glass fiber-reinforced polymer composites: fracture energy and failure mechanism. *Engineering Fracture Mechanics* 243:107522. [DOI], [OA].

Massart, T., R. Peerlings, and M. Geers (2007). An enhanced multi-scale approach for masonry wall computations with localization of damage. *International Journal for Numerical Methods in Engineering* 69(5):1022–1059. [DOI].

Matouš, K., M. Geers, V. Kouznetsova, and A. Gillman (2017). A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. *Journal of Computational Physics*
330:192-220. [DOI], [OA].

van der Meer, F. and L. J. Sluys (2009). A phantom node formulation with mixed mode cohesive law for splitting in laminates. *International Journal of Fracture* 158(2):107-124. [DOI], [OA].

van der Meer, F., L. J. Sluys, S. Hallett, and M. Wisnom (2012). Computational modeling of complex failure mechanisms in laminates. *Journal of Composite Materials* 46(5):603-623. [DOI].

van der Meer, F. (2012). Mesolevel modeling of failure in composite laminates: constitutive, kinematic and algorithmic aspects. *Archives of Computational Methods in Engineering* 19(3):381-425. [DOI].

van der Meer, F. (2016). Micromechanical validation of a mesomodel for plasticity in composites. *European Journal of Mechanics-A/Solids* 60:58-69. [DOI], [OA].

Melro, A., P. Camanho, F. Andrade Pires, and S. Pinho (2013). Micromechanical analysis of polymer composites reinforced by unidirectional fibres: Part I-Constitutive modelling. *International Journal of Solids and Structures*. [DOI], [OA].

Mercatoris, B., P. Bouillard, and T. Massart (2009b). Multi-scale detection of failure in planar masonry thin shells using computational homogenisation. *Engineering Fracture Mechanics* 76(4):479-499. [DOI].

Mercatoris, B. and T. Massart (2009a). Assessment of periodic homogenization-based multiscale computational schemes for quasi-brittle structural failure. *International Journal for Multiscale Computational Engineering* 7(2):153-170. [DOI], [HAL].

Miehe, C. and C. Bayreuther (2007). On multiscale FE analyses of heterogeneous structures: From homogenization to multigrid solvers. *International Journal for Numerical Methods in Engineering*. [DOI].

Moës, N. and T. Belytschko (2002). Extended finite element method for cohesive crack growth. *Engineering Fracture Mechanics* 69(7):813-833. [DOI], [HAL].

Neilsen, M. and H. Schreyer (1993). Bifurcations in elastic-plastic materials. *International Journal of Solids and Structures* 30(4):521-544. [DOI].

Nguyen, V. P., O. Lloberas-Valls, M. Stroeven, and L. J. Sluys (2012b). Computational homogenization for multiscale crack modeling. *Implementational and computational aspects. International Journal for Numerical Methods in Engineering* 89(2):192-226. [DOI], [OA].

Nguyen, V. P., M. Stroeven, and L. J. Sluys (2011). Multiscale Continuous and Discontinuous Modeling of Heterogeneous Materials: a Review on Recent Developments. *Journal of Multiscale Modelling* 03(04):229-270. [DOI].

Nguyen, V. P., M. Stroeven, and L. J. Sluys (2012a). An enhanced continuous-discontinuous multiscale method for modeling mode-I cohesive failure in random heterogeneous quasi-brittle materials. *Engineering Fracture Mechanics* 79:78-102. [DOI].

Oliver, J., M. Caicedo, E. Roubin, A. Huespe, and J. Hernández (2015). Continuum approach to computational multiscale modeling of propagating fracture. *Computer Methods in Applied Mechanics and Engineering* 294:384-427. [DOI], [OA].

Oliver, J., A. Huespe, J. Cante, and G. Díaz (2010). On the numerical resolution of the discontinuous material bifurcation problem. *International Journal for Numerical Methods in Engineering* 83(6):786-804. [DOI].

Rudnicki, J. and J. Rice (1975). Conditions for the localization of deformation in pressure-sensitive dilatant materials. *Journal of the Mechanics and Physics of Solids* 23(6):371-394. [DOI].

Simo, J., J. Oliver, and F. Armero (1993). An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. *Computational Mechanics* 12(5):277-296. [DOI].

van der Sluis, O., P. Schreurs, W. Brekelmans, and H. Meijer (2000). Overall behaviour of heterogeneous elastoviscoplastic materials: Effect of microstructural modelling. *Mechanics of Materials* 32(8):449-462. [DOI].

Song, J.-H., P. Areias, and T. Belytschko (2006). A method for dynamic crack and shear band propagation with phantom nodes. *International Journal for Numerical Methods in Engineering* 67(6):868-893. [DOI].

Souza, F. and D. Allen (2011). Modeling the transition of microcracks into macrocracks in heterogeneous viscoelastic media using a two-way coupled multiscale model. *International Journal of Solids and Structures* 48(22-23):3160-3175. [DOI], [OA].
Svenning, E., M. Fagerström, and F. Larsson (2017a). Localization aligned weakly periodic boundary conditions. *International Journal for Numerical Methods in Engineering* 111(5):493–500. [DOI], [OA].

Svenning, E., F. Larsson, and M. Fagerström (2017b). Two-scale modeling of fracturing solids using a smeared macro-to-micro discontinuity transition. *Computational Mechanics* 60(4):627–641. [DOI], [OA].

Svenning, E., F. Larsson, and M. Fagerström (2019). A two-scale modeling framework for strain localization in solids: XFEM procedures and computational aspects. *Computers & Structures* 211:43–54. [DOI].

Sánchez, P., P. Blanco, A. Huespe, and R. Feijóo (2013). Failure-Oriented Multi-scale Variational Formulation: Micro-structures with nucleation and evolution of softening bands. *Computer Methods in Applied Mechanics and Engineering* 257:221–247. [DOI].

Tadmor, E., M. Ortiz, and R. Phillips (1996). Quasicontinuum analysis of defects in solids. *Philosophical Magazine A: Physics of Condensed Matter, Structure, Defects and Mechanical Properties*. [DOI].

Toro, S., P. Sánchez, P. Blanco, E. De Souza Neto, A. Huespe, and R. Feijóo (2016a). Multiscale formulation for material failure accounting for cohesive cracks at the macro and micro scales. *International Journal of Plasticity* 76:75–110. [DOI], [OA].

Toro, S., P. Sánchez, J. Podestá, P. Blanco, A. Huespe, and R. Feijóo (2016b). Cohesive surface model for fracture based on a two-scale formulation: computational implementation aspects. *Computational Mechanics* 58(4):549–585. [DOI], [OA].

Turteltaub, S. and R. Suárez-Millán (2020). Energetically-consistent multiscale analysis of fracture in composites materials. *European Journal of Mechanics-A/Solids* 84:104079. [DOI], [OA].

Turteltaub, S., N. van Hoorn, W. Westbroek, and C. Hirsch (2018). Multiscale analysis of mixed-mode fracture and effective traction-separation relations for composite materials. *Journal of the Mechanics and Physics of Solids* 117:88–109. [DOI], [OA].

Verhoosel, C., J. Remmers, and M. Gutiérrez (2009). A dissipation-based arc-length method for robust simulation of brittle and ductile failure. *International Journal for Numerical Methods in Engineering* 77(9):1290–1321. [DOI].

Verhoosel, C., J. Remmers, M. Gutiérrez, and R. de Borst (2010). Computational homogenization for adhesive and cohesive failure in quasi-brittle solids. *International Journal for Numerical Methods in Engineering* 83(8-9):1155–1179. [DOI], [OA].

Wells, G. and L. J. Sluys (2001). A new method for modelling cohesive cracks using finite elements. *International Journal for Numerical Methods in Engineering* 50(12):2667–2682. [DOI].

Wu, J. Y. (2018). A geometrically regularized gradient-damage model with energetic equivalence. *Computer Methods in Applied Mechanics and Engineering* 328:612–637. [DOI].
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

<figure><img src="image_30.png" /></figure>

**Authors' contributions** L.K. carried out most of the study, performed numerical simulations, and drafted the manuscript. F.v.d.M. helped with implementation and numerical issues. Both authors developed the methodology, conceived of the study, and participated in its design, coordination, and critical review of the manuscript. Both authors read and approved the final manuscript.

**Supplementary Material** The original datasets used in the simulations of the present article can be found at the **permalink 10.5281/zenodo.6044934**. The shared material contains the inputs and outputs of the simulations. Detailed input parameters for each test can be found in the **.pro** file and the corresponding mesh in the **.msh** and **.geo** files. The Matlab script **postpro.m** post-processes the raw output files **.cracks** and **.dat**.

**Acknowledgements** Financial support from the Dutch Research Council NWO under Vidi grant 16464 is gratefully acknowledged.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.