Two-scale modeling of fracturing solids using a smeared macro-to-micro discontinuity transition

We consider multiscale modeling of fracturing solids undergoing strain localization, whereby Statistical Volume Elements (SVEs) are used to compute the homogenized macroscopic stresses and the eXtended Finite Element Method (XFEM) is used to represent macroscale displacement discontinuities. These discontinuities are imposed on the localized SVEs in a smeared sense, whereby the smearing width is related to the SVE size and the orientation of the macroscopic discontinuity. This smearing width relation, which is derived within the setting of Variationally Consistent Homogenization (VCH), prevents pathological dependence of the solution on the SVE size. The SVE size insensitivity is further improved by adopting the recently proposed localization aligned weakly periodic boundary conditions. Advantages of the proposed method are that it allows multiscale modeling of localized fracture without restrictive assumptions on the SVE size and without the need to explicitly track a localized region in the SVE.


Introduction
Multiscale modeling of fracturing solids is a challenging topic that has been the subject of considerable research efforts, cf. the reviews in [1,2]. A frequently used approach B Erik Svenning erik.svenning@chalmers.se 1 Division of Material and Computational Mechanics, Department of Applied Mechanics, Chalmers University of Technology, Gothenburg, Sweden is to perform numerical simulations on Statistical Volume Elements (SVEs), 1 whereby a fundamental challenge is that the SVE looses its representative character upon localization. This loss of representativeness leads, for standard first order computational homogenization, to pathological dependence of the numerical results on the macroscale mesh size and the SVE size. These problems have been addressed in the literature using several different schemes. A homogenization scheme that employs discrete cracks on both scales was proposed by Verhoosel et al. [5], where a homogenized traction-separation law is derived, assuming elastic material behavior prior to crack nucleation. An alternative scheme that assumes a weak discontinuity on the microscale was proposed by Coenen et al. [6]. Souza and Allen [7] instead considered visco-elastic materials, and developed a homogenized traction-separation law using interface elements on the microscale and XFEM on the macroscale. Methods employing XFEM were also proposed by Belytschko et al. [8] and Bosco et al. [9], whereas element embedded discontinuities were employed by Toro et al. [10].
The problem with pathological SVE size dependence can be addressed by separating the bulk deformation and the discontinuous deformation in a suitable way, thereby allowing incorporation of a length scale related to the fracture zone. Souza and Allen [7] compute the bulk and interface response from different SVEs, imposing only the discontinuous part of the strain on the localized SVE. Belytschko et al. [8] instead introduce a perforated unit cell, and separate the bulk deformation from the discontinuity by exclud- 1 In the literature, both Representative Volume Element (RVE), Statistical Volume Element (SVE) and Microstructural Volume Element [3] are used to denote a sample of the microstructure. To stress the fact that a sample of finite size will, in general, not be truly representative, we prefer the notion Statistical Volume Element (SVE), cf. Ostoja-Starzewski [4]. ing unstable subdomains. This method allows simulations involving macroscale localization, but requires that the SVE size is approximately equal to the macroscale element size. A scheme that allows wider scale separation was proposed in [3,6,9], where the effective microscale discontinuity is identified by means of image analysis techniques in order to separate the bulk and discontinuity parts of the deformation. In [11], the effective discontinuity was instead identified by monitoring the damage evolution in the SVE. Even though these schemes have successfully addressed several difficulties associated to multiscale localization, we note that previous work tends to focus on approaches where the topology of the microscale discontinuity is explicitly identified and tracked dynamically in different ways [3,6,[8][9][10][11]. Therefore, a simpler but more versatile approach, that does not require dynamic tracking of the microscale discontinuity, would be very valuable.
Regarding the choice of suitable boundary conditions (BCs) on the SVE, it is well known that both Dirichlet, strong periodic and Neumann BCs are inaccurate if localization occurs in the SVE. In fact, these BCs are inaccurate also during the early stages of crack propagation, i.e. prior to localization, if cracks intersect the SVE boundary, see e.g. the comparison between different BCs in [12]. Remedies proposed in the literature are to mix Dirichlet and Neumann BCs on different parts of the SVE boundary [8], to use Percolation Path Aligned BCs on strong form [3,6], or to use localization aligned periodic BCs on weak form [13].
In the present work, we choose to use localization aligned weakly periodic BCs on the SVE and adopt the concept of Variationally Consistent Homogenization (VCH) to develop a multiscale method that is computationally efficient when the macroscale deformation localizes in a narrow band. We allow for the use of any suitable fracture model (XFEM, interface elements, embedded discontinuities, nonlocal continuum damage) on the microscale, but we consider strong discontinuities on the macroscale, and impose these discontinuities on the SVE in a smeared sense, thereby avoiding the need to explicitly identify and track a localization region in the SVE. More precisely, the homogenized response obtained from the SVEs in the damaged zone provides an (implicit) cohesive zone law, relating the macro discontinuity to the effective macro traction computed from the SVE response. Using a smeared macro-to-micro discontinuity transition to impose the macro discontinuity on the SVE, the macro discontinuity is thereby implicitly driven by the microscale response.
The remainder of the paper is organized as follows: The strong and weak forms of the single scale model problem under consideration are established in Sect. 2, followed by the derivation of the two-scale formulation in Sect. 3, a few illustrative examples in Sect. 4, and some concluding remarks in Sect. 5.

Model problem
In this section, we will establish the fully resolved problem, i.e., the single scale problem prior to introduction of computational homogenization. To this end, let us consider a solid domain Ω with external boundary Γ ext and internal surfaces Γ int representing cracks in the material as indicated in Fig. 1. The internal boundaries Γ int have a predefined normal n int and consist of two-sided surfaces, as also indicated in Fig. 1. The external boundary consists of a part Γ ext,D with Dirichlet boundary conditions and a part Γ ext,N with Neumann boundary conditions, so that the boundary of Ω is decomposed as ∂Ω = Γ ext,D ∪ Γ ext,N ∪ Γ int . Letting superscripts + and − denote quantities on the positive side and the negative side of the internal boundaries Γ int , respectively, we define the normal as n int def = n int − . Since the cracks represented by Γ int may branch and intersect, n int is not necessarily continuous along a chosen part of Γ int . Considering small strains and quasistatic loading, and neglecting the body force, the strong form of the equilibrium equations is given by where σ = σ { } is the Cauchy stress, = ([u ⊗ ∇] sym ) is the strain, ∇ is the gradient operator, n is a unit normal vector, t is a prescribed traction andû is a prescribed displacement. The traction on the internal boundaries Γ int is described by a cohesive zone law specifying t + = −t − = t u in terms of the jump u def = u + − u − over the crack faces. Note that in the constitutive relations for the stress-strain relation σ = σ { } and the traction-separation law t = t u , the dependence of σ and t on internal variables has been omitted for brevity. Even though body forces have been neglected in Eq. (1), such forces can be included without fundamental difficulty.
The (displacement based, single field) weak solution corresponding to Eq. (1) is obtained by finding u ∈ U such that where d is the number of spatial dimensions in the problem, t u denotes the traction obtained from the cohesive zone constitutive law, and where Regarding the regularity requirements in U and U 0 , H 1 (Ω) denotes the space of square integrable functions with square integrable derivatives, and we recall that Ω is the "perforated" solid domain.
We remark that Eq. (2) represents the fully resolved problem in the sense that computational homogenization has not yet been introduced. In the next section, we will consider the situation that the fully resolved problem is intractable due to vast separation of length scales and derive a two-scale formulation of this equation.

Preliminaries
In this section, we aim to construct a two-scale model of the problem given by Eq. (2), where Γ int defines a large set of microcracks. On the macroscopic level, we consider the con-  Fig. 2 Macroscopic domainΩ divided into a partΩ d containing an effective macroscopic discontinuity surfaceΓ d , and a partΩ r that is free from macroscopic discontinuities. The width l d of the discontinuity regionΩ d is also shown part of the continuous bulk domainΩ where the microcracks are separated and the effective (macroscale) response shows no localization. In the discontinuity regionΩ d , on the other hand, we assume that microcracks have coalesced into macrocracks 3 causing a macroscale localization. We assume that the localization zoneΩ d is narrow, wherebyΩ d can be described accurately by the effective discontinuity surfacē Γ d and a width l d , cf. Fig. 2.

Variationally consistent homogenization for continua with macroscopic localization
To derive a two-scale model, we adopt the procedure of Variationally Consistent Homogenization (VCH), which was introduced by Larsson and Runesson [14] in the context of adaptive scale bridging, and later applied to transient heat flow in [15] and to microfracturing solids in [12]. The key ingredient in the present work, which is an extension of [12], is the introduction of macroscopic localization. As the first step of VCH, we introduce the decomposition U = U M ⊕ U s and consider a split of the solution into macro and micro parts according to u = u M + u s , where u M ∈ U M describes smooth (average) parts of the solution and u s ∈ U s are the microscale fluctuations. Assuming the same split for any function in U 0 , we can state the original problem in Eq.
(2) as the macro problem (6) and the micro problem where we recall that a(•, •), b(•, •) and l(•) were defined in Eqs. (3), (4) and (5), respectively. So far, we note that Eqs. (6) and (7) pertain to the Variational MultiScale method (VMS) introduced by Hughes et al. [16]. The second step of VCH is to restate the integrals in Eq. (6) using running averages in order to obtain a homogenized problem. To this end, we consider the split of the domain Ω =Ω r ∪Ω d as discussed in Sect. 3.1 and, for the regular domainΩ r , we introduce the approximation where Ω denotes a Statistical Volume Element (SVE), formally unique for any point inΩ.
Remark Note that, if we split the entire domain into adjacent SVEs so thatΩ r = ∪ i Ω ,i , and assume the SVE integral f to be piecewise constant in each Ω ,i , the approximation in Eq. (8) becomes an identity. For the discontinuity regionΩ d , we state the integrals based on the interfaceΓ d according to where the last approximation is based on the assumption of a narrow regionΩ d , and where f , that was defined in Eq. (8), is now evaluated in a point onΓ d .
In the presence of localization, it is clear that the quantity f in Eq. (9) will not be representative if evaluated on an arbitrarily sized SVE Ω . To overcome this issue, we con-siderΩ d entirely populated by SVEs through the (narrow) width l d , 4 cf. Fig. 2. Defining an effective discontinuity sur-faceΓ d, as the plane through the center of the SVE parallel toΓ d , we then obtain the geometric relation With this choice, it becomes evident that the homogenization in Eq. (9) is a homogenization onΓ d , whereas the direction perpendicular toΓ d is fully resolved. We also note that, since l d depends only on the SVE size and the orientation of the effective discontinuityΓ d , l d is independent of the loading direction.
Finally, in order to derive the homogenized problem, prolongation conditions need to be specified. More precisely, the macroscale part u M of the resolved solution inside Ω needs to be expressed in terms of the homogenized solution fieldū. To this end, we shall adopt different strategies forΩ r andΩ d . We start off with introducing the homogenized displacement ansatzū ∈Ū which is smooth onΩ def =Ω r ∪Ω d \Γ d but may be discontinuous acrossΓ d . For SVEs pertaining to the regular subdomainΩ r , we use the ansatz from conventional first order homogenization, where Ω (x) is the SVE centered aroundx. For SVEs located onΓ d , we adopt a smeared approach, whereby the average strain inside Ω is expressed as the sum of the average bulk strain and the strain contribution from the displacement jump smeared over l d . Hence, we state the macroscale displacement field as Here, we define the effective strain onΓ d as where we introduced the bulk strain onΓ d as 5 and where n is the normal ofΓ d atx. It is easily verified that the integration of Eq. (12) gives displacement compatibility betweenū and u M on the boundary ofΩ d . For the numerical implementation, we further introduce the approximation where the average of the strains replaces the integral expressions.
Using the integral transformations in Eqs. (8) and (9), together with the prolongation expressions in Eqs. (11) and (12), we may now restate Eq. (6) as the macroscale problem of findingū ∈Ū such that where¯ def = [ū ⊗ ∇] sym and the effective macroscale stress σ is obtained from the microscale solution according tō

Operational format of integral expressions
Solving the macroscale problem as given by Eq. (17) requires evaluation of an integral overΩ r . This evaluation is inconvenient from a computational perspective, because the finite element mesh will not matchΩ r . In order to restate is the discontinuity region excluding the sharp discontinuity surface. Then, we may employ the approximation where the discontinuity is removed from the argument (i.e. 0 and δ¯ 0 ) to account for the fact thatΩ d does not contain where we also used that Inserting the result from Eq. (20) in Eq. (17), the macroscale problem is then to findū ∈Ū such that wherē We note that I 1 represents the contribution from the effective bulk response, and that I 3 represents a contribution of cohesive zone type (homogenized from a damaging SVE) sincē Note, however, that the cohesive tractiont d depends also on the bulk strain becauseσ Furthermore, we remark that the novel term I 2 is scaled by l d , and therefore will be negligible for sufficiently small l d . To ensure energy consistency between the micro-and macroscales, the Hill-Mandel condition needs to be fulfilled. Following Larsson et al. [20] it is clear that a Hill-Mandel condition on the form is met. Utilizing the relation¯   6 With small, we here mean that l is small compared to the macroscale dimensions. Note that the SVE size cannot be chosen arbitrarily small: the SVE still needs to be sufficiently large to give a good statistical representation of the microstructure.
-Alt. II: l d is assumed to be small, so that -Alt. III: As Alt. II, but the bulk strain in the SVEs onΓ d is neglected:

Localization aligned weakly periodic boundary conditions
To solve the microscale problem obtained from Eq. (7) on an SVE Ω , the macroscopic part of the displacement given by Eq. (11) or Eq. (12) is imposed on the SVE boundary through suitable BCs. In the present work, we will use weakly periodic boundary conditions that are aligned to the direction of the effective discontinuity in the SVE [13]. To this end, we divide the SVE boundary into an image part Γ + and a mirror part Γ − as shown in Fig. 4a. Furthermore, we introduce a mapping ϕ per : Γ + → Γ − such that points on Γ + and Γ − are associated to each other according to x − = ϕ per (x + ).
For SVEs pertaining to bulk response, i.e. for quadrature points inΩ, the mapping is typically performed along horizontal and vertical lines between the surfaces of the square SVE domain. However, for SVEs on the localization band Γ d , the mapping is constructed such that the periodicity directions are aligned with the direction of the effective discontinuity, cf. Fig. 4a. We also define the jump 7 between a point x + on Γ + and the associated point x − = ϕ per (x + ) on Next, we impose localization aligned weakly periodic boundary conditions on the SVE by introducing an independent discretization for the boundary traction t λ and requiring u =¯ · x −x to hold in a weak sense on Γ + . The SVE problem is then to find u ∈ U and t λ ∈ T such that where we introduced the expressions and where L 2 Γ + denotes the space of square integrable functions on Γ + . See [13,20,21] for further details. We remark that, onΓ d , the weakly periodic boundary conditions are employed to impose the macroscopic strain d =¯ 0 + 1 l d u ⊗ n sym in a weak sense on the whole SVE boundary. Hence, there is no need to explicitly identify the location or width of the damaged zone in the SVE, only the localization direction needs to be determined in order to define the effective discontinuity normal direction n. Recall that n also defines the periodicity alignment. We also remark that l d cannot be chosen freely, it is related to the SVE size |Ω | according to Eq. (10). For the special case of 2D and a square SVE, we may compute l d explicitly as where l = √ |Ω | and α is the angle between the effective discontinuity and the x-axis.

Preliminaries
In this section, we demonstrate the performance of the proposed method with a few numerical examples. To assess the accuracy of the proposed method, we will present comparisons between fully resolved Direct Numerical Simulations (DNS) and two-scale simulations (FE 2 ). The open source software package OOFEM [22,23] has been used for the numerical implementation.
In Example 4.2, we present a comparison between Alt. I, Alt. II and Alt. III described in Sect. 3. In Examples 4.3 and 4.4, we use Alt. I and demonstrate that the proposed method accurately predicts homogenized stress-strain relations, provided that an appropriate kinematic ansatz is made for the macroscale displacement field. Finally, an outlook is presented in Example 4.5, indicating the application of the proposed method to cases where the macroscale localization pattern is not known a-priori. Here, Alt. III is used to reduce the computational cost.

Elastic specimen with vertical localization band
To investigate the SVE size dependence for different approximations to Eq. (21), we consider the specimen shown in Fig. 5. The specimen has length L = 10 mm, height H = 5 mm, notch diameter D = 0.5 mm and a thickness of 1 mm. The bulk material is isotropic and linear elastic, with Young's modulus E = 210 × 10 3 MPa and Poisson's ratio ν = 0.3. Plane strain is assumed. The interface is modeled by an isotropic elastic cohesive zone model with stiffness k = 1.0 × 10 5 N mm −3 .
The left edge of the specimen is clamped and the right edge is subjected to a prescribed displacement according to u x = 0.01y/H, u y = 0, see Fig. 6.
The response of the specimen is computed in two ways: i) with DNS and ii) with FE 2 for the interface response. For the DNS, we use a standard discretization with 6-node triangles in the bulk and model the interface by means of XFEM using the previously mentioned isotropic elastic interface law. 8 For the FE 2 , we also use 6-node triangles with isotropic linear elastic material model (i.e. no homogenization of bulk response in this example) and an XFEM discretization of the interface. However, the interface response is now computed from SVEs in the cohesive zone Gauss points.
As output from the simulations, we choose to monitor the normal traction along the cohesive interface (Fig. 7). The normal tractions predicted with Alt. I and Alt. II are shown in Fig. 7a, b, respectively. As can be seen, these approximations both agree very well with the DNS solution. Alt. III, shown in Fig. 7c, agrees very well with the DNS for the smallest SVE, but large discrepancies occur for the largest SVE. This is expected, because Alt. III neglects the bulk strain in the SVE, and the bulk strain is more important for large SVEs.
In summary, all three approximations are accurate for the smallest SVE (i.e. when strong scale separation holds), whereas Alt. I and Alt. II perform better than Alt. III when the SVE size approaches the size of the macroelements.

