Multi-scale continuous–discontinuous framework for computational-homogenization–localization

Journal of the Mechanics and Physics of Solids 60 (2012) 1486–1507

Contents lists available at SciVerse ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

Multi-scale continuous–discontinuous framework for computational-homogenization–localization E.W.C. Coenen a,b, V.G. Kouznetsova b,n, M.G.D. Geers b a

Materials innovation institute (M2i), P.O. Box 5008, 2600 GA Delft, The Netherlands Eindhoven University of Technology, Department of Mechanical Engineering, Section of Materials Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands

b

a r t i c l e in f o

abstract

Article history: Received 22 March 2011 Received in revised form 15 March 2012 Accepted 8 April 2012 Available online 16 April 2012

In this paper, a multi-scale technique is proposed for the modeling of microstructured materials up to the point of macroscopic failure. A continuous–discontinuous computational homogenization–localization framework is developed, which involves a discontinuity enhanced macroscale description. The underlying microstructural volume element (MVE) enables the incorporation of a band with high strains, i.e. a localization band. For the multi-scale coupling, special scale transition relations are established to handle the underlying material response of both the bulk and the localization zone. At a macroscopic (integration) point, the macroscale displacement jump and the deformation of the surrounding continuum material are used to formulate kinematical boundary constraints for the microscale MVE problem. Upon proper averaging of the MVE response, the macroscopic generalized stress is obtained. Simultaneously, a microscale effective displacement jump is recovered and returned to the macrolevel, based on the MVE deformation field. The equality of the macro- and effective microscale displacement jumps is enforced at the macroscale by Lagrange multipliers, being cohesive tractions at the interface. The applicability of the developed continuous–discontinuous computational homogenization framework is illustrated on two simple benchmark problems involving evolving macroscale discontinuities as a result of microstructural degradation by void growth. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Computational homogenization Multi-scale Localization Damage Cohesive crack

1. Introduction In the last decades, multi-scale models are developed and used in many application fields, e.g. advanced alloys, biological materials, functional and composite materials. In many cases, the key question is to establish a structure– property relation: how a particular microstructure governs the overall or effective material properties. At the same time, microstructural modeling can be used to assess material failure aspects and to investigate material strength and formability. Material failure at the macroscale results from a cascade of events at the microstructural level. Due to these events, eventually a zone with high strains, which crosses the microstructure, will form, e.g. by void growth and plastic strain localization in the ligaments between the voids, as is the case in the ductile damage of structural metals and metal matrix

n Correspondence to: Department of Mechanical Engineering, Eindhoven University of Technology, Den Dolech 2, 5612 AZ Eindhoven, The Netherlands. Tel.: þ 31 40 2475885; fax: þ31 40 2447355. E-mail address: [email protected] (V.G. Kouznetsova).

0022-5096/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jmps.2012.04.002

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1487

composites. This strain localization band leads to strain softening characterized by a decreasing load carrying capacity of the microstructure across the band. At the larger, macroscale, this phenomenon can be modeled by a cohesive crack, i.e. a discontinuity of zero width that can transmit a force, diminishing (eventually to zero) with increasing discontinuity opening. The initiation and development of a strain localization band for complex materials strongly depend on the loading path, non-linear material behavior of the microstructural constituents, and the microstructural evolution. In order to incorporate the micromechanical response in large scale structural analysis simulations at the engineering level, the mechanics at different length scales needs to be related. Multi-scale methods are a collection of different strategies to establish this coupling across the scales. Domain decomposition methods decompose the large scale domain in generally non-overlapping subdomains that can have different spatial resolutions (Zohdi and Wriggers, 1999; Farhat et al., 2000; Guidault et al., 2008). This makes it possible to resolve a strain localization band in the fine resolution domain and save computational effort in zones of less interest. A second class of multi-scale techniques are superposition (or multigrid) based methods that make use of a hierarchical decomposition of macro- and microscale effects (Fish et al., 1997; Hund and Ramm, 2007; Loehnert and Belytschko, 2007; Gitman et al., 2008). In the zone of interest the coarse scale discretization is complemented by a fine-scale overlay to enhance the resolution of the solution. In general, the zone of interest is not always known in advance and adaptive techniques may have difficulties incorporating the deformation history upon mesh refinement. Moreover, these techniques are limited to moderate separation of scales. In contrast to local mesh refinement techniques, multi-scale homogenization techniques rely on a separation of scales principle. In this class, computational homogenization approaches (also called FE2 in a particular form) are suitable to deal with non-linear and evolving microstructures. Computational homogenization is essentially based on the solution of two (nested) boundary value problems (BVPs), one for the macroscopic and one for the microscopic scales (Suquet, 1985; Ghosh et al., 1995; Feyel and Chaboche, 2000; Kouznetsova et al., 2001; Miehe et al., 2010). Though computationally expensive, the procedures developed allow to assess the macroscopic influence of microstructural parameters in a rather straightforward manner. Techniques of this type (i) do not require any constitutive assumption on the macrolevel; (ii) enable the incorporation of large deformations and rotations on both micro and macrolevels; (iii) are suitable for arbitrary material behavior, including physically non-linear and history dependent behavior; (iv) provide the possibility to introduce detailed microstructural information, including a physical and/or geometrical evolution of the microstructure, in the macroscopic analysis; (v) allow the use of any modeling technique at the microlevel for the solution of the corresponding BVP. Classical computational homogenization schemes are based on the full separation of scales and are therefore consistent with the standard local continuum mechanics concept, where the response at a (macroscopic) material point depends on the deformation (i.e. the first gradient of the displacement field) of that material point only. The (numerical) stress–strain response is obtained as the volume averaged stress–strain response of a finite microstructural volume (MVE) that is explicitly modeled. In classical computational homogenization, the MVE is assumed to be (statistically) representative for the macroscopic material point and therefore is called a representative volume element (RVE). The development of a strain localization band crossing the MVE violates the separation of scales principle and the RVE therefore loses its representative character (i.e. terminating the existence of an RVE). Moreover, the standard continuum formulation at the macroscale is unable to resolve the strain localization band or (cohesive) crack. In addition, due to the strain softening response of the microstructure, the macroscale boundary value problem becomes ill-posed. Clearly, the strain localization inevitably limits the concept of classical homogenization. Second-order computational homogenization regularizes the problem by using a second gradient continuum at the macroscale (Kouznetsova et al., 2002, 2004). This resolves the ill-posedness of the macroscale problem and allows moderate strain localization in a macroscopic zone larger than the MVE size. However, strong localization within a band smaller than the MVE size cannot be resolved at the macroscale. Consequently, the second-order approach cannot be used beyond the onset of localization, leading to crack propagation and fracture. The enrichment of the macroscale continuum with a discontinuity enables lumping the microscale strain localization band into a macroscale cohesive crack. Several authors have proposed methods to homogenize the micromechanics of an MVE within an interface (adhesive) layer towards a cohesive relation (Larsson and Zhang, 2007; Matous et al., 2008; Hirschberger et al., 2009). Multi-scale schemes that enrich the macroscale continuum with a discontinuity have recently been proposed by Massart et al. (2007), Belytschko et al. (2008), Verhoosel et al. (2010), Mercatoris and Massart (2011), Souza and Allen (2011), and Nguyen et al. (2011, 2012). Different techniques are used to extract an equivalent discontinuity and the corresponding macroscopic effective stress–strain relation from the microstructural response. This is done by either assuming the continuum surrounding the localization zone to behave linear elastically (Massart et al., 2007; Verhoosel et al., 2010; Mercatoris and Massart, 2011; Nguyen et al., 2011) or by using separate microscale volume for the bulk homogenization (Belytschko et al., 2008; Souza and Allen, 2011; Nguyen et al., 2012). In this paper, a multi-scale scheme is proposed that departs from a classical computational homogenization and which evolves gradually towards an enhanced scheme. Upon detection of localization (i.e. the onset of strain softening of the homogenized microstructural response), the macroscale is enriched by a discontinuity representing the macroscopic displacement jump equivalent to the displacement jump extracted from the microstructural response. The most essential feature distinguishing the framework proposed here is that both the discontinuity and the stress–strain response of the surrounding bulk material are extracted from the same MVE simulation, thus circumventing the need for different MVEs for the crack and for the bulk. Simultaneously, the evolving character of the fully coupled micro–macro problem is

1488

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

captured, providing details at both scales. This makes the framework applicable to a broad class of materials including arbitrary non-linear, history dependent and evolving microstructures. The paper is organized as follows, in Section 2 the governing relations of the enhanced macroscale formulation are derived. Section 3 focuses on the detailed derivation of the required scale transition relations. The scheme is illustrated in Section 4 on two numerical examples. ! The following notation for Cartesian tensors and tensor products is used throughout the paper: a , A, and nA denote, respectively, a vector, a second-order tensor, and a nth-order tensor. The following notation for vector and tensor ! ! ! ! ! ! operations is employed: the dyadic product a  b ¼ ai bj e i  e j , and the inner products A  B ¼ Aij Bjk e i  e k , ! c A : B ¼ Aij Bji , with e i , i ¼ 1; 2,3 the unit vectors of a Cartesian basis; conjugation of a tensor Aij ¼ Aji and left-conjugation ALC ijkl ¼ Ajikl . The subscript ‘‘M’’ refers to a macroscopic quantity, while the subscript ‘‘m’’ will denote a microscopic quantity. 2. Multi-scale framework The discontinuity enhanced computational homogenization scheme presented here departs from the classical computational homogenization scheme and is identical to it prior to localization. The classical scheme resolves the local stress balance at the macroscale, i.e. ! r0M  PcM ¼ 0 ð1Þ Here PM is the first Piola–Kirchhoff stress tensor and r0M is the gradient operator with respect to the reference macroscopic configuration. The (numerical) constitutive stress–strain relation of a macroscopic material point is obtained by imposing the macroscopic deformation gradient tensor FM in an averaged sense onto a microstructural volume element (MVE), see Fig. 1. After the microstructural boundary value problem is solved, the macroscopic stress PM is obtained as the volume average of the MVE stress field. In this way, the stress–strain relation is obtained for every macroscopic material point. In practice, depending on the employed macroscale solution technique, the stress–strain response is only evaluated in a finite number of material points, e.g. the Gauss points in the case of finite element modeling. In addition to the homogenized stress–strain response, the consistent tangent operator is returned to the macroscale in order to iteratively solve the macroscale boundary value problem. The tangent operator 4 CM can be defined as the linear map between the infinitesimal increments in the stress and deformation tensors according to dPM ¼ 4 CM : dFcM