Plate with inclined effective discontinuity surface
To investigate how the orientation of the effective discontinuity surface influences the macroscopic response, we consider a square plate with a circular hole and a band of porous material as shown in Fig. 8. The orientation of this softer band is described by the angle α between the band and the x-axis. We consider an isotropic and linear elastic bulk material in plane strain with Young's modulus E = 210 × 10 3 MPa and Poisson's ratio ν = 0.3. The material is homogeneous, except in the softer band, where small holes in the material are present as shown in Figs. 9 and 10. The side length of the plate is L = 100 mm and the diameter of the macroscopic hole is d = 60 mm. The small holes in the soft band have a diameter of 1 mm, and the spacing between the small holes (centercenter distance) is 1.25 mm. Three rows of small holes run across the entire plate as shown in Figs. 9 and 10.
We solve the problem using both DNS and FE 2 for different values of the angle α. DNS meshes for three different angles are shown in Fig. 9. For the DNS meshes, 3-node triangles with a side length of 1 mm are used, except in the region surrounding the small holes, where a side length of 0.025 mm is used. For the FE 2 simulations, 6-node triangles with an element side length of 5.0 mm is used on the macroscale. For the SVE meshes, 6-node triangles with a side length of 0.5 mm are used, except in the region surrounding the small holes, where a side length of 0.1 mm is used. Regarding the SVE size, we choose l = 20 mm.
The left edge of the specimen is clamped, and a uniform displacement of u 0 = 0.01 mm is applied in the x-direction on the right edge. As output from the simulations, we choose to monitor the reaction force in the x-direction.
The deformed shape of the structure, computed with DNS for α = 60 • , is shown in Fig. 11. As can be seen, a substantial amount of the deformation is localized to the soft band. This is also reflected in the stress field shown in Fig. 12: the effective stress is much higher around the holes than in the rest of the structure.
The reaction force on the clamped edge, computed with both DNS and FE 2 , is shown for different α in Fig. 13. The agreement between the DNS and FE 2 simulations is good, the largest difference in reaction force between the DNS and the FE 2 simulations is 7% (for α = 45 • ). The larger discrepancy close to α = 45 • is most likely due to the finite width of the soft region. More precisely, the DNS will capture the effect of the finite width soft region overlapping the upper right corner, whereas the FE 2 macroscale representation of the soft region is a strong discontinuity with zero (macroscale) thickness. Despite the discrepancy around α = 45 • , we conclude that the agreement is good for all values of α.

Nonlinear response of a plate with a softening region
To verify that the response during softening is insensitive to the SVE size, we consider a plate with a softening region containing a cohesive surface and holes as shown in Fig. 14. The plate has side length L = 100 mm and thickness 1 mm (plane strain is assumed). The cohesive band has holes with radius 1 mm and center-center distance 4 mm, and runs through the center of the plate at an angle of α = 25 • to the horizontal axis.
The cohesive interface is modeled by a standard bilinear cohesive law. In this example, we consider monotonic loading and thus disregard history effects in the cohesive zone material, thereby employing an "elastic" cohesive law with softening as indicated in Fig. 15. For the cohesive zone material parameters, we choose the initial stiffness as k = 1.0 × 10 5 N mm −3 , the peak stress as σ f n = 2.5 × 10 2 MPa and the fracture energy as G I = 2.0 × 10 1 N mm −3 .
For the DNS, the bulk material is discretized with 6-node triangles with a side length of 1 mm, except close to the interface, where the side length is 0.075 mm. The cohesive surface is modeled by interface elements. For the FE 2 simulations, XFEM is used on the macroscale, whereas interface elements are used to represent the microscale discontinuity. Regarding the bulk discretization, 6-node triangles with a side length of 1.25 mm are used on the macroscale. For the microscale, we use 6-node triangles with a side length of 0.36 mm, except close to the interface where the side length is reduced to 0.18 mm. Two SVE sizes are considered: l = 3.62 mm and l = 7.24 mm, cf. Fig. 16.
A prescribed displacement u 0 as indicated in Fig. 14 is added incrementally in 10 load steps. As output from the simulations, we choose to monitor the reaction force for the lower right support. The reaction force vs displacement curve for the lower right support, computed with DNS and FE 2 , is shown in Fig. 17. As can be seen, the FE 2 solution is insensitive to the SVE size and agrees very well with the DNS solution.