ð2Þ

As a result of microstructural heterogeneity, strain percolation paths (i.e. paths of strains higher than average) will form in the microstructure. Due to this, microstructural damage mechanisms may be activated, degrading the load carrying capacity of the microstructure and resulting in a decrease of the homogenized material stiffness. Ultimately, this strain softening leads to the loss of positive definiteness of the consistent tangent operator and the ill-posedness of the macroscale boundary value problem. Hill (1962) formulated the condition for the existence of a stationary discontinuity !d with normal N M in solids, see also Rice (1976) in the context of elasto-plastic materials. Such a stationary discontinuity can exist if the following condition is satisfied: !d !d detð N M  ð4 CM ÞLC  N M Þ r 0

ð3Þ

This is used as the criterion for enabling the transition from the classical continuum computational homogenization to an enhanced multi-scale scheme enriched by the discontinuity. Simultaneously, condition (3) also provides the (effective, averaged) !d orientation N M of the microscale strain localization band to be adopted as the orientation of the macroscale discontinuity.

Fig. 1. Classical computational homogenization framework, valid prior to strain localization.

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1489

Fig. 2. Continuum macrostructure with displacement discontinuity (center) and heterogeneous microstructure associated with a macroscopic material bulk point (left) and microscopic material volume surrounding a point on the macroscopic discontinuity (right).

The enrichment of the macroscale by a discontinuity introduces a new degree of freedom in the deformation field, i.e. a displacement jump or crack opening. In the remainder of this section the governing equations of this discontinuity enriched macroscale boundary value problem will be derived. The required constitutive relations will be extracted from the mechanical response of the underlying microstructure by extending the multi-scale scheme and corresponding scale transition relations. 2.1. Macroscopic boundary value problem !d Consider a macroscopic body O0M with a discontinuity Gd0M with local normal N M , Fig. 2 (middle). The deformation gradient tensor FM in the continuum (bulk) part O0M \Gd0M of the body is conventionally defined as the gradient with respect to the reference configuration of the current position vector field, i.e. ! FM ¼ ðr0M x M Þc in O0M \Gd0M ð4Þ The displacement jump at the discontinuity is defined as the difference between the current position vector from the two sides (‘‘ þ’’ and ‘‘ ’’) of the discontinuity ! þ ! ! ½½ u M  ¼ ð x M  x M Þ

on Gd0M

ð5Þ

The stress state in the body will be described by the first Piola–Kirchhoff stress tensor PM . To satisfy the stress balance on the internal boundary created by the displacement discontinuity, the tractions across the discontinuity have to be continuous and in balance with the ‘‘bulk’’ stress, i.e. d !d þ ! PM  N M ¼ P M  NM

!d !d þ !d with N M ¼ N M ¼  N M

ð6Þ

In the multi-scale framework the discontinuity displacement jump will be extracted from the microstructure, as will be ! explained in Section 3.4. The effective microstructural displacement jump is denoted by ½½ v M . For the converged solution ! of the macroscale boundary value problem the macroscale displacement jump ½½ u M  and the effective microscale ! displacement jump ½½ v M  should be equal. In the solution procedure, this equality will be enforced by a Lagrange ! multiplier l M . With the above definitions, the weak form of the macroscopic boundary value problem can be written as Z Z Z ! ! ! ! ! PM : dFcM dO0M þ dð l M  ð½½ u M ½½ v M ÞÞ dGd0M ¼ t 0M  d x M d@O0M ð7Þ O0M \Gd0M

Gd0M

@O0M

where the left hand side represents the internal virtual work. The right hand side of (7) expresses the external virtual work, ! with t 0M tractions (per unit reference undeformed area) acting on the external boundary @O0M of the body. The body forces have been omitted for brevity. Expanding the brackets in the second term of (7) leads to Z Z Z Z ! ! ! ! ! ! ! PM : dFcM dO0M þ d l M  ð½½ u M ½½ v M Þ dGd0M þ l M  d½½ u M  dGd0M ¼ t 0M  d x M d@O0M ð8Þ O0M \Gd0M

Gd0M

! 8d x M 2 U,

! ! U ¼ fd x M 2 C0 ðO0M \Gd0M Þ : d x M ¼ 0

! 8d½½ u M  2 V, ! 8d l M 2 L,

Gd0M

@O0M

on @Ou0M g

! V ¼ fd½½ u M  2 C1 ðGd0M Þg

! L ¼ fd l M 2 C1 ðGd0M Þg

where @Ou0M denotes the part of the external boundary on which displacement boundary conditions are applied. In (8) it ! has been taken into account that ½½ v M  is a constitutive quantity that will be obtained through the micro-to-macro scale transition (together with the stress PM ) and thus there is no corresponding test function.

1490

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

To obtain the corresponding strong form of the macroscopic boundary value problem, the chain rule and the divergence theorem can be applied to the first term in (8) resulting in Z Z !O0M !d ! ! N M  PcM  d x M d@O0M  N M  PcM  d½½ u M  dGd0M Gd0M

@O0M

Z  O0M \Gd0M

Z þ

! ðr0M  PcM Þ  d x M dO0M þ

! Gd0M

!

l M  d½½ u M  dGd0M ¼

Z

Z

!

Gd0M

!

!

d l M  ð½½ v M ½½ u M Þ dGd0M

! ! t 0M  d x M d@O0M

ð9Þ

@O0M

The standard argument that the above expression has to hold for all admissible test functions leads to the strong form of the macroscopic boundary value problem for a body with a displacement discontinuity ! r0M  PcM ¼ 0 local equilibrium in O0M \Gd0M ð10aÞ !O0 ! N M  PcM ¼ t 0M

tractions on the external boundary @O0M

ð10bÞ

! ! ! ½½ u M ½½ v M  ¼ 0 displacement  jump relation on Gd0M

ð10cÞ

!d

!

!

traction equilibrium on discontinuity Gd0M ð10dÞ ! From (10d) it is clear that the Lagrange multiplier l M is actually the cohesive traction which is in balance with the bulk stress PM .

l M  N M  PcM ¼ 0

2.2. Macroscopic constitutive relations To complete the above macroscopic boundary value problem, constitutive closure is required for the stress PM and the ! effective microscale displacement jump ½½ v M  as functions of the macroscale ‘bulk’ deformation field FM and the ! displacement jump ½½ u M  at the discontinuity. The macroscopic constitutive description will be obtained numerically by upscaling the response of the underlying microstructure through dedicated scale transition relations (to be derived in the next section). Here, the general form of the required constitutive relations will be established. For a macroscopic material point in the continuous part of volume O0M \Gd0M , the stress–strain response is derived according to classical computational homogenization approach, as schematically shown in Fig. 1. In this approach a macroscopic point is coupled to a nonlocalizing microstructural volume V 0m , see Fig. 2 (left). The mechanical response of the discontinuity is described by the stress PM at the discontinuity and the effective ! microscopic displacement jump ½½ v M . Both variables are obtained by modeling a single finite microstructural volume V 0m . This microstructural volume is crossed by a strain localization band, which is the microscopic physical realization of the macroscopic displacement discontinuity on Gd0M . In Fig. 2 (right), such a coupling of a macroscopic volume crossed by a discontinuity to an underlying microstructure crossed by a strain localization band is schematically depicted. In this work it has been assumed that a single localization band prevails within the microstructural volume, which is justified for random microstructures or periodic microstructures with an imperfection. Mathematically this means, that for each macroscopic point, condition (3) may only have one solution. Furthermore, it is assumed that within the microstructural volume the strain localization band can be centered on average on Gd0m , a straight line in 2D and a plane in 3D. The normal !d to this line (or plane), N M , corresponds to the local normal of the macroscale discontinuity Gd0M . e M of a small material volume V 0m containing a discontinuity Gd From a macroscopic point of view, the deformation F 0m is given by !d e M ¼ FM þ d0M ½½! F u M   N M

ð11Þ

where FM is, as before, the continuous bulk deformation and the second term accounts for the contribution of the !d ! discontinuity; d0M is the Delta-Dirac function centered at the discontinuity. Assuming that FM , ½½ u M  and N M are constant e M is computed as within the volume V 0m , the volume average of deformation F 1 V 0m

Z

d

!d e M dV 0m ¼ FM þ JG0m J ½½! F u M   N M V 0m V 0m