Plate with localizing microcracks
In this example, we consider a plate with a circular hole and a microstructure consisting of a matrix with circular inclusions as shown in Fig. 18. The plate, which has length L = 50 mm and hole diameter D = 20 mm, is clamped at the left edge. The right edge is fixed in the horizontal direction and a uniform displacement is applied in the vertical direction.
The matrix material is isotropic and linear elastic with Young's modulus E m = 210 × 10 3 MPa and Poisson's ratio ν m = 0.3. The inclusions are also isotropic and linear elastic, with Young's modulus E i = 10.0 × 10 3 MPa and Poisson's ratio ν i = 0.3. The inclusions have diameter d = 0.5 mm and are randomly distributed in the plate. For the microscale problem in the FE 2 simulations and the DNS, cracks nucleate if the highest principal stress exceeds 400 MPa. Crack propagation is modeled using the concept of material forces, whereby the cracks propagate in the direction of the material force when the magnitude of the material force exceeds 0.25 N. On the macroscale, we consider two different meshes as shown in Fig. 19. For the macroscale crack initiation, we choose for simplicity 9 to nucleate cracks when the highest principal strain exceeds 1 × 10 3 in any integration point. A macroscopic crack segment is thereby introduced in the direction perpendicular to the corresponding principal strain direction. Such macroscopic new segment is also allowed to intersect already existing cracks on the macroscale. Existing macroscale cracks propagate in the direction perpendicular to the highest principal strain when the highest principal strain in the integration point closest to the crack tip exceeds 1 × 10 −4 .
The resulting effective stress predicted with DNS is shown in Fig. 20, and the displacement field, computed with DNS and FE 2 , is shown in Fig. 21. In the DNS, the propagating microcracks interact with the microstructure to form an irregularly shaped macrocrack. The FE 2 simulation does not capture the small scale irregularities in the crack pattern on the macroscale, but the overall crack pattern is in qualitative agreement with the DNS. The reaction forces computed with DNS and FE 2 are shown in Fig. 22. Despite the crude macroscale crack initiation criterion adopted in the present example, the FE 2 solution agrees fairly well with the DNS.

Conclusions
In the present work, we model fracturing solids undergoing strain localization using a two-scale approach based on Variationally Consistent Homogenization (VCH). We propose a continuous-discontinuous homogenization scheme, where the eXtended Finite Element Method (XFEM) is used to represent narrow macroscale localization bands. The displacement discontinuities across the macroscopic discontinuity surfaces are transferred to the microscale in a smeared sense, using a smearing width l d that is related to the orientation of the effective discontinuity and the SVE size. Using a correct smearing width turns out to be crucial for obtaining accurate results without pathological dependence on the SVE size. Adopting the procedure of VCH and assuming that the macroscale localization region is well approximated by a strong discontinuity, we derive a two-scale scheme containing a conventional bulk term, a term of cohesive zone type and a correction term. The latter term has negligible effect on the accuracy if strong scale separation holds, but is indeed needed to retain any (possible) symmetry of the problem. A major benefit of the proposed scheme is that it does not require dynamic tracking of an evolving damage region on the microscale, and that it does not explicitly assume a particular constitutive behavior on the microscale. In the numerical examples, we show that the proposed scheme is insensitive to the SVE size if all terms obtained in the derivation are included, and that approximations are possible for sufficiently small SVE sizes l (i.e. in the case of strong scale separation). Furthermore, it is demonstrated that the proposed scheme, in combination with localization aligned weakly periodic BCs, is accurate for varying crack orientations. We also illustrate the potential of the proposed scheme with an FE 2 simulation of a localizing microstructure. In the last example presented, the potential of the proposed method has been demonstrated using sim-  plified models for crack initiation and propagation. Future work therefore involves more accurate models for macroscale crack propagation, i.e. a proper strategy for insertion and propagation of (possible) macroscale discontinuities. In summary, a multiscale modeling scheme that allows microscopic and macroscopic localization is proposed. The good performance of the proposed scheme is demonstrated by several numerical examples.