! where the following property of the Delta-Dirac function dGd , valid for a continuous function f ð x Þ, has been used Z Z ! ! ! dGd ð x Þf ð x Þ dV ¼ f ð x Gd Þ dGd V

Gd

ð12Þ

ð13Þ

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1491

Fig. 3. Multi-scale homogenization scheme, incorporating localization and damage.

! The geometric ratio appearing in (12) in front of the discontinuity ½½ u M  crossing the microstructural volume is the ratio between the measure of the support of the effective strain localization band JGd0m J (a length in 2D and area in 3D) and the microstructural volume V 0m . This ration will be denoted

bM ¼ JGd0m J=V 0m

ð14Þ

The whole volume averaged deformation will be called the macro-homogeneous deformation F bM !d ! F bM ¼ FM þ bM ½½ u M   N M

ð15Þ

where the superscript b indicates that this deformation is related to a certain discontinuity size to volume ratio bM . The constitutive relations for a macroscopic infinitesimal material volume containing the discontinuity Gd0M are written in a general form as ! ! ð16Þ PM ¼ GPMM fF bM ð X , tÞ, t 2 ½0,t, X 2 Gd0M g ! ! !½½ v M  b ! ! fF M ð X , tÞ, t 2 ½0,t, X 2 Gd0M g ½½ v M  ¼ G M

ð17Þ

while for the material points away from the discontinuity the constitutive relation is obtained through a classical computational homogenization scheme and has a general form. ! According to (10d) the continuum stresses PM and the traction at the discontinuity l M will automatically be in equilibrium for a converged state of the microstructural volume. Therefore, no independent constitutive relation is required for the tractions at the discontinuity. Hence, equilibrium relation (10d) is a priori satisfied. To solve the macroscopic boundary value problem iteratively, the following tangent operators are required: dPM ¼ 4 CM : ðdF bM Þc

ð18Þ

! d½½ v M  ¼ 3 DM : ðdF bM Þc

ð19Þ

The enriched multi-scale scheme is schematically shown in Fig. 3. Upon the moment of enrichment, indicated by the fulfillment of condition (3), the macro–micro coupling switches from a continuum deformation FM into the macrohomogeneous deformation F bM . The effective deformation is imposed onto the underlying microstructure through kinematical constraints, which will be discussed in more detail in Section 3.1. Hereby, the microscale boundary value problem remains a standard quasi-static equilibrium formulation. The micro–macro scale transition, as before, transfers ! the homogenized stress PM . In addition, the microscopic displacement jump ½½ v M  is derived from the microscale deformation and returned to the macrolevel. In order to solve the macroscale boundary value problem iteratively, the consistent tangents 4 CM and 3 DM are provided as well. Note, that in the case when there is no microscale localization band and, consequently, no macroscopic discontinuity, the ratio bM equals zero, i.e. the macro-homogeneous deformation F bM is equal to the bulk deformation FM . In other words, classical computational homogenization is a special case of the homogenization strategy proposed here. Obviously, the upscaling of the ! microscale effective displacement jump ½½ v M  is not applicable in the case of classical computational homogenization. 3. Microstructural modeling In this section, the microscale boundary value problem and the required scale transitions for a microstructural volume element crossed by a strain localization band are derived.

1492

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

The finite microstructural volume V 0m is in a state of equilibrium. In terms of the microscale first Piola–Kirchhoff stress tensor Pm the local momentum equilibrium in the absence of body forces is written as !

r0m  Pcm ¼ 0

ð20Þ

here r0m is the gradient operator with respect to the initial configuration of the microstructural volume. The material behavior of each microstructural constituent a (e.g. matrix, inclusion, etc.) is described by its proper constitutive law, specifying a time and history dependent stress–deformation relation for each microstructural constituent ðaÞ Pm ðaÞ ðaÞ ¼ Gm fFm ðtÞ, t 2 ½0, tg Pm

ð21Þ

It is emphasized that the present framework is not restricted to any particular constitutive law; the microstructural material behavior may be arbitrary and include a physical and/or geometrical evolution of the microstructure (e.g. phase transitions, decohesion, delamination, inter or intragranular fracture, etc.). The MVE problem is completed by the essential boundary conditions for the external boundaries of the microstructural volume. The boundary conditions are kinematical constraints that impose the macroscopic deformation onto the microstructural volume. After definition of the microscale kinematical quantities in Section 3.1, the macro-to-micro scale relation will be derived in Section 3.2. Upon solution of the microstructural boundary value problem, micro-to-macro scale ! transitions are required to upscale the stress PM and the effective displacement jump ½½ v M  based on the microstructural response. These scale transitions are detailed in Sections 3.3 and 3.4, respectively. 3.1. Microscale kinematics In Fig. 4, a microstructural volume V 0m and its corresponding macroscale volume are schematically shown. At the macroscale, the discontinuity divides the volume into two regions. The displacement field along the line perpendicular to !d the discontinuity (along the direction N M ) can be modeled by a Heaviside function H0M , which is a discontinuous step function. Here, this step function is defined to be zero on one side of the discontinuity and one on the other side. This is indicated in Fig. 4 by the white and black colors, respectively. The spatial derivative of the Heaviside function in the !d N M -direction is the Delta-Dirac function d0M . In Fig. 4, the profiles of H0M and d0M are plotted for the cross section a0 –a00 . Note that these functions are defined with respect to the coordinates of the undeformed configuration. At the microscale, the displacement field does not include a sharp discontinuity but a strain localization band. This implies that at the microscale the sharp macroscopic displacement jump is spread out over a certain width, i.e. a band with high strains. The displacement field resulting from such a smeared displacement jump can be modeled by a regularized or !d smooth step function H0m , see Fig. 4. This function is differentiable and therefore its gradient r0m H0m ¼ d0m N M is defined in the full domain V 0m , where d0m is a smooth equivalent of the Delta-Dirac function, which has the following property: Z 1 JGd J d0m dV 0m ¼ 0m ¼ bM ð22Þ V 0m V 0m V 0m where the definition (14) for bM has been used. For further derivations it is convenient to introduce another continuous scalar valued function DB0m defined on the microstructural volume according to !d

!

DB0m ¼ DH0m bM N M  D X m

Fig. 4. Enrichment of the displacement field in the macroscopic (top) and microscopic (bottom) domains.

ð23Þ

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1493

In Fig. 4 the profile of the function B0m in the cross section b0 -b00 is shown. The gradient of this function is given by !d

r0m DB0m ¼ ðd0m bM Þ N M

ð24Þ

where the property (22) has been used. The functions D H0m and DB0m depend on the position and the orientation of the strain localization band within an MVE. To estimate this position and orientation in the actual MVE simulations, an image analysis technique, called Hough transform, is used, see Appendix A. ! ! The material line elements in the initial D X m and the current D x m configurations of the microstructural volume characterized by ratio bM are described by ! ! ! D x m ¼ F bM  D X m þ D w m ð25Þ where the finite material vectors are defined with respect to a reference point r within the microstructural volume ! ! !r ! ! !r ! ! !r element, i.e. D X m ¼ X m  X m , D x m ¼ x m  x m , and D w m ¼ w m  w m . The first term in the right hand side of (25) describes the homogeneous deformation expressed in terms of the macro-homogeneous deformation tensor F bM , i.e. the ! homogeneous deformation of a microstructural volume with a localization band ratio bM . The vector field D w m in the second term represents the deviation from that uniform deformation inside the volume V 0m . It contains two contributions, both resulting from the heterogeneous microstructure. One part reflects the ‘collective’ microstructural strain response leading to the strain localization band development, the other part includes all remaining microfluctuations, denoted by ! D z m . The collective microstructural strain response will be described using the function DB0m , scaled by the value of the ! ! effective microscale displacement jump ½½ v M . The vector field D w m can thus be expressed as the superposition of the displacements resulting from the strain localization and the microfluctuation field, i.e. ! ! ! D w m ¼ DB0m ½½ v M  þ D z m ð26Þ In Fig. 5, a one-dimensional interpretation of Eqs. (25) and (26) is given. Fig. 5(a) shows the displacement field ! ! ! !

D u m ¼ D x m D X m as a function of the initial position D X m . The macro-homogeneous deformation F bM is imposed in a

weak sense on the microstructural volume V 0m . This implies that only the relative displacements of the left and right materials point are imposed. The dotted black line indicates a homogeneous deformation of the microstructural volume. The bold black line is the resulting fine scale displacement field. The difference between the dotted line and the bold black line is the result of the heterogeneous deformation in the microstructure. This field is plotted separately in Fig. 5(b) and can be split into the contribution from the strain localization and an additional microfluctuation field. The strain ! localization part is shaded with light gray and follows the shape of the function DB0m scaled with the magnitude of ½½ v M . ! ! The remaining microfluctuation field D z m is shaded in dark gray. The field D w m is the combination of the two gray shaded areas, i.e. the difference between the real microstructural deformation and the homogeneous deformation of the microstructural volume. 3.2. Macro-to-micro scale transition: microscale boundary conditions Based on the microscale position vector field in the current configuration, as defined in Eq. (25), the microscopic deformation gradient Fm is calculated as ! ! Fm ¼ ðr0m D x m Þc ¼ F bM þ ðr0m D w m Þc ð27Þ

Fig. 5. Different contributions to the microscale displacement field for a one dimensional case with strain localization in Gd0m .

1494

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

The scale transition relation postulated here requires the volume average of the microstructural deformation tensor Fm to be equal to the macro-homogeneous deformation F bM , i.e. Z 1 Fm dV 0m ¼ F bM ð28Þ V 0m V 0m Substitution of Eq. (27) in relation (28) with application of Gauss’ theorem results in the following constraint on the ! microscale D w m : Z Z ! ! ! ðr0m D w m Þc dV 0m ¼ D w m  N m d@V 0m ¼ 0 ð29Þ V 0m

@V 0m

! here @V 0m is the external boundary of the microstructural volume V 0m with outward normal N m . For the continuous function DB0m the following result holds: Z Z Z ! ! !d DB0m N m d@V 0m ¼ r0m DB0m dV 0m ¼ ðd0m bM Þ dV 0m N M ¼ 0 V 0m

@V 0m

ð30Þ

V 0m

where first Gauss’ theorem has been applied, followed by substitution of (24) and (22). With (30), substitution of the ! definition (26) into constraint (29) allows to reformulate the constraint on the field D w m into an equivalent constraint on ! the microfluctuation field D z m Z Z ! ! ! ðr0m D z m Þc dV 0m ¼ D z m  N m d@V 0m ¼ 0 ð31Þ V 0m

@V 0m

! Thus the constraints (29) and (31) are equivalent. Since D w m is the primary unknown in the MVE problem, while the ! ! displacement jump ½½ v M  and the microfluctuation field D z m are only separated in a post-processing step at the end of ! ! each increment, the constraint (29) on D w m is preferred in the implementation over the constraint (31) on D z m . ! ! Note that prior to localization, the first term in (26) vanishes and the microfluctuations D z m and D w m are identical. The constraint (29) is then identical to the constraint used in classical homogenization on representative volume elements. The constraint (29) can be imposed on an MVE in many alternative ways. This constraint should be applied in such a way that it minimally interferes with the emerging strain localization and, at the same time, it needs to provide a good estimate for the effective stress–strain response. The boundary conditions that are most often used in MVE analyses are summarized in Table 1. In the listed order, the effective stiffness obtained by applying the different boundary conditions will, in general, decrease from over- to underestimated. This is due to the decreased amount of kinematical constraints imposed onto the microstructural volume. On the other hand, with decreasing number of constraints, i.e. with increasing kinematical freedom, the freedom for a strain localization band to develop increases. However, the sensitivity to spurious localization modes also increases. This influences both the mode and the initiation (where and when) of the predicted macroscale localization, as defined in Section 2. Recently, Coenen et al. (2012, 2011) proposed so-called percolation path aligned boundary conditions. Like periodic boundary conditions, these constraints kinematically tie material points located at opposing boundaries, which provides a good effective stress–strain estimate. In the new boundary conditions, the opposing boundaries do not have to coincide with the opposite sides of a (rectangular) microstructural volume element, because this would only allow for a horizontal or vertical strain localization bands (and some combinations of periodic bands). Therefore, in the new approach the opposing boundaries are defined according to the direction of the (developing) strain localization band (strain percolation path). The detailed description of these boundary conditions is beyond the scope of this paper and can be found in Coenen et al. (2012, 2011). The multi-scale framework presented in this paper is as such independent of the particular choice for the boundary ! conditions. It is only required that kinematical constraints in terms of D w m are imposed on the MVE boundary such that constraint (29) is satisfied. 3.3. Micro-to-macro scale transition: macroscopic stress PM After solving the microstructural boundary value problem, the microscale response needs to be upscaled to the macrolevel. In order to define the macroscopic stress, the Hill–Mandel condition is typically used (Hill, 1963; Suquet, Table 1 Classical boundary conditions for microstructural volume elements. BC type

! Constraint on D w m

Taylor Linear Periodic Minimal (traction)

Effective stiffness

Localization band

D w m ¼ 0 8 X m 2 V 0m

Strong overestimation

Fully suppressed

! ! ! D w m ¼ 0 8 X m 2 @V 0m !þ ! !þ ! ! D w m ¼ D w m 8 X m 2 @V 0m s.t. N m ¼  N m R ! ! @V 0m D w m  N m d@V 0m ¼ 0

Overestimation

Constrained at the boundaries

Good estimation

Respecting periodicity

Underestimation

Sensitive to spurious localization

!

! !

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1495

1985). This condition requires the variation of the microstructural work per unit volume to be equal to the local variation of work at the macroscale. Here, the Hill–Mandel condition is modified, such that it requires the variation of macroscopic work for a volume with localization ratio bM to be equal to the microstructural virtual work volume averaged over a microstructural volume element that is consistent with this ratio bM , i.e. Z 1 Pm : dFcm dV 0m ¼ PM : ðdF bM Þc 8dF bM ð32Þ V 0m V 0m Applying the divergence theorem and taking into account the equilibrium of the microstructure (20) results in the following expression for the variation of volume specific microstructural work: Z Z 1 1 ! ! Pm : dFcm dV 0m ¼ p  dD x m d@V 0m ð33Þ V 0m V 0m V 0m @V 0m m ! ! where p m ¼ N m  Pcm are the tractions on the boundary of the microstructural volume. Substitution of the variation of Eq. (25) results in   Z Z Z ! 1 1 1 ! ! ! ! ! p m  dD x m d@V 0m ¼ p m  D X m d@V 0m : ðdF bM Þc þ p  dD w m d@V 0m ð34Þ V 0m @V 0m V 0m @V 0m V 0m @V 0m m Since the boundary tractions resulting from the imposed constraint (29) do not contribute to the external work, the last term in the above expression vanishes Z ! ! p m  dD w m d@V 0m ¼ 0 ð35Þ @V 0m

Comparing relations (32) and (34) and taking into account (35), the macroscopic stress measure PM is obtained as Z Z ! 1 1 ! PM ¼ p m  D X m d@V 0m ¼ Pm dV 0m ð36Þ V 0m @V 0m V 0m V 0m In the last term the divergence theorem has been used to transform the boundary integral expression for the macroscopic stress to a volume integral. From (36) it is clear that the macroscopic stress equals the volume average of the stress in the microstructure (like for the standard computational homogenization scheme). ! 3.4. Micro-to-macro scale transition: displacement jump ½½ v M  For the macroscopic boundary value problem, not only the homogenized stress response needs to be returned to the ! macrolevel, but also the displacement jump ½½ v M . This requires a scale transition that discriminates ‘bulk’ from ‘localization’ deformation within the microstructural volume, i.e. between the dark and light gray areas in Fig. 5. Prior to detection of macroscale localization, the deformation of the microstructure consists of bulk macroscopic deformation only, with a superimposed microfluctuation field. At the onset of macroscale localization, a displacement discontinuity is ! ! introduced at the macroscale, with an initial zero displacement jump, i.e. ½½ u M  ¼ 0 . To preserve consistency with the ! macroscale, the corresponding microscale displacement jump ½½ v M  should be such that upon the moment of macroscale !n localization it yields a zero vector. To this purpose, we first introduce the relative displacement field D u m relating the !L current position vectors to those at the onset of macroscale localization D x m (the superscript L is used to mark this onset). !n ! !L Based on the position vector field defined in Eqs. (25) and (26) this displacement field D u m ¼ D x m D x m is written as ! !n !n ! D u m ¼ F bMn  D X m þ DB0m ½½ v M  þ D z m ð37Þ !n Here the superscript n indicates that a quantity is taken relative to its value at the onset of localization, i.e. D z m ¼ ! ! !L ! D z m D z m , F bMn ¼ F bM F bML , and ½½ v M L ¼ 0 . !n After the microstructural boundary value problem is solved, the displacement field D u m is known. The scalar-valued !n function DB0m is defined in Eq. (23). Then, the realized microfluctuation field D z m can be expressed in terms of the known ! ! fields and the yet unknown displacement jump ½½ v M . The realized displacement jump ½½ v M  is then resolved through a !n least square minimization of the realized microfluctuation field D z m . In the standard least squares terminology, the ! !n !n ! microfluctuation field D z m is the residual, D u m F bMn  D X m is the data set and ½½ v M  is the desired parameter. This leads to the following minimization problem: Z Z !n ! @D z m !n !n !n min D z m  D z m dV 0m ) ð38Þ !  D z m dV 0m ¼ 0 ! V 0m V 0m @½½ v M  ½½ v M  The derivative of the microfluctuation field with respect to the displacement jump following from (37) equals !n @D z m ¼ DB0m I ! @½½ v M 

ð39Þ

1496

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

!n Substitution of the microfluctuation D z m defined by Eq. (37) and its derivative (39) into Eq. (38) results in Z ! ! !n ! DB0m ðD u m F bMn  D X m DB0m ½½ v M Þ dV 0m ¼ 0

ð40Þ

V 0m

! ! Since (40) is linear in ½½ v M , it allows for a closed-form solution given by the following expression for ½½ v M : Z ! 1 !n ! DB0m ðD u m F bMn  D X m Þ dV 0m ½½ v M  ¼ S V 0m here the scalar S is defined as Z S¼ DB20m dV 0m

ð41Þ

ð42Þ

V 0m

This completes the formulation of the continuous–discontinuous computational homogenization–localization framework. Some essential details on the implementation can be found in Appendix B. 4. Numerical examples The proposed computational homogenization–localization scheme is illustrated in this section on two simple benchmark examples that demonstrate the key features of the multi-scale framework. The first considered example is quasi two-dimensional, allowing for a comprehensive interpretation and evaluation of the results, as well as the comparison with the fully resolved direct numerical simulation. The second example illustrates the applicability of the framework to fully two-dimensional problems. 4.1. Quasi two dimensional example 4.1.1. Problem description ! The problem consists of a slab of voided elasto-plastic material of length LM ¼ 150 mm in e 1 direction and infinitely ! ! ! ! large in e 2 and e 3 -direction s, see Fig. 6(a). In e 2 and e 3 directions, a plane strain state is assumed. The left edge of the ! sample is fixed. The displacement DLM of the right edge in e 1 -direction is prescribed. The macroscale domain is discretized by Nelem one dimensional elements with linear displacements and constant approximations for the displacement jump and the traction Lagrange multiplier. One Gauss integration point per element is used, indicated by ‘‘  ’’ in Fig. 6(a). With each macroscale Gauss point an MVE is associated, which has dimensions 1  0.5 mm2, see Fig. 6(b). Note that the MVE can be considered as an RVE (i.e. representative) prior to localization only. The MVE contains 41 voids with a radius of 0.03 mm. The matrix material is an isotropic von Mises elasto-plastic material with linear hardening. The elastic stiffness and Poisson’s ratio are E¼200 GPa and n ¼ 0:3, respectively. The plastic properties are the initial yield stress soy ¼ 500 MPa and hardening modulus H¼60 MPa. All MVE’s are identical, except for one corresponding to the macroscopic point indicated by a ‘‘%’’ in Fig. 6(b). The effective thickness of this MVE in the ! e 3 -direction is decreased by 1% to create an imperfection that triggers localization. In the macroscopic 1D continuum, only a single direction of the discontinuity is realizable that is known a priori, i.e. !d ! ! along the e 2 -direction (with normal N M ¼ e 1 ) in this case. For this reason, the standard periodic boundary conditions on the MVE have been used for this problem. This, combined with the imposed uniaxial deformation loading condition, ensured that only vertical localization bands exist within the MVE. ! In the following, first the extraction of the effective microstructural discontinuity jump ½½ v M  will be illustrated on an MVE analysis. After that the results of the coupled multi-scale simulation of the above problem will be presented. 4.1.2. Extraction of the effective displacement jump ! To illustrate the procedure to extract the displacement discontinuity ½½ v M  from the microscale MVE response, a uniaxial deformation state consistent with the macroscopic loading has been applied on the MVE shown in Fig. 6(b). The

Fig. 6. (a) Schematic representation of the macroscale problem and loading conditions and (b) underlying microstructural volume element.

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1497

Fig. 7. (a) Homogenized stress–deformation curve of the MVE (the cross indicates moment of macroscale localization). (b) Deformed MVE geometry at an ! b overall deformation of F 11 M ¼ 1:08. Indicated by the gray scale are the nodal displacements (in the e 1 -direction) relative to their values at the onset of macroscale localization.

prescribed macroscale deformation tensor is given by ! ! ! ! ! b! F bM ¼ F 11 M e 1  e 1þ e 2  e 2þ e 3  e 3

ð43Þ

b F 11 M

where the component contains both the bulk deformation and the macroscale displacement jump. The parameter bM for this problem is calculated as

bM ¼

JGd0M J H ¼ 1 mm1 ¼ V 0m WH

ð44Þ

Fig. 7(a) shows the homogenized stress response as a function of the prescribed deformation. Indicated with a marker is b 11 the onset of localization at an overall deformation of F 11 M ¼ FM ¼ 1:071. The moment of localization is determined by 11

11

evaluating Eq. (3), which in a 1D case reduces to dPM =dFM o0. Fig. 7(b) shows the deformed MVE geometry at an overall b F 11 M

deformation of ¼ 1:08, i.e. after the onset of strain softening. The strain localization band within the MVE is clearly visible. At the onset of localization, a discontinuity will be introduced at the macroscale and the effective displacement jump ! ½½ v M  calculated from the microscale displacement field is required. The displacement jump is calculated according to relation (41) which makes use of the scalar field DB0m . This field, defined in (23), can be expressed as a function of the local !d coordinate rm along the normal N M with respect to the geometrical center c of the MVE, defined by !d

!

!c

rm ¼ N M  ð X m  X m Þ In order to calculate the scalar field DB0m the following smooth step function H0m has been used:   1 1 r rLm H0m ðrm Þ ¼ þ tanh m 2 2 lm

ð45Þ

ð46Þ

with rLm the position of the strain localization band within the MVE, which is determined by the Hough transform, as described in Appendix A. The difference rm rLm is in fact the signed distance function of an MVE material point with respect to the localization band. The width of the microscale localization band is described by the parameter lm ; a procedure to identify this parameter will be discussed in the following. ! In Fig. 8(a) the nodal displacements in e 1 -direction of all nodes in the MVE are plotted (as dots) against their position b rm at an overall MVE strain deformation of F 11 M ¼ 1:08. Like in Fig. 7(b), the nodal displacements shown are relative to the corresponding values at the onset of macroscale localization, i.e. u1mn . Based on this displacement field, the effective displacement jump ½½v1M  can be computed using relation (41). In order to identify the localization band width parameter lm , the jump values determined for different values of lm are shown in Fig. 8(b). It can be seen that for the values of lm between very small (that would correspond to a sharp Heaviside function (46)), to approximately 0.1 mm the computed jump is almost constant. Based on this, the value lm ¼ 0:05 mm, here resulting in ½½v1M  ¼ 0:0086 mm, has been selected to be used in the rest of the simulations. Note, that the identified value of lm is in the range of the void size in this example (0.06 mm diameter). The solid line in Fig. 8(a) represents the coarse scale displacement field given by (37). In Fig. 9 the numerically obtained relation between the homogenized stress component P11 M and the effective displacement jump ½½v1M  is plotted. The stress decreases with increasing displacement jump, which implies that the load carrying capacity of the microstructure is decaying. The stress-discontinuity opening relation in Fig. 9 can be interpreted as a homogenized cohesive relation.

1498

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

Fig. 8. (a) MVE nodal displacements (dots) relative to their values at the onset of macroscale localization. The solid line represents the resulting ! ‘homogenized’ displacement field of the MVE. (b) The calculated jump ½½ v M  for different values of the width parameter lm . The dot corresponds to the value lm ¼ 0:05 mm, used in the simulations.

Fig. 9. Homogenized stress–displacement jump response.

4.1.3. Multi-scale simulations In this section, the multi-scale analysis of the problem schematically shown in Fig. 6 will be presented. The results obtained with the classical first-order homogenization and the newly proposed discontinuity enriched scheme are confronted and verified against a full scale direct numerical simulation (DNS) solution (which is feasible for the simple academic problem considered here). The full-scale DNS solution has been obtained by directly discretizing the heterogeneous macroscale problem ! consisting of 150 MVEs stacked next to each other in e 1 direction. One MVE has an out-of-plane thickness decreased ! by 1%. In the e 2 -direction periodic boundary conditions have been used to mimic the ‘infinite’ dimension. Displacements of the left edge are constrained and an indirect displacement control (de Borst, 1987) has been used to follow the snapback behavior and obtain a fully converged solution after the point of localization. Macroscale results obtained with a classical (first-order) computational homogenization scheme are presented in Fig. 10(a). The tractions tM as a function of the elongation DLM =LM are shown for different discretizations of the macroscale problem (with Nelem the number of elements used). It is clear that prior to localization the results are mesh independent. After localization the classical scheme predicts an increasingly brittle macroscopic response upon mesh refinement. At first glance, upon mesh refinement, the classical computational homogenization solution seems to converge to the reference solution (dashed line in Fig. 10(a)). This is due to the fact that in the particular case when the macroscopic element length equals the MVE width, the result of the classical scheme is indeed close to reference response. In this case however, the advantage in computational effort of computational homogenization is completely lost since the separation of scales is not present anymore. The situation then resembles a multigrid or a coupled volumes technique. Moreover, further macroscopic mesh refinement results in an overestimation of the snap-back compared to the reference solution. Fig. 10(b) shows the macroscopic traction–elongation response obtained with the enriched multi-scale scheme. Prior to localization the obtained response is exactly on top of the reference solution and identical to the one obtained with the classical homogenization scheme. After localization, in contrast to the classical scheme, the results obtained with the enriched scheme are in close agreement with the reference solution. Upon mesh refinement, the slope of the unloading branch after localization is unaltered and all discretizations predict the same traction for the onset of localization. In terms of elongation, a small difference is visible, which is only due the fact that the imperfection was applied to a single element, the size of which varies for each discretization.

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1499

Fig. 10. Macroscale traction–elongation response for different macroscale mesh sizes obtained by (a) classical and (b) discontinuity enriched computational homogenization; the reference DNS solution is also shown.

Fig. 11. Microscale response obtained by the multi-scale modeling compared to the DNS reference full scale solution.

In addition to the macroscale response, computational homogenization also provides detailed information on the microscale, e.g. the local stress–strain response, microstructural deformation, and distributions of microscale fields. In Fig. 11, the macroscale mesh with N elem ¼ 3 for the multi-scale analysis is shown with the three MVEs corresponding to the Gauss points. This is done for two instances: prior to localization (denoted by J) and after localization (). The equivalent von Mises contour plots are shown on the deformed MVE geometries. In addition, the stress versus deformation MVE responses of the Gauss points are shown. The central MVE contains the strain localization band which results in a softening response of this MVE and leads to the unloading of the other MVEs. The shaded gray area in the stress–strain curve of the softening MVE is the deformation attributed to the discontinuity opening ½½v1M . The left edge of this shaded gray area reflects the unloading curve of the bulk material within the MVE. (Note, however, that this MVE does not unload but follows the solid line of increasing deformation and decreasing stress.)

1500

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

At the bottom of Fig. 11 details of the equivalent von Mises stress distribution of the DNS full-scale reference solution are shown. These volumes are cropped from the long heterogeneous strip modeled for the reference solution and correspond to the respective MVEs. As can be noticed, the microstructural responses obtained by the full-scale and multiscale solution procedures show a very good agreement, both prior (J) and after () the localization. 4.2. Two-dimensional example The second numerical example illustrates the applicability of the framework in a fully two-dimensional setting. Although a relatively simple example will be considered, it allows the illustration of the initiation and propagation of a non-trivial two-dimensional macroscopic crack, governed by the microstructural strain localization. 4.2.1. Problem description The considered macro-geometry is a rectangular area of dimensions 25 mm  20 mm, shown in Fig. 12(a), where also the applied boundary conditions are indicated. The bottom edge is fixed in the vertical direction. To the top edge a vertical displacement is imposed, which increases linearly from 12U y at the left corner to Uy at the right corner. The remaining nodes on the left and right edges are free vertically, while the horizontal displacements of these edges are constrained to be point-wise equal (through tying constraints). The macroscopic domain is discretized by 16 linear triangular plane strain elements, with one integration point per element, as shown in Fig. 12(a). After strain localization, the localized element is enriched with a displacement jump discontinuity. In the present implementation this enrichment is based on the embedded discontinuity concept (Simo et al., 1993; Armero and Garikipati, 1996; Wells and Sluys, 2000; Linder and Armero, 2007), extended to accommodate the micro–macro scale coupling proposed here. The details of this numerical implementation are beyond the scope of this paper and will be presented elsewhere. To each macroscale integration point an MVE is assigned, shown in Fig. 12(b). In the reference configuration the MVE dimensions are 1 mm  1 mm. It contains 80 cylindrical voids of 60 mm diameter. The matrix material is modeled as elasto-plastic according to an isotropic logarithmic hyperelasto-plastic law (Simo, 1992). The elastic material parameters are the shear modulus 80.19 GPa and the bulk modulus 164.21 GPa. The plastic flow is described by an isotropic von Mises model with exponential hardening. The involved parameters are the initial flow stress of 0.45 GPa, the residual flow stress 0.715 GPa, linear hardening coefficient 0.129 GPa, and saturation exponent 16.93. The MVE is discretized by 4484 bi-linear plane strain quadrilateral finite elements (average element size 0.013 mm), with selective reduced integration. The MVE FE mesh is not shown here for clarity. The macro-to-micro scale transition, i.e. the imposition of the macro-homogeneous deformation on the MVE, involves the percolation path aligned boundary conditions, as briefly discussed in Section 3.2 and presented in detail in Coenen et al. (2012). These boundary conditions allow for the development of a microscale localization band in a direction determined by the interaction between the microstructure and the macroscopic loading, rather than the initial periodicity directions (e.g. horizontal and vertical). Note that for this non-linear problem the direct numerical simulation is hardly feasible and would require massive (parallel) computing facilities. 4.2.2. Results of multi-scale simulations Fig. 13 shows in the center the deformed macroscale mesh with the enriched elements (white colored) at the end of the simulation (U y ¼ 0:55 mm); the thick black lines indicate the embedded cohesive crack discontinuities; the arrows represent the direction and the amount of the computed discontinuity openings. The crack initiates at the lower right element and propagates towards the upper left corner. Also shown in Fig. 13 are the deformed geometries of several MVEs

Fig. 12. (a) The geometry of the macroscale problem with the boundary conditions, the finite element mesh and integration points. (b) Undeformed geometry of a microstructural volume element.

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1501

Fig. 13. Center: the deformed (deformation  5) macroscale enriched elements (white) and the embedded discontinuities (thick black lines) at the end of the simulation (U y ¼ 0:55 mm). The black arrows indicate the displacement discontinuity jump; the length of the largest arrow is 0.14 mm, the smallest 0.05 mm. Insets show the deformed geometries (real displacements) of several MVEs corresponding to the indicated macroscopic points.

! ! Fig. 14. The undeformed macroscale mesh with the field of the macroscopic stress component in N direction. Arrows indicate the tractions l M at the discontinuity (shown by the thick black line); the length of the largest arrow is 234 MPa, the smallest 162 MPa.

corresponding to different macroscopic integration points. The strain localization bands within MVEs ‘a’–‘f’ are clearly visible, corresponding to the macroscopic discontinuity. The deformation of the MVEs associated to the macroscale elements not crossed by the crack is relatively small. ! ! The macroscopic normal stress orthogonal to the localization band PM : N  N is plotted on the undeformed ! macroscale mesh in Fig. 14. The direction N represents the average direction of the normals to the realized discontinuities and is shown in Fig. 14. The stress field reveals a pattern that is to be expected in the presence of a cohesive softening crack, i.e. as the crack propagates, the stress in the affected elements reduces. The arrows in Fig. 14 represent the ! computed Lagrange multipliers l M at the discontinuity. Indeed, in accordance with Eq. (10d), for the converged numerical ! solution, l M equals the tractions at the discontinuity. In Fig. 15, the solid black line shows the total macroscale reaction force on the top edge of the macroscopic domain as a function of the prescribed displacement Uy. The total force is a sum of the individual reaction forces, indicated by the gray shaded areas, at the top edge nodes 1–5. The dashed vertical lines in Fig. 15 mark the instances, at which the effective strain localization of an MVE has been detected (according to condition (3)) and the corresponding macroscopic element has been enriched with the embedded discontinuity. It can be seen that as the macroscopic cohesive crack propagates, and thus more macroscale elements are being enriched, the total force–displacement curve turns into a negative slope. Considering the individual nodal reaction forces, it can be seen that as the crack propagates from right to left, first the

1502

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

Fig. 15. Macroscale total reaction force at the top edge (black solid line) and the reaction forces at the individual nodes (gray shaded areas) versus prescribed displacement of the top edge. The dashed vertical lines indicate moments at which macroscale elements localized.

contribution of the right node (F 5y ) starts to decrease, followed by the decay of the reaction forces at the nodes more to the left. This pattern is indeed to be expected for the considered multi-scale example. 5. Conclusion This paper presented a multi-scale technique which kicks off from a classical computational homogenization scheme and gradually evolves towards a scheme that consistently upscales the effect of microscale strain localization towards a macroscale cohesive crack, thereby enabling the multi-scale modeling of localization and damage up to the point of macroscopic fracture. Like other computational homogenization schemes this techniques is especially useful for complex materials that undergo complex loading paths. The modeling departs from a classical computational homogenization scheme. The scheme relies on a microstructural volume element MVE, which was initially locally representative for the microstructure considered. The local macroscopic deformation gradient tensor is transferred to the microscale MVE using appropriate scale transition relations, leading to kinematical constraints on the boundary of the MVE. These extended constraints enable the combination of a reliable effective stress–strain response and the development of a strain localization band governed by the underlying microstructure. Upon solution of the MVE boundary value problem, the resulting stress tensor is recovered according to the Hill–Mandel condition. The macroscopic consistent tangent is extracted from the MVE tangent stiffness matrix. During the loading process a strain localization band is allowed to develop in the microstructure. Strain softening in the band degrades the load carrying capacity of the microstructure across the band, which is detected by an acoustic analysis of the homogenized tangent stiffness tensor. Strain localization limits the representative character of the MVE, thereby violating the assumed separation of scales, inevitably limiting classical homogenization approaches. To overcome this limitation, the macroscale continuum is enriched with a displacement discontinuity and consequently the scale transition in the multi-scale scheme is enhanced at the onset of localization. The kinematics of the macroscale domain is described in terms of the local deformation gradient tensor in the bulk of the material and a displacement discontinuity. From these two quantities a macro-homogeneous deformation tensor is composed that characterizes an MVE crossed by a discontinuity with a certain volume to discontinuity size ratio. Kinematical constraints are used to impose the effective deformation on the MVE. After the MVE analysis, the resultant stress and microscale displacement jump (lumping the microscale strain localization into a crack opening at the macroscale) as well as their respective consistent tangent operators are returned to the macrolevel. The local equilibrium remains applicable in the continuum part of the macroscale domain. At the discontinuity, the equality between the macroscale displacement and effective microscale displacement jump is enforced by Lagrange multipliers, equal to the cohesive tractions. Scale transition relations in macroscale material points and associated MVEs (RVEs) that are not localized remain classical. Two numerical examples have been presented illustrating the essential features of the proposed framework: a simple quasi-2D problem enabling a comparison with the reference DNS results and a fully 2D problem, featuring a non-trivial propagating cohesive crack. Although the framework assumes a single localization band crossing the MVE at the microscale (which applies to random or imperfect periodic microstructures), at the macroscale, multiple cracks originating from different underlying MVEs can run through the macroscale domain, capturing even more complex crack paths. This will be illustrated in the future work. The numerical examples shown in this paper are representative for ductile failure mechanism, i.e. the void

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1503

growth and plastic strain localization in the ligaments. The presented framework, however, also allows to include discontinuous failure events at the fine scale, e.g. micro-cracks. Finally, it is noted that the proposed multi-scale scheme can be combined with different finite element technologies for the discretization of the discontinuity enriched macroscale, e.g. embedded localization zones, extended finite element method (XFEM), etc.

Acknowledgment This research was carried out under project number MC2.05205b in the framework of the Research Program of the Materials innovation institute M2i (www.m2i.nl). Appendix A. Identification and quantification of strain localization using the Hough transform The Hough transform is a methodology to identify lines or other parameterized shapes in noisy images and provide an estimate of the associated parameters (Hough, 1959; da Fontoura Costa and Marcondes Cesar, 2009). The procedure makes ! ! use of a mapping of the image values (e.g. gray scale) from the image space ð e 1 , e 2 Þ into a parameter space. In this work the Hough transform is used to determine the location of the localization band within the MVE. The ‘image’ to be analyzed will consists of the values of some microscale field, e.g. an incremental equivalent strain, plotted on the MVE geometry. The feature to be identified is the band of the points with high strain values, i.e. the strain localization band. The identification of lines generally uses the (a, r) parameter space. Here r is the normal distance between the line ! and the coordinate origin and the angle a is defined between the line normal and the  e 2 -direction, see Fig. 16. All lines passing through a point ðX 1 ,X 2 Þ in the image space obey the following relation:

r ¼ X 1 cos a þ X 2 sin a

ðA:1Þ

This relation transforms each image point into a sinusoidal curve in the parameter space. In such a mapping, the sinusoidals of the points belonging to the same line in the image intersect at a single point in the parameter space. The coordinates ðaL , rL Þ of this intersection point provide the parameters of the line. In the implementation, both the image and the parameter space are not continuous. The parameter space is sampled with N a  Nr bins, which form the accumulator array Qij with Dr and Da the width of the bins in r and a directions, respectively. The center of a bin is given by ðab , rb Þ. The discrete image space consists of a finite number of points NX. In the accumulator Qij the image value ak is summed according to Q ij ¼

NX X

Ak ak

if ðrjb DrÞ r ðX 1k cos aib þ X 2k sin aib Þ r ðrjb þ DrÞ

ðA:2Þ

k¼1

where Ak is the weight of image point k. For an image consisting of a rectangular grid of points, generally, the weight of each point is chosen to be equal and Ak ¼1. However, if the image consists of an irregularly distributed set or image points, Ak can be chosen equal to the area represented by the image point k. If the image is such that it contains a line consisting of high intensity image points, these line points will create high value sinusoidals in the accumulator array. Where these sinusoidals cross, they will create a peak in the accumulator array Qij at the position of the parameters aib and rjb . These parameters represent the line with highest intensity in the image. In Fig. 17(a) the b incremental equivalent strain field for the overall deformation of F 11 M ¼ 1:08 of the MVE used in the example in Section 4.1 is shown. The corresponding Hough transform based on this strain field is displayed in Fig. 17(a). It shows a clear maximum at !d ! position rL ¼ 0:076 and aL ¼ p=2. These parameters represent a strain localization band with the normal N M ¼ e 1 and positioned to the left of the MVE center, which corresponds to the visually observed band in Fig. 17(a). In a cross section of the Hough transform surface at a ¼ p=2, shown in Fig. 17(c), the position of the localization band shows even more clearly.

Fig. 16. Hough transform mapping.

1504

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

Fig. 17. (a) Equivalent incremental strain plot on which the Hough transform is based; (b) the full (r, a) parameter space and (c) a cross section perpendicular to the strain localization band (aL ¼ p=2).

Appendix B. Numerical implementation The equilibrium equation (20), constitutive relations (21), and kinematical boundary conditions derived from the constraints (29) together form a standard quasi-static non-linear boundary value problem on the MVE. In this work, a conventional finite element procedure is used to solve this microstructural problem. In the following, some essential details of the numerical implementation of the microscale problem and scale transition relations are given. Here, a matrix and a column are denoted by A and a, respectively. Superscript T is used to indicate the transposed of matrices and columns. The individual component I of a ecolumn and IJ of a matrix are denoted by fagI and fAgIJ , respectively. B.1. Kinematical boundary conditions ! The current position vector field (25) is expressed in terms of a displacement field D u m . To eliminate the unknown ! !r !r displacement u m of the reference point, one node on the boundary is fixed, i.e. u m ¼ 0 , thereby suppressing rigid body translations of the MVE. The unknown displacement field is written as ! ! ! u m ¼ ðF bM IÞ  D X m þ D w m ðB:1Þ ! This expression is used to transform constraint (29) on D w m derived in Section 3.2, into the constraint on the ! displacement field u m Z Z ! ! ! ! u m  N m d@V 0m  D X m  ðF bM IÞc  N m d@V 0m ¼ 0 ðB:2Þ @V 0m

@V 0m

For example, for the periodic boundary conditions, as given in Table 1, this results in !þ ! ! þ ! u m ¼ u m þðD X m D X m Þ  ðF bM IÞc

ðB:3Þ b

For the numerical implementation of these constraints, the term containing the macroscale deformation F M I is rewritten by making use of control vectors. In case of a two dimensional analysis, two control vectors are necessary to contain all degrees of freedom of F bM . It is convenient to choose the control vectors in the undeformed reference configuration to coincide with the orthonormal coordinate base vectors, collected in a column AM AM 



! ! ¼ ½ e 1 , e 2 T

ðB:4Þ

The control vectors defining the current configuration are given by aM 

-

¼ ðF bM IÞ  AM 

ðB:5Þ

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1505

The macroscopic deformation F bM I can then be written as -T

-





F bM I ¼ aM  AM

ðB:6Þ

When implementing in a general purpose finite element code, these control vectors are easily introduced through control or ‘dummy’ nodes (i.e. not connected to the FE mesh). In some cases it is also possible to make the dummy nodes coincide with physical nodes, which allows to reformulate the periodicity equation (B.3) in a set of homogeneous tying relations. This is typically done in the classical computational homogenization scheme. The displacement field is discretized and represented by a column of nodal displacements um . Due to the finite number of  boundary nodes, the kinematical boundary conditions (B.2) result in a finite number of constraints. Each constraint eliminates the degrees of freedom of one dependent boundary node (denoted by d) from the system. The remaining boundary nodes are denoted by i (for independent), which also includes the reference node r. These degrees of freedom will be retained in the system. The discretization of constraint (B.2), with account for (B.6), results in a set of linear relations -d um 

-i

-d

-i

-T

-

-i

-















¼ D um þ ðD Xm D D Xm Þ  AM  aM ¼ D um þB aM

ðB:7Þ

Here D is the dependency matrix containing the coefficients of the linear constraints between the boundary nodes; matrix B contains the coefficients related to the control vectors -d

-i

-T







B ¼ ðD Xm D D Xm Þ  AM -i Xm 

ðB:8Þ

-d Xm 

Columns D and D contain the initial position vectors of the independent and dependent boundary nodes, respectively, relative to the reference point. The exact form of matrices D and B depends on the type of the boundary conditions used and, obviously, the geometry and discretization of the MVE. B.2. Macroscopic stress In order to calculate the macroscopic stress from the converged solution of the MVE problem, the boundary integral (36) is discretized resulting in PM ¼

1 -T g m  D Xm  V 0m 

ðB:9Þ

-

where g m are the forces acting on the boundary nodes combined into a column. This column can be partitioned into the  -a

constraint forces acting on the independent and dependent nodes and the forces f M working on the control nodes. Using 

the standard argument that no additional work should be added to the system by the tying constraints, the constraint forces can be eliminated from (B.9) reducing it to PM ¼

1 -aT f  AM  V 0m M

ðB:10Þ

Thus, only the variables related to the control nodes are involved in the final expression of PM , making it computationally significantly more efficient than the direct evaluation of the volume or surface integrals (36). B.3. Displacement jump ! For the computation of the effective microscale displacement jump ½½ v M , Eq. (41) is next evaluated in a discrete FE setting nnodes -n 1X - nn ! an DBn0m ð um D Xm  ðF bM *Þc Þ ½½ v M  ¼   S n m

ðB:11Þ

where nnodes is the number of nodes in the mesh; anm is a coefficient following from the discretized form of the integral, e.g. containing the effective area associated with node n. In (B.11) it has been taken into account that the displacement of the - nn reference node has been fixed, which allows to omit D in front of um . The volume integral in S (Eq. (42)) can be evaluated  either analytically or numerically. ! For the subsequent calculation of the consistent tangent associated to ½½ v M  and F bM , it is convenient to eliminate the contribution of the dependent nodes from (B.11). The nodes nnodes in the mesh can be subdivided into the dependent boundary nodes d, the reference node r and the remaining ‘free’ nodes f. The free nodes f include the independent boundary nodes and all interior nodes. The dependent degrees of freedom can be eliminated from Eq. (B.11), which results in -f 1 - nf ! ½½ v M  ¼ MT ð um D Xm  ðF bMn Þc Þ  S e 

ðB:12Þ

1506

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

where the element I of the column M is computed as e dX nodes J fMgI ¼ aIm DBI0m þ aJm DB0m fDgJI

ðB:13Þ

J¼1 - nd

- nf

- nf







with fDgJI the components of the extended dependency matrix defined as um ¼ D um ¼ ½D,0 um , with D the dependency matrix relating the boundary nodes only. Note that in (B.12) the contribution of the reference node has been omitted, since ! !r !nr ! u m ¼ 0 and D X m ¼ 0 .

B.4. Macroscopic tangents To iteratively solve the macroscopic boundary value problem, the consistent tangents of the macroscopic stress PM and ! the effective microscale displacement jump ½½ v M  with respect to the macroscopic deformation F bM are required. The consistent tangents 4 CM and 3 DM are defined in Eqs. (18) and (19), respectively. The tangents can be obtained numerically by perturbation of the components of F bM . Computationally more efficient, however, is to determine the tangents by condensation of the stiffness matrix resulting from the converged microscopic boundary value problem. To this aim, the ! linearized microscale system of equations relating the variation of nodal displacements d u to the variation of external ! forces d f is partitioned as 2 - 3 2 3 2 - 3 0 f d ff 2 3 d um 6 m 7 6 e 7 K ff K fa K fr 6  7 6 7 6 - 7 6 -a 7 6K 7 6 - 7 7 ¼ 6d fa 7 6 d fM 7 ðB:14Þ 4 af K aa K ar 5  6 7 d a 6 M 7 6 M 7 6  7 4 5 6  7 6 4 K rf K ra K rr 4 !r 5 -r !r 5 d um d fm df m

The matrices K contain second-order tensors relating nodal displacements to nodal external forces. Note that, in this linearized system of equations the dependent degrees of freedom d have already been eliminated. The system (B.14) is -f

evaluated in a converged state at the end of an increment. Therefore, the variation of the forces acting on the free nodes f m 

are close to zero within the convergence tolerance. -f Condensing d um out of the system (B.14), making use of (B.5) and substituting the result in the variation of the  macroscopic stress (B.10), allows the identification of the expression for the constitutive tangent 4

CM ¼

1 ðð AM ÞT  K aa  AM ÞLC  V 0m  %

ðB:15Þ

with K aa ¼ K aa K af  ðK ff Þ1  K fa %

ðB:16Þ

which are the same expressions as those used in the first-order computational homogenization. ! To derive the constitutive tangent associated with the displacement jump ½½ v M , the linearized system (B.14) is exploited as well. From this system together with 2 - 3 2 - 3 da A 6  M 7 6 M 7 b ðB:17Þ 4 !r 5 ¼ 4 !r 5  ðdF M Þc dum D Xm - nf

the variations of the displacement of the free nodes d um is derived  2 - 3 A - nf 6 M 7 d um ¼ ðK ff Þ1  ½K fa K fr   4 : ðdF bM Þc !r 5  D Xm

ðB:18Þ

! !r Substitution of this linear relation in the variational form of Eq. (B.12), exploiting that D X m ¼ 0 and taking into account a consistent sequencing of all arrays and matrices yields the following expression: -f 1 ! d½½ v M  ¼  MT ððK ff Þ1  K fa  AM þ I  D Xm Þ : ðdF bM Þc   S e

ðB:19Þ

E.W.C. Coenen et al. / J. Mech. Phys. Solids 60 (2012) 1486–1507

1507

Comparing the result of the latter Eq. (B.19) with Eq. (19)results in the following expression for the constitutive tangent: 3

DM ¼ 

-f 1 T M ððK ff Þ1  K fa  AM þ I  D Xm Þ   S e

ðB:20Þ

References Armero, F., Garikipati, K., 1996. An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids. Int. J. Solids Struct. 33, 2863–2885. Belytschko, T., Loehnert, S., Song, J.H., 2008. Multiscale aggregating discontinuities: a method for circumventing loss of material stability. Int. J. Numer. Methods Eng. 73, 869–894. de Borst, R., 1987. Computation of post-bifurcation and post-failure behavior of strain-softening solids. Comput. Struct. 25, 211–224. Coenen, E.W.C., Kouznetsova, V.G., Geers, M.G.D., 2011. Enabling microstructure-based damage and localization analyses and upscaling. Modelling Simulation Mater. Sci. Eng. 19, 074008. Coenen, E.W.C., Kouznetsova, V.G., Geers, M.G.D., 2012. Novel boundary conditions for strain localization analysis in microstructural volume elements. Int. J. Numer. Methods Eng. 90, 1–21. Farhat, C., Pierson, K., Lesoinne, M., 2000. The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems. Comput. Methods Appl. Mech. Eng. 184, 333–374. Feyel, F., Chaboche, J.L., 2000. FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fiber SiC/Ti composite materials. Comput. Methods Appl. Mech. Eng. 183, 309–330. Fish, J., Suvorov, A., Belsky, V., 1997. Hierarchical composite grid method for global–local analysis of laminated composite shells. Appl. Numer. Math. 23, 241–258. Luciano da Fontoura Costa, Roberto Marcondes Cesar, Jr., 2009. Shape Classification and Analysis: Theory and Practice, CRC Press, London. Ghosh, S., Lee, K., Moorthy, S., 1995. Multiple scale analysis of heterogeneous elastic structures using homogenization theory and Voronoi cell finite element method. Int. J. Solids Struct. 32, 27–62. Gitman, I.M., Askes, H., Sluys, L.J., 2008. Coupled-volume multi-scale modelling of quasi-brittle material. Eur. J. Mech.—A/Solids 27, 302–327. Guidault, P., Allix, O., Champaney, L., Cornuault, C., 2008. A multiscale extended finite element method for crack propagation. Comput. Methods Appl. Mech. Eng. 197, 381–399. Hill, R., 1962. Acceleration waves in solids. J. Mech. Phys. Solids 10, 1–16. Hill, R., 1963. Elastic properties of reinforced solids: some theoretical principles. J. Mech. Phys. Solids 11, 357–372. Hirschberger, C.B., Ricker, S., Steinmann, P., Sukumar, N., 2009. Computational multiscale modelling of heterogeneous material layers. Eng. Fract. Mech. 76, 793–812. Hough, P.V.C., 1959. Machine analysis of bubble chamber pictures, In: Proceedings of the International Conference on High Energy Accelerators and Instrumentation. Hund, A., Ramm, E., 2007. Locality constraints within multiscale model for non-linear material behaviour. Int. J. Numer. Methods Eng. 70, 1613–1632. Kouznetsova, V.G., Brekelmans, W.A.M., Baaijens, F.P.T., 2001. An approach to micro–macro modeling of heterogeneous materials. Comput. Mech. 27, 37–48. Kouznetsova, V.G., Geers, M.G.D., Brekelmans, W.A.M., 2002. Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme. Int. J. Numer. Methods Eng. 54, 1235–1260. Kouznetsova, V.G., Geers, M.G.D., Brekelmans, W.A.M., 2004. Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy. Comput. Methods Appl. Mech. Eng. 193, 5525–5550. Larsson, R., Zhang, Y., 2007. Homogenization of microsystem interconnects based on micropolar theory and discontinuous kinematics. J. Mech. Phys. Solids 55, 819–841. Linder, C., Armero, F., 2007. Finite elements with embedded strong discontinuities for the modeling of failure in solids. Int. J. Numer. Methods Eng. 72, 1391–1433. Loehnert, S., Belytschko, T., 2007. A multiscale projection method for macro/microcrack simulations. Int. J. Numer. Methods Eng. 71, 1466–1482. Massart, T.J., Peerlings, R.H.J., Geers, M.G.D., 2007. An enhanced multi-scale approach for masonry wall computations with localization of damage. Int. J. Numer. Methods Eng. 69, 1022–1059. Matous, K., Kulkarni, M.G., Geubelle, P.H., 2008. Multiscale cohesive failure modeling of heterogeneous adhesives. J. Mech. Phys. Solids 56, 1511–1533. Mercatoris, B.C.N., Massart, T.J., 2011. A coupled two-scale computational scheme for the failure of periodic quasi-brittle thin planar shells and its application to masonry. Int. J. Numer. Methods Eng. 85, 1177–1206. Miehe, C., Dettmar, J., Zah, D., 2010. Homogenization and two-scale simulations of granular materials for different microstructural constraints. Int. J. Numer. Methods Eng. 83, 1206–1236. Nguyen, V.P., Lloberas-Valls, O., Stroeven, M., Sluys, L.J., 2011. Homogenization-based multiscale crack modelling: from micro-diffusive damage to macrocracks. Comput. Methods Appl. Mech. Eng. 200, 1220–1236. Nguyen, V.P., Stroeven, M., Sluys, L.J., 2012. An enhanced continuous.discontinuous multiscale method for modeling mode-I cohesive failure in random heterogeneous quasi-brittle materials. Eng. Fract. Mech. 79, 78–102. Rice, J.R., 1976. The localization of plastic deformation. In: Koiter, W.T. (Ed.), Theoretical and Applied Mechanics, North-Holland Publishing Company, pp. 207–220. Simo, J.C., 1992. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Comput. Methods Appl. Mech. Eng. 99, 61–112. Simo, J.C., Oliver, J., Armero, F., 1993. An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. Comput. Mech. 12, 277–296. Souza, F.V., Allen, D.H., 2011. Modeling the transition of microcracks into macrocracks in heterogeneous viscoelastic media using a two-way coupled multiscale model. Int. J. Solids Struct. 48, 3160–3175. Suquet, P.M., 1985. Local and global aspects in the mathematical theory of plasticity. In: Sawczuck, A., Bianchi, G. (Eds.), Plasticity Today: Modelling, Methods and Applications, Elsevier Applied Science Publishers, London, pp. 279–310. Verhoosel, C.V., Remmers, J.J.C., Gutie´rrez, M.A., 2010. A partition of unity-based multiscale approach for modelling fracture in piezoelectric ceramics. Int. J. Numer. Methods Eng. 82, 966–994. Wells, G.N., Sluys, L.J., 2000. Application of embedded discontinuities for softening solids. Eng. Fract. Mech. 65, 263–281. Zohdi, T., Wriggers, P., 1999. A domain decomposition method for bodies with heterogeneous microstructure based on material regularization. Int. J. Solids Struct. 36, 2507–2525.