A micromorphic computational homogenization framework for heterogeneous materials

Accepted Manuscript

A micromorphic computational homogenization framework for heterogeneous materials R. Biswas, L.H. Poh PII: DOI: Reference:

S0022-5096(16)30895-X 10.1016/j.jmps.2017.02.012 MPS 3069

To appear in:

Journal of the Mechanics and Physics of Solids

Received date: Revised date: Accepted date:

8 December 2016 26 February 2017 27 February 2017

Please cite this article as: R. Biswas, L.H. Poh, A micromorphic computational homogenization framework for heterogeneous materials, Journal of the Mechanics and Physics of Solids (2017), doi: 10.1016/j.jmps.2017.02.012

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

R. Biswas, L. H. Poh∗

CR IP T

A micromorphic computational homogenization framework for heterogeneous materials Department of Civil and Environmental Engineering, National University of Singapore 1 Engineering Drive 2, E1A-07-03, Singapore 117576

AN US

Abstract

AC

CE

PT

ED

M

The conventional first-order computational homogenization framework is restricted to problems where the macro characteristic length scale is much larger than the underlying Representative Volume Element (RVE). In the absence of a clear separation of length scales, higher-order enrichment is required to capture the influence of the underlying rapid fluctuations, otherwise neglected in the first-order framework. In this contribution, focusing on matrix-inclusion composites, a novel computational homogenization framework is proposed such that standard continuum models at the micro-scale translate onto the macro-scale to recover a micromorphic continuum. Departing from the conventional FE2 framework where a macroscopic strain tensor characterizes the average deformation within the RVE, our formulation introduces an additional macro kinematic field to characterize the average strain in the inclusions. The two macro kinematic fields, each characterizing a particular aspect of deformation within the RVE, thus provide critical information on the underlying rapid fluctuations. The net effect of these fluctuations, as well as the interactions between RVEs, are next incorporated naturally into the macroscopic virtual power statement through the Hill-Mandel condition. The excellent predictive capability of the proposed homogenization framework is illustrated through three benchmark examples. It is shown that the homogenized micromorphic model adequately captures the material responses, even in the absence of a clear separation of length scales between macro and micro. ∗

Corresponding author. Telephone: +65 6516 4913. Email address: [email protected] (L. H. Poh)

Preprint submitted to Elsevier

March 2, 2017

ACCEPTED MANUSCRIPT

CR IP T

Keywords: Computational homogenization, Multi-scale modelling, Micromorphic continuum, Higher order continuum, Matrix-inclusion composites 1. Introduction

CE

PT

ED

M

AN US

Multi-phase or heterogeneous materials boast enhanced mechanical properties compared to their respective constituents, and have been the subject of increasing interest in recent years. The mechanical response of a heterogeneous material is highly dependent on the microstructural characteristics e.g. morphology, spatial distribution, material and interface properties. A predictive analysis of such materials can be achieved by performing direct numerical simulations (DNS), where the microstructural details are explicitly modelled with a mesh discretization of sufficient resolution. However, such comprehensive analyses are seldom practical for a typical engineering problem due to their inherent high computational costs. Motivated by the need for predictive models with a more amenable computational cost, the first-order computational homogenization approach was proposed (Feyel and Chaboche, 2000; Ghosh et al., 1995; Kouznetsova et al., 2001; Matsui et al., 2004; Miehe et al., 2010; Smit et al., 1998; Suquet, 1985; Terada and Kikuchi, 2001). Without making any constitutive assumptions at the macro-scale, the homogenized material response is determined on the fly by solving a micro-scale Boundary Value Problem (BVP) associated with each macroscopic material point. Arbitrary non-linear and history-dependent material behaviour of the micro-constituents can thus be incorporated naturally into the framework. The first-order computational homogenization approach is built upon a separation of scales assumption. This assumption, however, is violated in several situations, for instance, when

AC

• the characteristic structural size is of the same order of magnitude as the underlying microstructure, e.g. in Micro-Electro-Mechanical Systems (MEMS); • the characteristic loading wavelength spans across only a few microstructural cells, e.g. in micro-indentation experiments;

2

ACCEPTED MANUSCRIPT

• the macroscopic deformation localizes into narrow bands that transverse across only a few microstructural cells, e.g. during strain softening.

AC

CE

PT

ED

M

AN US

CR IP T

In the absence of a clear separation of length scales, higher-order enrichment is required to capture the propagation of the underlying rapidly fluctuating responses onto the macro-scale. To this end, a second-order computational homogenization method was proposed by Kouznetsova et al. (2004, 2002), where application of the Hill-Mandel condition results in a second-gradient continuum at the macro-scale. Therewith, a length scale parameter associated with the second gradient of displacement characterizes the size of an underlying RVE, through which the micro-fluctuation effect is captured and translated to the macro-scale. Kaczmarczyk et al. (2008) later extended the second-gradient framework to incorporate traction and displacement boundary conditions on the RVEs, along with the periodic ones in a generalized unified manner. An alternative second-gradient framework is proposed by Bacigalupo and Gambarotta (2010), with an appropriate representation of the micro-displacement field to ensure its continuity across the interfaces between adjacent unit cells. In general, the numerical framework for the second-gradient approach requires C 1 continuity. This requirement can be relaxed to C 0 continuity by adopting mixed formulations with additional degrees of freedom involving an assumed macro strain field, as well as a Lagrange multiplier enforcing the kinematic constraints (e.g. Kouznetsova et al., 2004, 2002). Note that the latter can be removed via a static condensation procedure. Alternatively, the Discontinuous Galerkin (DG) method can be adopted at the macro-scale to weakly impose the C 1 continuity requirement (e.g. Nguyen et al., 2013; Nguyen and Noels, 2014). The generalized micromorphic theory (Eringen, 1968; Forest, 2009, 2016; Germain, 1973) constitutes another broad class of higher-order continua, where additional morphic kinematic fields are introduced to characterize the energetic contribution of the rapidly fluctuating micro-processes, otherwise not picked up in a standard model. In the most general formulation, the morphic variable involved is a second-rank non-symmetric tensor, from which other higher-order theories can be recovered as special cases (Forest, 2009; Forest and Sievert, 2003, 2006) e.g. a Cosserat continuum is retrieved by constraining the morphic field to be purely rotational; the second gradient theory is obtained by setting the morphic field as the deformation gradient via a penalty constraint. 3

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

In literature, there are several computational homogenization methods developed within the generalized continuum framework. Feyel (2003) adopted the Cosserat homogenization framework proposed by Forest and Sab (1998), to develop a bottom-up technique for fiber-matrix composite systems. Larsson and Diebels (2007) reformulated the second-order computational homogenization framework by placing a constraint on the second-gradient deformation tensor in terms of the micropolar stretch, which leads to a micropolar continuum. J¨anicke and Diebels (2010); J¨anicke et al. (2009) implemented a numerical homogenization scheme incorporating a full micromorphic continuum at the macro-scale which captures the size-dependent boundary layer effect in cellular materials. However, the associated kinematic constraints could not be prescribed solely at the RVE boundary, since the volume integral characterizing the morphic variable cannot be readily transformed into a corresponding surface integral. To this end, J¨anicke and Steeb (2012) proposed the mimimal loading conditions to impose the macro kinematic variables onto the RVE in terms of volume average constraints. A critical review of the non-homogeneous boundary conditions adopted in the generalized computational homogenization approaches can be found in Forest and Trinh (2011). Recently, a homogenization strategy was presented in H¨ utter (2016) such that a micromorphic continuum was recovered at the macro-scale based on a Cauchy micro-continuum. Therewith, the equilibrium condition at the micro-scale has to be appended with additional (non-standard) body forces, which propagate onto the macro-scale as generalized morphic stresses. In this contribution, a novel computational homogenization framework is proposed to obtain the macroscopic response of matrix-inclusion composites via a micromorphic continuum. The proposed approach addresses the limitations of a conventional FE2 framework, and avoids the deficiencies of existing second-gradient schemes, i.e. a dependency on the choice of RVE and an over-constrained fluctuation field as highlighted in Forest and Trinh (2011). Herewith, departing from the second-gradient frameworks, an additional macro kinematic field is introduced to characterize the average deformation in the inclusions. As per the conventional FE2 framework, the macro-distortion tensor is defined as the average distortion within the RVE. The underlying rapidly fluctuating response is thus captured by these two macro kinematic fields, which is next propagated onto the macro-scale via the Hill-Mandel condition. The resulting formulation at the macro-scale recovers a micromorphic continuum, instead of being postulated a priori. The length scale parameter associated with the higher-order term thus characterizes the 4

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

nonlocal interaction between neighbouring micro-mechanisms within a RVE, which in turn provides a regularization effect and enables an accurate prediction of the size-effect. An elasto-plastic behaviour can be adopted for the micro-constituents easily, in contrast to some existing approaches restricted to elastic analyses. The numerical implementation of the proposed framework is straightforward since it only requires C 0 continuity in the primary variables. Furthermore, a clear physical interpretation is attached to the morphic variable, which mostly remains elusive in a conventional micromorphic model. Note also that our framework starts with standard continuum models at the micro-scale, a departure from another class of homogenization theories where higher-order formulations are utilized at the small-scale, e.g. in Li et al. (2014) where a macro gradient Cosserat continuum is obtained based on a (Cosserat continuum equivalent) discrete particle assembly at the meso-scale; the micromorphic continuum in Poh (2013) based on the homogenization of a meso-scale gradient plasticity model; the computational homogenization strategy in Hirschberger et al. (2008) where a micromorphic continuum at the micro-scale translates onto the macro-scale as the interfacial response of a material layer. The paper is organized as follows. In Section 2, the micromorphic computational homogenization strategy is outlined. A decomposition of micro kinematic fields in terms of the macro-scale variables is proposed in Section 3, which is then utilized in the Hill-Mandel condition to derive the governing equations leading to a micromorphic continuum in Section 4. Section 5 describes the details of the microstructural modelling. Finally, the excellent predictive capability of the proposed homogenization framework is illustrated in Section 6 with three numerical examples. Indicial notations are mostly used for clarity of presentation. To perform the scale transition consistently, a two-scale coordinate system is adopted with x and y denoting the macro and micro coordinates respectively. To keep track of quantities at the two scales, those existing at the micro-scale are denoted with a ˆ symbol. Differentiation with respect to macro and ˆ ∂A ∂A micro coordinates are written as ∂xijk = ∇k Aij and ∂ykij = Aˆij,k respectively. R The volume average over the RVE is denoted as hˆ·i = ˆˆ1 Vˆ ˆ· dVˆ , where ˆ· V represents a generic field at the micro-scale.

5

ACCEPTED MANUSCRIPT

2. Homogenization strategy

CR IP T

In this contribution, we focus on composite materials with soft inclusions. During deformation, a polarizing effect emanates from each soft inclusion to induce a rapidly fluctuating response within a RVE. To adequately capture the overall interactions between the micro-mechanisms, and their propagation to the macro-scale, a bottom-up strategy is proposed such that standard formulations at the micro-scale translate consistently onto the macro-scale to recover a micromorphic continuum. The homogenization framework is based on the concurrent coupling of micro and macro-scale BVPs, as depicted in Fig 1, which entails the following:

AC

CE

PT

ED

M

AN US

(a) As per the standard homogenization approach, a distortion tensor Hij at the macro-scale captures the average distortion within the RVE. Note that the distortion tensor is non-symmetric in general since it corresponds to the displacement gradient. An additional macro kinematic e ij is furthermore introduced to characterize the average distortion field H within the inclusions of a RVE. The micro-scale kinematic fields are next e ij . decomposed in terms of Hij and H (b) Applying the Hill-Mandel condition to determine how the micro-stresses, arising from different mechanisms within a RVE, propagate and manifest themselves at the macro-scale. The governing equations for the resulting micromorphic continuum are next derived from the principle of virtual power. e ij ) from macro to micro, in (c) Translating the kinematic fields (Hij and H order to define the BVP of a RVE with the associated constraint conditions. The relevant macro stresses and tangent operators are extracted directly from the RVE, which are translated back to the macro-scale to complete the coupling between the two scales.

6

ACCEPTED MANUSCRIPT

Decomposition of kinematic fields in terms of

(b)

and

Deriving the macro-scale governing equations Macro-scale BVP

tangents

AN US

(c)

CR IP T

(a)

M

Micro-scale BVP

Figure 1: Micromorphic computational homogenization strategy.

ED

3. Decomposition of kinematic fields

PT

The first step of the homogenization framework in Fig.1 entails the definition of two macro kinematic fields : • A macro distortion tensor Hij to characterize the average deformation within the RVE as per standard homogenization approach.

AC

CE

e ij to characterize the average defor• An additional (morphic) variable H mation of the inclusions within a RVE, which in turn captures the fluctuating response induced by the polarizing effect of the softer phase. e ij characterizes a different deformation mechanism comNote that H pared to Hij , thus enabling the framework to extract additional information on the underlying processes, in the spirit of a generalized micromorphic continuum.

The homogenization framework requires a decomposition of the micro kinematic variables into slow (macro) and fast (micro) varying fields. To 7

ACCEPTED MANUSCRIPT

Macro continuum

x2

x x1 y2

y1

AN US

x

CR IP T

this end, a two-scale coordinate system is adopted with x and y representing the macro and micro coordinates respectively, as illustrated in Fig. 2. For clarity of presentation, the kinematic fields are written together with their arguments in this section.

Figure 2: A RVE centered about a macroscopic point x.

PT

ED

M

Since the rapid fluctuations within the RVE are mostly induced by the softer inclusions, we decompose the micro kinematic fields in terms of the e ij (x) and its macro spatial gradient, superimposed with a slow varying H rapidly fluctuating field uˆ∗i (x, y) that captures the extent of heterogeneities at the micro-scale o 1n e e e ∇k Hij (x) + ∇k Hji (x) yj yk + uˆ∗i (x, y) , (1) uˆi (x, y) = Hij (x)yj + 4 where uˆi (x, y) denotes the micro displacement field. The micro distortion ˆ ij (x, y) is next obtained by differentiating (1) with respect to y field H

AC

CE

ˆ ij (x, y) = uˆi,j (x, y) H

n o e ij (x) + ∇k H e ji (x) + ∇j H e ik (x) + ∇j H e ki (x) yk e ij (x) + 1 ∇k H =H 4 (2) + uˆ∗i,j (x, y) .

The macro-micro transition is effected via kinematic constraints associe ij (x), which are defined as ated with the macro variables Hij (x) and H Z 1 ˆ ij (x, y)dVˆ , Hij (x) = H (3) Vˆ Vˆ 8

ACCEPTED MANUSCRIPT

Z

o 1n e e ji (x) + ∇j H e ik (x) + ∇j H e ki (x) y¯k |I , ∇k Hij (x) + ∇k H uˆ∗i,j (x, y)dVˆ = − 4 VˆI (6)

AN US

1 VˆI

CR IP T

Z 1 e ˆ ij (x, y)dVˆ , Hij (x) = H (4) VˆI VˆI with Vˆ and VˆI representing the RVE and inclusions volume respectively. e ij (x) characterize different aspects of the same local Note that Hij (x) and H ˆ ij (x, y) within a RVE. Substituting (2) into (3) and (4), the distortion field H following kinematic constraints are obtained Z 1 e ij (x) , (5) uˆ∗i,j (x, y)dVˆ = Hij (x) − H ˆ V Vˆ R where the fact that Vˆ y dVˆ = 0 has been utilized. Furthermore,

PT

ED

M

where y¯k |I denotes the overall centroid of the inclusions about the RVE origin. Note from (5) and (6) that the average rapid fluctuation within the RVE is characterized through the two macro kinematic fields, while a non-local interaction induced by the polarizing effect of inclusions is incorporated via e ij (x). the macro spatial gradient of H To facilitate subsequent derivations and drawing reference to (5), the micro fluctuation field is rewritten as   e kl (x) , uˆ∗i,j (x, y) = Pˆijkl (y) Hkl (x) − H (7)

AC

CE

where Pˆijkl (y) is a polarization tensor characterizing the rapid fluctuations. In general, the polarization tensor associated with an elasto-plastic composite is history and state dependent, which is extracted numerically from the RVE at each deformation step. This is in contrast to the state-independent 4th order elastic strain concentration tensor utilized in a linear homogenization theory (e.g. Klusemann and Svendsen, 2010; Pierard et al., 2004). In view of (5), a constraint on the polarization tensor is thus Z 1 Pˆijkl (y)dVˆ = δik δjl . (8) Vˆ Vˆ

9

ACCEPTED MANUSCRIPT

4. Homogenization towards a micromorphic continuum

CR IP T

We next proceed to step (b) of the homogenization strategy in Fig.1, where the homogenized model is determined through the Hill-Mandel condition which imposes an equivalence of power density across the two scales. Correspondingly, the power density at a macroscopic point P den (x) is given by the average power expended within the underlying RVE P den (x) = hPˆ (x, y)i .

(9)

AN US

Since a standard continuum description is adopted at the micro-scale, the average internal power expended within a RVE is given by Z 1 ˆ˙ ij dVˆ ˆ (10) σ ˆij H hP (x, y)i = ˆ ˆ V V where σ ˆij denotes the micro stress field. Utilizing the kinematic decomposition (2) and polarization tensor defined in (7), the power density at a macro point is obtained as

ij

ED

M

˜˙ ij + 1 hˆ ˜˙ ij + ∇k H ˜˙ ji + ∇j H ˜˙ ik + ∇j H ˜˙ ki ) + hˆ P den (x) = hˆ σij iH σij yk i(∇k H σij uˆ˙ ∗i,j i 4 (11a) ˜˙ , ˜˙ + ζ ∇ H (11b) = σ H˙ + ξ H ij

ij

ij

ijk

k

ij

PT

where use has been made of (9). The stress tensors in (11b) can be identified as σij = hˆ σkl Pˆklij i ,

(12a)

CE

ξij = hˆ σij i − hˆ σkl Pˆklij i , 1 ζijk = hˆ σij yk + σ ˆji yk + σ ˆik yj + σ ˆjk yi i , 4

(12b) (12c)

AC

with σij denoting the macro stress, ξij and ζijk the generalized stresses emanating from interactions between the two phases within a RVE. Note the departure from many homogenization approaches that define the macro stress as the average micro stress within a RVE. In our formulation, e ij captures the redistributhe kinematic decomposition corresponding to H tion of stress field at the micro-scale. Through the Hill-Mandel condition, this 10

ACCEPTED MANUSCRIPT

AN US

CR IP T

polarizing effect propagates and manifests itself in σij via (12a). Note also that the work contribution by the fast fluctuating field uˆ∗i in (11a) is incorporated into the formulation through the polarization tensor. In contrast, the average work done by the rapid fluctuations is assumed to vanish in most computational homogenization frameworks, including those enhanced with higher order terms (e.g. Kouznetsova et al., 2004). It is furthermore highlighted that (11b) recovers a micromorphic continuum (e.g. Eringen, 1968; Forest, 2009, 2016; Forest and Sievert, 2003) based on a RVE description with standard continuum models, a departure from many enhanced homogenization methods where higher order terms are postulated a priori at the macro-scale. We next consider an arbitrary macroscopic region V which is much larger than the RVE size, i.e. V >> l. The total power expended in V can thus be obtained by integrating the macro power density (11a) in space, Z e˙ ij + ζijk ∇k H e˙ ij )dV Pint = (σij H˙ ij + ξij H VZ Z ˙ ˙ e e e˙ ij )dS , = − (σij,j u˙ i − ξij H ij + ζijk,k H ij )dV + (σij nj u˙ i + ζijk nk H

M

V

∂V

(13)

PT

ED

where ∂V represents the external boundary of the considered domain, with nj as the corresponding outward unit normal. Ignoring body forces, the power expended by the external forces is postulated as Z e˙ ij )dS , Pext = (ti u˙ i + mij H (14) ∂V

AC

CE

where ti and mij are the standard and higher order tractions respectively. The method of virtual power (Germain, 1973) is next applied to derive the governing equations for the micromorphic continuum. By imposing Pint = Pext for arbitrary kinematic fields, a system of governing equations is obtained as ∇j σij = 0i , ξij = ∇k ζijk ,

(15a) (15b)

where (15a) is the standard equilibrium condition, while (15b) is a homoge11

ACCEPTED MANUSCRIPT

nized micro-force balance. The tractions acting at the external surfaces are also determined as (16a) (16b)

CR IP T

ti = σij nj , mij = ζijk nk .

4.1. Comparison with second gradient computational homogenization Referring to the definitions of generalized stress tensors in (12b) and (12c), the homogenized micro-force balance (15b) can be written as

(17)

AN US

1 hˆ σkl Pˆklij i = hˆ σij i − ∇k hˆ σij yk + σ ˆji yk + σ ˆik yj + σ ˆjk yi i , 4

ED

M

which describes the interactions between stress quantities arising from different micro-mechanisms. Specifically, the average polarized stress characterizing the rapidly fluctuating deformation within a RVE has to be balanced by the corresponding average Cauchy stress, augmented with an average moment stress characterizing the nonlocal interactions between the two phases. Substituting (17) into (15a) the macro equilibrium condition reduces to   1 σij yk + σ ˆji yk + σ ˆik yj + σ ˆjk yi i = 0i , (18) ∇j hˆ σij i − ∇k hˆ 4

AC

CE

PT

which has the same form as that in the second-order computational homogenization framework (Kouznetsova et al., 2004). Note that the resemblance between the two higher-order continua has also been reported in literature, i.e. the balance equations of a micromorphic continuum can be combined to recover the governing equation of a strain gradient continuum, independent of any homogenization schemes (Forest and Sievert, 2003). In our micromorphic approach, where the (higher-order) equilibrium condition is decomposed into a system of balance equations, the resulting numerical implementation requires only C 0 continuity (see Appendix A). Our formulation thus can be implemented more easily, compared to the C 1 continuity requirement, or C 0 continuity with additional numerical treatments as in the second-order computational homogenization framework. Moreover, despite the similarity in the governing equation, subtle differences exist between the two formulations, e.g. a non-vanishing contribution to the homogenized power statement by the rapidly fluctuating processes in (11a). 12

ACCEPTED MANUSCRIPT

5. Multiscale computational homogenization framework

CR IP T

To complete the homogenization framework, we now proceed to step (c) of the strategy listed in Section 2. Without making any assumptions at the macro-scale, the relevant material responses are extracted from the BVP solution of the RVE embedded at each macro material point. The multiscale framework thus requires a consistent transfer of information between the two scales, as highlighted in Fig. 1. In this section, the coupling between the macro and micro scales are expressed in a form that facilitates subsequent finite element implementation.

AN US

5.1. Macro-micro transition Focusing on matrix-inclusion composites, a square RVE of size l is considered within the context of a 2D analysis, as depicted in Fig. 3. Standard continuum models are adopted at the micro-scale to describe the behaviour of the RVE. In the absence of body forces, the state of equilibrium is given by

ED

2

M

σ ˆij,j = 0i .

Top

3

y2

Left

Right

PT CE AC

(19)

l

y1

1

4

Bottom l

Figure 3: Schematic representation of a RVE. Symbols × and O indicate analogous points at the opposite surfaces of the RVE. The four master nodes are labelled as 1,2,3,4.

The constituent behavior is described with arbitrary history dependent constitutive laws. The macro-micro transition is achieved by enforcing the 13

ACCEPTED MANUSCRIPT

CR IP T

kinematic constraints discussed in Section 3. As the micro-scale BVP is solved using standard finite element method, constraints (5) and (6) are imposed in terms of the primary variable uˆi , i.e. the continuous displacement field within the RVE. Applying the divergence theorem, (5) is rewritten as Z 1 e ij uˆ∗i n ˆ j dSˆ = Hij − H (20) Vˆ ∂ Vˆ

AN US

where ∂ Vˆ and n ˆ j represent the external boundary of the RVE and the corresponding outward normal respectively, as shown in Fig. 3. In literature, the periodic boundary condition has been widely adopted for the constraint involving fluctuation field (e.g. Kouznetsova et al., 2004), which have been reported to provide better estimates of the effective properties (Kaczmarczyk et al., 2008, 2010). Here, a generalized periodic boundary condition is adopted to effect (20), with the fluctuation field at analogous points on the opposite edges of the RVE boundary related according to

(21a) (21b)

M

R L e uˆ∗R ˆ∗L i −u i = (Hij − Hij )(yj − yj ) , T B e uˆ∗T ˆ∗B i −u i = (Hij − Hij )(yj − yj ) ,

ED

where the superscripts L, R, T and B denote the left, right, top and bottom edges of the RVE respectively. Drawing reference from the decomposition in (1), (21a) and (21b) are rewritten as

CE

PT

1 e ij + ∇k H e ji )(yjR − yjL )(ykR − ykL ) , (22a) uˆR ˆLi = Hij (yjR − yjL ) + (∇k H i −u 4 1 T B T B T B e e uˆTi − uˆB i = Hij (yj − yj ) + (∇k Hij + ∇k Hji )(yj − yj )(yk − yk ) . (22b) 4

To enforce these constraints within the finite element framework, displacements at the four corner (master) nodes are prescribed as

AC

1 (c) (c) e ij + ∇k H e ji )y (c) y (c) uˆi = Hij yj + (∇k H j k 4

for c = 1, 2, 3, 4

(23)

with superscript c denoting the corner node location as depicted in Fig. 3.

14

ACCEPTED MANUSCRIPT

Rearranging the terms, (22a) and (22b) reduce to (24a) (24b)

CR IP T

1 1 (4) (1) (3) (2) ˆLi = (1 − h)(ˆ uˆR ui − uˆi ) + (1 + h)(ˆ ui − uˆi ) , i −u 2 2 1 1 (2) (1) (3) (4) uˆTi − uˆB ui − uˆi ) + (1 + g)(ˆ ui − uˆi ) , i = (1 − g)(ˆ 2 2

where ∂ Iˆ =

r S

AN US

with g = 2yl 1 , h = 2yl 2 representing the normalized coordinates of the pair of analogous points at the top-bottom and left-right boundaries respectively. The constraint in (6) is next written in terms of the displacement at the matrix-inclusion interfaces as Z 1 e ij , uˆi n ˆ Ij dS = H (25) ˆ VI ∂ Iˆ ∂ Iˆk , with r denoting the total number of inclusions inside the

k=1

k=1

e ij , Wjk uˆki = H

(26)

ED

N X

M

RVE. The outward normal of an inclusion is denoted by n ˆ Ij . Numerically, (25) is discretized as

CE

PT

where Wj is a N × 1 vector of coefficients that account for the weighted contribution of each boundary node towards the surface integral in (25) for e ij (Larsson and Runesson, 2007; Miehe and Koch, 2002), each component of H ˆ with N as the total number of nodes along ∂ I. The macro-micro transition is thus realized through (23),(24), and (26). Together with (19), and the constitutive laws for the micro-constituents, the micro-scale BVP is fully described.

AC

5.2. Micro-macro transition After the micro-scale BVP is solved, the material response propagates to the macro-scale via the macro-stress (σij ) and generalized stresses (ξij , ζijk ). For computational efficiency and ease of numerical implementation, the volume integrals defined in (12a)-(12c) are converted to surface integrals using the divergence theorem (see Appendix A.4), leading to the following

15

ACCEPTED MANUSCRIPT

expressions Z 1 σij = tˆi yj dSˆ , ˆ ˆ V ∂I VˆI ˜ ξij = λij , Vˆ Z
1 tˆi yj yk + tˆj yi yk dSˆ , ζijk = 4Vˆ ∂ Iˆ

CR IP T

(27a) (27b)

(27c)

AN US

˜ ij with tˆi denoting the traction on the external boundaries of the RVE, and λ representing the Lagrange multiplier associated with the kinematic constraint ˜ ij characterizes the redistribution of stresses within the in (6). Physically, λ RVE resulting from the nonlocal interactions among inclusions, as well as between the two phases. The final step in the multi-scale framework involves the computation of tangent operators for the micromorphic continuum. The macro-scale stress increments are linearized in terms of the kinematic field variables as e kl + C e dσij = Cijkl dHkl + Cijkl dH ijklm d∇m Hkl , (3) (4) e (2) e kl , dξij = C dHkl + C dH d∇m H kl + C (1)

ijkl ijkl ijklm (3) (4) e lm + Cijklmn d∇n H e lm Cijklm dHlm + Cijklm dH

(28a) (28b)

,

(28c)

ED

dζijk =

(2)

M

(1)

PT

where C (i) , C(j) , and C (i, j = 1, 2, 3, 4) represent the 4th , 5th , and 6th order tangent stiffness respectively. Here, the tangent operators are computed via a static condensation of the RVE structural stiffness matrix – see Appendix A.5.

CE

6. Numerical examples

AC

The numerical framework for the micromorphic computational homogenization approach is elaborated in Appendix A. The overall response of a composite material having a random distribution of inclusions, can be approximated reasonably well with the regular periodic stacking of RVEs (Van der Sluis et al., 2000). In this section, the predictive capabilities of the proposed formulation are illustrated with there benchmark examples under plane strain condition. In order to investigate the adequacy of the micromorphic framework in capturing the morphological influence, two RVE morpholo16

ACCEPTED MANUSCRIPT

CR IP T

gies with the same volume fraction of inclusions are considered, as shown in Fig. 4. To furthermore investigate the influence of RVE selection for a given morphology, four variations of RVE 1 depicted in Fig. 5 are considered.

RVE 2

RVE 1

M

AN US

Figure 4: Two RVE morphologies having 30% volume fraction of inclusions.

RVE 1a

RVE 1c

RVE 1d

ED

RVE 1b

Figure 5: Different RVEs considered for a given morphology.

PT

The von Mises plasticity model is adopted for describing the material behaviour of both the inclusion and matrix phases. The hardening rule is given by

CE

σy = σ0 + hp (εp )m ,

(29)

AC

with σy and σ0 denoting the current and initial yield strengths respectively, εp the effective plastic strain, hp and m the hardening modulus and exponent respectively. The problems considered are analysed with both conventional first-order and the micromorphic computational homogenization approaches, to illustrate the improved performance of the proposed formulation when benchmarked against the reference solutions from DNS. For the DNS, the mi-

17

ACCEPTED MANUSCRIPT

crostructural details are modelled explicitly with a mesh discretization of sufficient resolution.

AN US

CR IP T

6.1. Pure bending of a composite plate The first example considers the pure bending of a composite plate. Taking into account the antisymmetric deformation in pure bending, the analyses are restricted to the upper half of the plate. Furthermore, the structural axial response is periodic over a RVE size, due to the regular stacking of RVEs. In order to study the influence of different RVE choices, and also to eliminate any possible surface effect, we consider a plate segment comprising of 6 × 10 RVEs, as depicted in Fig. 6. For a given bending curvature (κ), the right edge of the plate segment is subjected to a linearly varying displacement profile (u1 = bκx2 ). An elastic-perfectly plastic von Mises model is adopted for both inclusion and matrix phases with the same elastic constants E = 200 GPa and ν = 0.3. The yield strength of the inclusion is assumed as 4 MPa, which is 100 times lower than that of the matrix.

M

b = 6l

AC

CE

x2

PT

2

RVE 1

ED

h = 10l

1

u1= b x2

or

x1

RVE 2

(a)

(b)

Figure 6: (a) A segment of the composite plate in pure bending. (b) Finite element discretization of the plate segment, and the underlying RVEs. The discretization at the macro-scale is such that a Gauss Point corresponds to each of the two locations indicated in (a), to facilitate a direct comparison with DNS solutions.

18

ACCEPTED MANUSCRIPT

1.4

1.4

DNS First-order Micromorphic

DNS First-order Micromorphic

1.2

Moment ( ×10 6 Nm)

AN US

1.2

Moment ( ×10 6 Nm)

CR IP T

As illustrated in Fig. 7, both the first-order and micromorphic computational homogenization approaches accurately predict the moment-curvature responses of the composite plate for different morphologies. Considering RVE 1, the axial strain profiles within a RVE at the two reference locations, indicated in Fig. 6a, are depicted in Fig. 8. It is easily observed that the predictions by the micromorphic framework compare well with the reference DNS solutions. In contrast, the bending mode contribution is lacking in the first-order approach. This illustrates the need for the higher order terms to adequately capture the higher-order deformation modes.

1 0.8 0.6 0.4

1

0.8 0.6 0.4 0.2

0.2

0

0 0.01

0.02

0.03

Curvature (m

-1

0.05

ED

(a)

0.04

)

M

0

0

0.01

0.02

0.03

Curvature (m

-1

0.04

0.05

)

(b)

AC

CE

PT

Figure 7: Moment-curvature responses obtained from the first-order and micromorphic approaches, benchmarked against the reference DNS solutions with (a) RVE 1 and (b) RVE 2.

19

ACCEPTED MANUSCRIPT

CR IP T

(a)

AN US

(b)

DNS

Micromorphic

1.4e-2

7.0e-3

First-order 0

ED

M

ˆ 11 within a RVE at locations (a) 1 and (b) 2 as indicated in Figure 8: Contour plots for H Fig. 6a, obtained from the micromorphic and the first-order homogenization framework, benchmarked against the reference DNS solution (deformation × 25).

AC

CE

PT

The same problem is next analysed with the different RVEs considered in Fig. 5 for the same microstructural morphology. The moment-curvature responses shown in Fig. 9 are almost identical for all RVEs considered, hence illustrating that the micromorphic solution is independent on the choice of RVE.

20

ACCEPTED MANUSCRIPT

1.4 RVE 1a RVE 1b RVE 1c RVE 1d

CR IP T

1

0.8

0.6

0.4

AN US

Moment ( ×10 6 Nm)

1.2

0.2

0 0

0.01

0.02

0.03

Curvature (m

-1

0.04

0.05

)

M

Figure 9: Moment-curvature responses obtained with the micromorphic approach for the different choices of RVE depicted in Fig. 5.

AC

CE

PT

ED

6.2. Shear deformation of a composite plate The applicability of the proposed micromorphic framework in situations without a strong separation of length scales is investigated in this section. We consider an infinitely long composite plate of height h subjected to a sinusoidal shear deformation, as illustrated in Fig. 10. In the following, a plate segment comprising of 10 RVEs is considered, with periodic boundary conditions imposed at the vertical surfaces. At its top and bottom surfaces, the displacement boundary conditions are prescribed as   2πx1 u1 = 0 , u2 = A sin , (30) λ where A and λ represent the loading amplitude and the wavelength of the prescribed deformation respectively. For the given deformation, a continuous model at the macro-scale will experience a shear distortion given by   2π 2πx1 H21 = ∇1 u2 = A cos . (31) λ λ 21

ACCEPTED MANUSCRIPT

CR IP T

In the following, we investigate the regime that violates the strong scale separation assumption inherent in a first-order homogenization approach. Specifically, we consider two load cases where the characteristic micro and macro length scales are of a similar order of magnitude, i.e. O(λ) ≈ O(l), where • Case 1: A = 0.01l, λ = 10l , • Case 2: A = 0.005l, λ = 5l .

Note that the amplitudes of the shear distortion (H21 ) are identical in the two load cases for the values of A and λ adopted.

PBC

AN US

λ

10 l

ED

M

x1

h

PT

η

CE

Figure 10: A composite plate segment comprising of 10 RVEs is subjected to a sinusoidal shear load, with periodic boundary conditions imposed at the vertical surfaces. Three different microstructural phases (η = 0, 3l , 2l 3 ) are considered.

AC

For the matrix, the material properties adopted are σ0 = 400 MPa, hp = 2 GPa, and m = 0.5, E = 200 GPa, and ν = 0.3. The inclusions are assumed to be elastic-perfectly plastic with E = 10 GPa, ν = 0.3, and σ0 = 20 MPa. A DNS of the problem involving full microstructural details is first conducted. In general, the actual microstructural arrangement of the composite plate is not known a priori. We thus consider three possible phase shifts with η = 0, 3l and 2l3 , as shown in Fig. 10. For all cases considered, 22

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

the average shear stress within each RVE is extracted from the DNS, which serve as the reference solution. The same problem is next solved with the conventional first-order, as well as the micromorphic computational homogenization framework. Finally, the macro shear stress profiles determined from the multi-scale approaches are compared against the reference DNS solutions. For load case 1 with a weak separation between the two characteristic length scales (λ = 10l), the first-order computational homogenization method provides a reasonable estimate with an error in the range of ∼ 5 % for the two RVEs considered, as depicted in Fig. 11. In both cases, a small, but observable improvement is obtained with the micromorphic framework. The results suggest that in the regime of a weak scale separation between micro and macro, the nonlocal interactions between the phases already have a non-negligible influence on the overall response. This nonlocal effect is further amplified in load case 2, where a scale separation no longer exists with λ = 5l. Accordingly, the first-order approach becomes inadequate, with a prediction error of ∼ 16 %, as shown in Fig. 12 for the two RVEs considered. In contrast, the micromorphic approach continues to predict an accurate structural response for the two morphologies, due to an adequate characterization of the underlying rapid fluctuations via constraints (5) and (6).

50

CE

-100

AC

0.2

0 -50 -100 -150

-150

0

50

12

σ

σ

-50

DNS First-order Micromorphic

100

PT

0

12

(MPa)

100

150

(MPa)

DNS First-order Micromorphic

ED

150

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

x /L

x /L

1

1

(b)

(a)

Figure 11: Macro shear stress (σ12 ) obtained from the first-order and micromorphic approaches, benchmarked against the reference DNS solutions, for load case 1 (λ = 10l) with (a) RVE 1 and (b) RVE 2.

23

ACCEPTED MANUSCRIPT

150

150

(MPa)

50

50 0

12

0

σ

σ

DNS First-order Micromorphic

100

12

(MPa)

100

-50

-50 -100

-100

-150

-150 0

0.2

0.4

0.6

0.8

0

1

0.2

CR IP T

DNS First-order Micromorphic

0.4

0.6

0.8

1

x /L

x /L

1

1

(b)

AN US

(a)

Figure 12: Macro shear stress (σ12 ) obtained from the first-order and micromorphic approaches, benchmarked against the reference DNS solutions, for load case 2 (λ = 5l) with (a) RVE 1 and (b) RVE 2.

AC

CE

PT

ED

M

To better illustrate the size effect, the macro shear stress determined from the first-order and micromorphic approaches are plotted against the normalized coordinate xλ1 for RVE 1, as shown in Fig 13. Recall that both load cases have identical strain amplitudes and that a standard von Mises model is adopted at the micro-scale. Therefore, any difference in the global responses is solely induced by the interactions between the micro-mechanisms characterized by the (macro) loading wavelength and the RVE size. The firstorder approach, assuming a separation of length scales, cannot distinguish between the two load cases, as evident from Fig 13(a). The micromorphic e ij , where an framework resolves this limitation with the enrichment of ∇k H associated length scale parameter characterizes the size of the underlying RVE. The extent of fluctuations within a RVE is thus scaled consistently, through this length scale parameter onto the macro-scale, which interacts with the prescribed loading to predict the correct size effect.

24

ACCEPTED MANUSCRIPT

150

150

Load case 1 Load case 2

0

12

0 -50 -100

-50 -100

-150

-150 0

0.2

0.4

0.6

0.8

1

0

x /λ

0.2

CR IP T

50

(MPa)

50

σ

σ

Load case 1 Load case 2

100

12

(MPa)

100

0.4

0.6

0.8

1

x /λ

1

1

(b)

AN US

(a)

Figure 13: Size effect induced by the interactions between loading wavelength (λ) and RVE size (l) is not captured by the (a) first-order approach, but accurately reproduced with the (b) micromorphic framework. Analyses based on RVE 1.

M

The problem is next analysed with RVE 1b in Fig. 5. Consistent with the bending example in Section 6.1, the macroscopic response for this problem is independent of the RVE definition for the same underlying microstructure, as shown in Fig 14.

ED

150

RVE 1a RVE 1b

100

0

CE

σ12 (MPa)

PT

50

-50

AC

-100 -150 0

0.2

0.4

0.6

0.8

1

x 1 /L

Figure 14: Global responses obtained from the micromorphic approach with different RVEs for load case 2.

25

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

6.3. Shear localization This section considers the localized deformation of a perforated plate induced by an embedded material imperfection. The plate is subjected to tensile loading in the x2 direction. Taking into account symmetry properties, only one quarter of the plate consisting of 5 × 10 RVEs is modelled, as shown in Fig. 15. The microstructure is represented by RVE 1 with a central void. The matrix material is assumed to be elastic-perfectly plastic, with E = 200 GPa, ν = 0.3, and σ0 = 400 MPa. Within the imperfection zone of size 2l, the yield stress of the matrix material is reduced by 20%. This imperfection region thus triggers the development of a localized shear band beyond a threshold loading. The problem is first analysed using DNS to establish the reference solution, and solved again with the first-order and micromorphic computational homogenization approaches.

Figure 15: A perforated plate with an imperfection zone subjected to tensile loading.

26

ACCEPTED MANUSCRIPT

12 DNS First-order Micromorphic

CR IP T

10

P (MN)

8

6

4

0 0

0.5

1

AN US

2

1.5

2

2.5

3

3.5

u 2 ( × 10 -4 m)

Figure 16: Load-displacement curves obtained form the DNS, first-order computational homogenization and micromorphic approach for the localization problem.

AC

CE

PT

ED

M

As illustrated in Fig. 16, the load-displacement response of the plate is well predicted by both the first-order and micromorphic computational homogenization approaches. Due to the presence of the imperfection region, however, a shear band initiates therewith and develops across the specimen. The shear bandwidth thus characterizes the localized deformation at the macro-scale, which is of the same order of magnitude as the size of a RVE. The development of a shear band thus violates the scale separation assumption in the first-order approach. To illustrate the inadequacy of the first-order approach, despite the match in global response against the reference solution, the volume average of axial strain (H22 ) within each RVE in the DNS is computed, and the corresponding contour plot depicted in Fig. 17. As observed, the axial strain profile obtained from the micromorphic approach match closely the reference profile. In contrast, the first-order approach predicts an inaccurate localization phenomenon with significantly higher strain values within the imperfection zone.

27

DNS

First-order

Micromorphic 1.25%

0.90%

0.75%

M

1.50%

AN US

CR IP T

ACCEPTED MANUSCRIPT

0.50%

0.25%

0

ED

Figure 17: Contour plots for H22 obtained form DNS, first-order computational homogenization and micromorphic approach at δ = 0.0035 m.

AC

CE

PT

Beyond a threshold deformation, the nonlinear geometrical effect will trigger a strain softening structural response. In such situations, the first-order approach loses its ellipticity in a quasi-static problem and suffers from the same mesh sensitivity issues that plague a standard continuum model. The second-gradient framework (Kouznetsova et al., 2004) was shown to resolve this issue by incorporating a RVE-related length scale parameter associated with the strain gradient term, which restores the well-posedness of the problem during strain softening. Physically, this length scale parameter characterizes the rapidly response within a RVE and propagates it to the macro process zone within the shear band. The overall structural behaviour is next obtained by accounting for the non-local interactions at the two-scales, thus naturally addressing the limitation of the first-order approach. Since the proposed micromorphic approach boasts a similar mechanism with the secondgradient framework as shown in (18), we expect it to resolve the same issue 28

ACCEPTED MANUSCRIPT

AN US

CR IP T

in presence of the geometrical strain softening, a glimpse of which is provided in Fig. 17 with a regularized shear band at infinitesimal deformation. Note, however, that both second-gradient and micromorphic framework are inadequate when the deformation is highly localized, e.g. a macroscopic crack that transverse across a RVE, such that the material responses fluctuate beyond a quadratic nature (Geers et al., 2010). For such cases, multi-scale schemes with discontinuity enrichments at the macro-scale should be adopted (e.g., Bosco et al., 2014, 2015; Coenen et al., 2012). To complete the analysis, the localization problem is next solved with RVE 1d in Fig. 5. Consistent with the other two examples, the choice of RVE does not significantly influence the global response, as shown in Fig 18.

12 RVE 1a RVE 1d

10

M

6

ED

P (MN)

8

4

PT

2

0

CE

0

0.5

1

1.5

2

u 2 ( × 10

-4

2.5

3

3.5

m)

AC

Figure 18: Global responses obtained from the micromorphic approach with different RVEs.

7. Conclusion This paper presents a novel computational homogenization approach that consistently translates standard continuum models at the micro-scale onto the macro-scale, where a generalized micromorphic continuum is recovered 29

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

without any a priori assumptions. The proposed framework introduces an additional macro kinematic variable to characterize the micro-mechanism underlying a rapidly fluctuating response. Although we have focused only on composites with soft inclusions in this paper, the proposed strategy can be adopted for other heterogeneous materials. Accordingly, the morphic field should describe the deformation in the matrix phase when considering a composite with hard inclusions. In another instance, an additional field characterizes the underlying displacement jump for the analytical homogenization of intergranular fracture (Sun and Poh, 2016). The net effect of resulting microstructural fluctuations is propagated to the macro-scale via the HillMandel condition, which in turn provides the definition of macro stress and the generalized stress tensors. The numerical framework of the micromorphic approach requires only C 0 continuity for the basic field variables. For the composites considered in this paper, the excellent predictive capability of the proposed homogenization model is demonstrated with three numerical examples. Considering the pure bending problem of a composite plate, it is shown that the micromorphic computational homogenization approach adequately predicts moment-curvature responses, as well as provides an accurate description of the second-order deformation within a RVE. The second example demonstrates the applicability of the proposed homogenization framework even in the absence of a clear separation between the loading wavelength and the RVE size. In this regime considered, the nonlocal interactions between different underlying mechanisms induce an overall size effect at the macro-scale. This is captured naturally by the proposed micromorphic approach via the internal length scale parameter associated with the higher-order morphic field. Finally, the shear localization problem illustrates a regularizing effect of the micromorphic approach. In the presence of large spatial gradients spanning across only a few RVEs, an accurate shear band is predicted, in contrast to an erroneous localized deformation obtained with the first-order approach. For all the examples considered, it is shown that the micromorphic solutions are independent of the RVE choice, for a given microstructure. Appendix A. Numerical implementation In the following, essential details of the numerical framework are presented.

30

ACCEPTED MANUSCRIPT

AN US

CR IP T

Appendix A.1. Nested solution scheme (F E 2 ) The proposed homogenization framework is implemented numerically by adopting a nested solution scheme, where BVPs at the macro and microscales are solved simultaneously using the finite element method. The steps involved are illustrated in Fig. A.19 and summarized as below. The macroscopic domain is discretized with finite elements and the defore ij and ∇k H e ij , is determation at each Gauss point, characterized by Hij , H e mined from the discretized primary variables ui and Hij . This macro-scale information is next translated to the micro-scale to enable the imposition of kinematic constraints on the RVEs, as described in Section 5.1. Upon resolving for the resulting micro-scale BVP, macroscopic stress tensors (σij , ξij , ζijk ) are extracted from the RVE and propagated back to the macro Gauss Points. Additionally, the tangent operators are computed through a static condensation of RVE stiffness matrix and passed onto the macro-scale. The macro-scale BVP is solved with an implicit iterative framework. At each iteration, the element stiffness matrices (ele K) and internal element

and

ele

M

nodal force vectors (ele Fint ) are computed based on the stress tensors and tangent operators extracted from the RVE analyses. The assembly of ele K Fint leads to the global stiffness matrix (K) and the vector of internal

AC

CE

PT

ED

forces (Fint ) respectively.

31

ACCEPTED MANUSCRIPT

Macro

Micro

Macro Element: i=1,2, ...n Gauss Point: j=1,2, ...m

CR IP T

Start RVE analysis

Apply kinematic constraints

No j=j+1 Store stress tensors, tangent stiffness matrices

tangents

No i=i+1

Yes Compute eleK and eleFint

Yes

M

i = n

AN US

j = m

Extract tangent siffness matrices and stress tensors

ED

Assembly of K and F int

Check for convergence

PT

Figure A.19: Nested solution scheme

AC

CE

Appendix A.2. Weak formulation of macro-scale BVP The governing equations for the micromorphic continuum (15a) and (15b) are satisfied simultaneously in the weak form with the respective kinematic variables as the weight functions Z δui σij,j dV = 0 , (A.1) Ω Z e ij (ξij − ζijk,k ) dV = 0 . δH (A.2) Ω

32

ACCEPTED MANUSCRIPT

Ω

Ω

Γ

(A.3) (A.4)

CR IP T

Applying the divergence theorem, Z Z δHij σij dV = δui ti dS , Ω Z Γ Z Z e e e ij mij dS . δ Hij ξij dV + δ Hij,k ζijk dV = δ H

AN US

The stress quantities are next linearized such that Z Z Z δHij dσij dV = − δHij σij dV + δui ti dS , (A.5) Ω Ω Γ Z Z Z Z e e e e ij,k ζijk dV δ Hij dξij dV + δ Hij,k dζijk dV = − δ Hij ξij dV − δ H Ω Ω Ω Ω Z e ij mij dS . + δH (A.6) Γ

M

Appendix A.3. Spatial discretization To facilitate finite element implementation, the stress tensors and kinematic variables are written in Voigt notations, where a single and a double underscore represent a vector and a matrix respectively. The basic variables are discretized as e = Ne ae , H H H

,

(A.7)

ED

u = Nu au

where Nu and NHe are quadratic and linear shape functions, au and aHe the

PT

nodal degrees of freedom. The gradient terms are similarly discretized as ,

CE

H = Bu au

e = Be ae , ∇H H H

(A.8)

with Bu and BHe as the corresponding gradient operator matrices.

AC

Referring to (28), the incremental stress quantities are discretized as

e + C(1) d∇H e , dσ = C (1) dH + C (2) dH

e + C(2) d∇H e , dξ = C (3) dH + C (4) dH e + C d∇H e . dζ = C(3) dH + C(4) dH

(A.9) (A.10) (A.11)

33

ACCEPTED MANUSCRIPT

CR IP T

Utilizing the spatial discretizations in (A.7)–(A.11), the linearized system of equations in (A.5) and (A.6) are assembled to form the numerical framework as   K11 K12  da   F  u 1   = , (A.12) daHe K21 K22 F2 with

Ω

K21 =

K22 =

Z

BuT C (2) NHe dV

Ω

NHTe C (3) Bu dV

Ω

NHTe C (4) NHe dV

Z

+

Z

Ω

F1 =

Z

F2 =

Z

Γ

+

NuT t

dS −

NHTe m

Z

Ω

Ω

+

Z

Ω

T BH e dV , e C BH

Z

BuT C(1) BHe dV ,

T (3) BH e C Bu dV ,

NHTe C(2) BHe dV

dS −

+

Z

BuT σdV ,

Ω

Z

PT

Γ

+

Z

AN US

K12 =

Z

(A.13)

M

Ω

BuT C (1) Bu dV ,

ED

K11 =

Z

Ω

NHTe ξdV

Ω

T (4) BH e dV e C NH

(A.14) (A.15)

(A.16)

(A.17) −

Z

Ω

T BH e ζdV .

(A.18)

AC

CE

Appendix A.4. Macroscopic stress tensors For ease of numerical implementation, the macroscopic stress tensors defined in (12) are written in terms of equivalent surface integrals. To this end, we first define the total potential energy of a RVE as 1 Π = 2 ∗

Z

Vˆ

ˆ ij dVˆ −λij σ ˆij H



1 Vˆ

Z

∂ Vˆ

  Z  1 ˜ ˆ ˆ e uˆi n ˆ j dS − Hij −λij uˆi n ˆ j dS − Hij , VˆI ∂ Iˆ (A.19)

34

ACCEPTED MANUSCRIPT

˜ ij are the Lagrange multipliers associated with constraints where λij and λ (3) and (4) respectively. Satisfying δΠ∗ = 0 , the Lagrange multipliers are identified as (A.20)

CR IP T

on ∂ Vˆ , on ∂ Iˆ ,

λij n ˆj = σ ˆij n ˆj I I ˜ ij n λ ˆ j = (ˆ σij − σ ˆijM )ˆ nIj

(A.21)

AN US

with superscripts M and I denoting the matrix and inclusion phases respectively. Note that condition (3) is satisfied via the periodic boundary conditions described in Section 5.1. The tying forces arising from constraints (24a) and (24b) thus correspond to the Lagrange multiplier λij in (A.19). Applying the divergence theorem, (10) is next rewritten as Z Z 1 1 ˙ ˜ ij n ˆ hP (x, y)i = σ ˆij n ˆ j uˆi dS + λ ˆ Ij uˆ˙ i dS , (A.22) Vˆ ∂ Vˆ Vˆ ∂ Iˆ

ED

M

where use has been made of (A.21). Substituting (1) into (A.22) we obtain  Z Z 1 I ˜ ˆ ˆ ˆ e˙ ij hP (x, y)i = σ ˆik n ˆ k y j dS + λik n ˆ k y j dS H ˆ ∂ Vˆ ∂ Iˆ  Z V Z 1 I I ˜ ˜ ˆ ˆ e˙ ij,k + (ˆ σiq n ˆ q yj yk + σ ˆjq n ˆ q yi yk )dS + (λiq n ˆ q yj yk + λjq n ˆ q yi yk )dS H 4VˆZ ∂ Vˆ ∂ Iˆ Z 1 1 ˜ ij n σ ˆij n ˆ j uˆ˙ ∗i dSˆ + λ ˆ Ij uˆ˙ ∗i dSˆ . (A.23) + ˆ ˆ V ∂ Vˆ V ∂ Iˆ

PT

Considering (A.20) and the constraints in (5) and (6), (A.23) can be further simplified into Z

∂ Vˆ

Z VˆI ˜ e˙ 1 ˆ ˙ e˙ ij,k . ˆ ti yj dS Hij + λij H ij + (tˆi yj yk + tˆj yi yk )dSˆ H ˆ ˆ V 4V ∂ Vˆ (A.24)

AC

CE

1 hPˆ (x, y)i = ˆ V

35

ACCEPTED MANUSCRIPT

(A.25a)

CR IP T

Comparing (A.24) with (11b), the macro stresses are obtained as Z 1 σij = tˆi yj dSˆ , Vˆ ∂ Vˆ VˆI ˜ λij , ξij = Vˆ Z
1 tˆi yj yk + tˆj yi yk dSˆ . ζijk = 4Vˆ ∂ Vˆ

(A.25b)

(A.25c)

σij =

AN US

The surface integrals in (A.25) are evaluated numerically by utilizing the reaction forces at the boundary nodes. For the periodic boundary conditions in (24), only the corner nodes contribute towards (A.25). Accordingly, the macroscopic stress tensors at a macro point are obtained as 4 1 X (c) (c) fi yj , Vˆ c=1

(A.26b) (A.26c)

ED

M

VˆI ˜ λij , ξij = Vˆ 4  1 X  (c) (c) (c) (c) (c) (c) ζijk = fi y j y k + fj y i y k , 4Vˆ c=1

(A.26a)

where the superscript c denotes the location of corner nodes depicted in Fig. 3.

AC

CE

PT

Appendix A.5. Macroscopic tangent operators The linearized system of equations at the micro-scale for each RVE can be written as      ˆ id  ˆ ii K K ui dˆ ri ˆ dˆ  dˆ  u = dˆ r, (A.27) = i.e. K ˆ di K ˆ dd dˆ ud dˆ rd K

where subscript d denotes the dependent Degrees of Freedom (DOF) involved in the periodic boundary conditions (24a) and (24b), which are related to the independent DOFs uˆi in the system as dˆ ud = Cˆ dˆ ui , with Cˆ the associated coefficient matrix. The vector dˆ r represents the residual forces. Elimination of the dependent DOFs leads to the following condensed system of equations 36

ACCEPTED MANUSCRIPT

(Kouznetsova et al., 2004) ˆ ∗ dˆ K ui = dˆ r∗ , where ˆ∗ = K ˆ ii + K ˆ id Cˆ + Cˆ T K ˆ di + Cˆ T K ˆ dd Cˆ , K T

rd . dˆ r∗ = dˆ ri + Cˆ dˆ

CR IP T

(A.28)

(A.29a)

(A.29b)

AN US

˜ are next introduced in order to Dummy DOFs (denoted by subscript λ) impose the constraint (26) such that e = uˆ ˜ , W uˆI = H λ

(A.30)

ED

M

with W and uˆI representing the coefficient matrix in constraint (26) and DOFs along the matrix-inclusion interface respectively. The components of e are prescribed directly on the dummy nodes. Rearranging the terms, H (A.30) is written as h i uˆ  a ˆ ˆ uˆId = C1 C2 , (A.31) uˆλ˜

AC

CE

PT

with uˆId denoting the dependent DOFs on the inclusion-matrix interface, and uˆa the rest of the retained DOFs in (A.28) respectively. ˆ ∗ ) is further partitioned and augThe condensed RVE stiffness matrix (K mented to accommodate the dummy nodes as  ∗  ˆ ˆ ∗∗   K 0 K   aa aId dˆ ra∗ ua   dˆ  0 rλ∗˜  I 0  uλ˜  =  (A.32)    dˆ  dˆ  ,  ∗  ∗ ∗ dˆ uId ˆ ˆ dˆ rId K Id a 0 KId Id where I is the identity matrix. Subsequently, the dependent DOFs (ˆ uId ) are

removed from (A.32) by making use of the condensation procedure adopted

37

ACCEPTED MANUSCRIPT

where

ˆ ∗∗ = K ˆ∗ + K ˆ ∗ Cˆ1 + Cˆ1 T K ˆ ∗ + Cˆ1 T K ˆ ∗ Cˆ1 , K aa aa aId Id a Id Id ˆ ∗ Cˆ2 + Cˆ1 T K ˆ ∗ Cˆ2 , ˆ ∗∗˜ = K K aId Id Id aλ

ˆ ˜∗∗˜ = I + Cˆ2 T K ˆ ∗ Cˆ2 , K Id Id λλ T dˆ ra∗∗ = dˆ ra∗ + Cˆ1 dˆ rI∗d ,

M

T

rλ∗˜ + Cˆ2 dˆ dˆ rλ∗∗ rI∗d . ˜ = dˆ

PT

ED

The system (A.33) is next partitioned as   ∗∗ ˆ pp ˆ ∗∗˜ K ˆ ∗∗ K K pf   ∗∗    pλ   u  dˆ rp df p  ˆ ∗∗ ˆ ∗∗ ˆ ∗∗  dˆ K˜ K˜ ˜ K˜  dˆ  ∗∗  ˜ ,   dˆ r u ≈ = λp λ λ λf ˜ dλ    λ˜  λ   ∗∗ 0 uf  ˆ ∗∗ ˆ ∗∗ ˆ ∗∗  dˆ dˆ rf K K ˜ K fλ

(A.34c) (A.34d) (A.34e) (A.34f)

(A.35)

ff

CE

fp

(A.34a)

(A.34b)

AN US

ˆ ˜∗∗ = Cˆ2 T K ˆ ∗ + Cˆ2 T K ˆ ∗ Cˆ1 , K Id a Id Id λa

(A.33)

CR IP T

in (A.28). The reduced system of equations is expressed as   ∗∗ ˆ ∗∗˜  ˆ aa  " ∗∗ # K K aλ dˆ ra dˆ u   a = ,  ˆ ∗∗ ˆ ∗∗  dˆ dˆ rλ∗∗ uλ˜ Kλa Kλ˜λ˜ ˜ ˜

AC

with subscript p and f denoting the corner nodes where displacement values are prescribed and internal free nodes respectively. By performing static condensation on (A.35), the system of equations is reduced to the final form (Kouznetsova et al., 2004) ¯ ¯ ˜     Kpp K pλ dˆ up df   = ¯˜ K ¯ ˜˜ ˜ . K dˆ uλ˜ dλ λp λλ

(A.36)

38

ACCEPTED MANUSCRIPT

In view of (23), the vector dˆ up can be extracted from the linearization prescribed displacements at the corner nodes for c = 1, 2, 3, 4 (A.37)

CR IP T

(c) (c) 1 e ij +d∇k H e ji ) y (c) y (c) , dˆ ui = dHij yj + (d∇k H j k 4

e Drawing reference to (A.26), the whereas dˆ uλ˜ is directly related to dH. variation in the macro stress tensors can be identified as dσij =

4 1 X (c) (c) dfi yj Vˆ

(A.38a)

c=1

AN US

VˆI ˜ dξij = dλij Vˆ 4  1 X  (c) (c) (c) (c) (c) (c) dζijk = dfi yj yk + dfj yi yk 4Vˆ c=1

(A.38b)

(A.38c)

AC

CE

PT

ED

M

The macro tangent operators, defined in (28), can be computed next utilizing

39

ACCEPTED MANUSCRIPT

(A.36) and (A.38) as 4 4 1 X X ¯ (ab) = Kpp Vˆ

(3)

Cαjk = (4)

Cαβ

4 ˜ VˆI X ¯ (λa) Kλp ˜ Vˆ a=1

(A.39b)

, y (a) (k)

(A.39c)

(α,j)

VˆI ¯ Kλ˜λ˜ = , Vˆ (α,β)

(1)

Cijklm

(A.39a)

(i,k) a=1 b=1 4 X (aλ) ¯ ˜ y (a) , K ˜ (j) pλ a=1 (i,β)

CR IP T

1 Vˆ

(2)

Cijβ =

y (a) y (b) , (j) (l)

" 4 4 1 X X ¯ (ab) = Kpp 4Vˆ a=1 b=1



(ab) ¯ pp y (a) y (b) y (b) + K (j (l) (m) (i,k)

(A.39d)

AN US

(1) Cijkl

y (a) y (b) y (b) (j) (k) (m)

(i,l)

#

, (A.39e)



PT

ED

M

4 ˜ ˜ VˆI X  ¯ (λa) (a) ¯ (λa) Kλp y (a) y + K y (a) y (a)  , (A.39f) ˜ ˜ (k) (l) (j) (l) λp ˆ 4V a=1 (α,j) (α,k) " # 4 4 X X 1 (3) ¯ (ab) y (a) y (a) y (b) + K ¯ (ab) y (a) y (a) y (b) , K Cijklm = pp pp (j) (k) (m) (i) (k) (m) ˆ 4V a=1 b=1 (i,l) (j,l) (A.39g)   4 ˜ ˜ 1 X  ¯ (aλ) (4) ¯ (aλ) Cijkβ = Kpλ˜ y (a) y (a) + K y (a) y (a)  , (A.39h) ˜ (j) (k) (i) (k) pλ ˆ 4V (2)

Cαjkl =

a=1

1 16Vˆ

CE

(1)

AC

Cijklmn =

(i,β)

4 X 4  X a=1 b=1

¯ ab +K pp

(j,β)

(ab) ¯ pp y (a) y (a) y (b) y (b) + K (j) (k) (m) (n)

(ab) ¯ pp K

(i,l)

y (a) y (a) y (b) y (b) (i) (k) (m) (n)

y (a) y (a) y (b) y (b) (j) (k) (l) (n) (i,m)

¯ (ab) +K pp

(j,l)

y (a) y (a) y (b) y (b) (i) (k) (l) (n) (j,m)



,

(A.39i)

(ab) ¯ pp ¯ pp at the location of rows and columns where K refers to the submatrix of K ˜ ¯ (aλ) corresponding to the corner nodes a and b respectively. Similarly, K and ˜ pλ

40

ACCEPTED MANUSCRIPT

˜ ¯ ˜ and K ¯ ˜ associated with the corner node ¯ (λa) denote the submatrices of K K ˜ λp pλ λp

CR IP T

a and dummy DOFs. The matrix and vector components are indexed using the subscripts i, j, k, l, m, n = 1, 2 and α, β = 1, 2, 3, 4. Note that the Voigt notation (subscript α and β ) is used for the stress tensor ξij and the morphic e ij to facilitate numerical computation. field H Acknowledgements

AN US

The work presented in this paper is supported by the MOE Tier 2 Grant R302000139112. References

Bacigalupo, A., Gambarotta, L., 2010. Second-order computational homogenization of heterogeneous materials with periodic microstructure. ZAMMJournal of Applied Mathematics and Mechanics/Zeitschrift f¨ ur Angewandte Mathematik und Mechanik 90 (10-11), 796–811.

ED

M

Bosco, E., Kouznetsova, V., Coenen, E., Geers, M., Salvadori, A., 2014. A multiscale framework for localizing microstructures towards the onset of macroscopic discontinuity. Computational Mechanics 54 (2), 299–319.

PT

Bosco, E., Kouznetsova, V., Geers, M., 2015. Multi-scale computational homogenization–localization for propagating discontinuities using x-fem. International Journal for Numerical Methods in Engineering 102 (3-4), 496–527.

CE

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

AC

Eringen, A. C., 1968. Mechanics of micromorphic continua. In: Mechanics of generalized continua. Springer, pp. 18–35. Feyel, F., 2003. A multilevel finite element method (fe 2) to describe the response of highly non-linear structures using generalized continua. Computer Methods in applied Mechanics and engineering 192 (28), 3233–3244.

41

ACCEPTED MANUSCRIPT

CR IP T

Feyel, F., Chaboche, J.-L., 2000. Fe 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials. Computer methods in applied mechanics and engineering 183 (3), 309– 330. Forest, S., 2009. Micromorphic approach for gradient elasticity, viscoplasticity, and damage. Journal of Engineering Mechanics 135 (3), 117–131. Forest, S., 2016. Nonlinear regularization operators as derived from the micromorphic approach to gradient elasticity, viscoplasticity and damage. In: Proc. R. Soc. A. Vol. 472. The Royal Society, p. 20150755.

AN US

Forest, S., Sab, K., 1998. Cosserat overall modeling of heterogeneous materials. Mechanics Research Communications 25 (4), 449–454. Forest, S., Sievert, R., 2003. Elastoviscoplastic constitutive frameworks for generalized continua. Acta Mechanica 160 (1-2), 71–111. Forest, S., Sievert, R., 2006. Nonlinear microstrain theories. International Journal of Solids and Structures 43 (24), 7224–7245.

ED

M

Forest, S., Trinh, D. K., 2011. Generalized continua and non-homogeneous boundary conditions in homogenisation methods. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift f¨ ur Angewandte Mathematik und Mechanik 91 (2), 90–109.

PT

Geers, M. G., Kouznetsova, V. G., Brekelmans, W., 2010. Multi-scale computational homogenization: Trends and challenges. Journal of computational and applied mathematics 234 (7), 2175–2182.

CE

Germain, P., 1973. The method of virtual power in continuum mechanics. part 2: Microstructure. SIAM Journal on Applied Mathematics 25 (3), 556–575.

AC

Ghosh, S., Lee, K., Moorthy, S., 1995. Multiple scale analysis of heterogeneous elastic structures using homogenization theory and voronoi cell finite element method. International Journal of Solids and Structures 32 (1), 27– 62. Hirschberger, C. B., Sukumar, N., Steinmann, P., 2008. Computational homogenization of material layers with micromorphic mesostructure. Philosophical Magazine 88 (30-32), 3603–3631. 42

ACCEPTED MANUSCRIPT

H¨ utter, G., 2016. Homogenization of a cauchy continuum towards a micromorphic continuum. Journal of the Mechanics and Physics of Solids.

CR IP T

J¨anicke, R., Diebels, S., 2010. Numerical homogenisation of micromorphic media. Technische Mechanik 30 (4), 364–373.

J¨anicke, R., Diebels, S., Sehlhorst, H.-G., D¨ uster, A., 2009. Two-scale modelling of micromorphic continua. Continuum Mechanics and Thermodynamics 21 (4), 297–315.

AN US

J¨anicke, R., Steeb, H., 2012. Minimal loading conditions for higher-order numerical homogenisation schemes. Archive of Applied Mechanics 82 (8), 1075–1088.

Kaczmarczyk, L., Pearce, C. J., Bi´cani´c, N., 2008. Scale transition and enforcement of rve boundary conditions in second-order computational homogenization. International Journal for Numerical Methods in Engineering 74 (3), 506–522.

M

Kaczmarczyk, L., Pearce, C. J., Bi´cani´c, N., 2010. Studies of microstructural size effect and higher-order deformation in second-order computational homogenization. Computers & structures 88 (23), 1383–1390.

ED

Klusemann, B., Svendsen, B., 2010. Homogenization methods for multi-phase elastic composites: Comparisons and benchmarks. Technische Mechanik 30 (4), 374–386.

CE

PT

Kouznetsova, V., Brekelmans, W., Baaijens, F., 2001. An approach to micromacro modeling of heterogeneous materials. Computational Mechanics 27 (1), 37–48.

AC

Kouznetsova, V., Geers, M., Brekelmans, W., 2004. Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy. Computer Methods in Applied Mechanics and Engineering 193 (48), 5525–5550. Kouznetsova, V., Geers, M. G., Brekelmans, W. M., 2002. Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme. International Journal for Numerical Methods in Engineering 54 (8), 1235–1260. 43

ACCEPTED MANUSCRIPT

Larsson, F., Runesson, K., 2007. Rve computations with error control and adaptivity: the power of duality. Computational Mechanics 39 (5), 647– 661.

CR IP T

Larsson, R., Diebels, S., 2007. A second-order homogenization procedure for multi-scale analysis based on micropolar kinematics. International journal for numerical methods in engineering 69 (12), 2485–2512.

AN US

Li, X., Liang, Y., Duan, Q., Schrefler, B., Du, Y., 2014. A mixed finite element procedure of gradient cosserat continuum for second-order computational homogenisation of granular materials. Computational Mechanics 54 (5), 1331–1356.

Matsui, K., Terada, K., Yuge, K., 2004. Two-scale finite element analysis of heterogeneous solids with periodic microstructures. Computers & Structures 82 (7), 593–606.

M

Miehe, C., Dettmar, J., Z¨ah, D., 2010. Homogenization and two-scale simulations of granular materials for different microstructural constraints. International Journal for Numerical Methods in Engineering 83 (8-9), 1206– 1236.

ED

Miehe, C., Koch, A., 2002. Computational micro-to-macro transitions of discretized microstructures undergoing small strains. Archive of Applied Mechanics 72 (4-5), 300–317.

CE

PT

Nguyen, V.-D., Becker, G., Noels, L., 2013. Multiscale computational homogenization methods with a gradient enhanced scheme based on the discontinuous galerkin formulation. Computer Methods in Applied Mechanics and Engineering 260, 63–77.

AC

Nguyen, V.-D., Noels, L., 2014. Computational homogenization of cellular materials. International Journal of Solids and Structures 51 (11), 2183– 2203. Pierard, O., Friebel, C., Doghri, I., 2004. Mean-field homogenization of multiphase thermo-elastic composites: a general framework and its validation. Composites Science and Technology 64 (10), 1587–1603.

44

ACCEPTED MANUSCRIPT

Poh, L., 2013. Scale transition of a higher order plasticity model–a consistent homogenization theory from meso to macro. Journal of the Mechanics and Physics of Solids 61 (12), 2692–2710.

CR IP T

Smit, R., Brekelmans, W., Meijer, H., 1998. Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi-level finite element modeling. Computer Methods in Applied Mechanics and Engineering 155 (1), 181–192.

AN US

Sun, G., Poh, L., 2016. Homogenization of intergranular fracture towards a transient gradient damage model. Journal of the Mechanics and Physics of Solids 95, 374–392.

Suquet, P., 1985. Local and global aspects in the mathematical theory of plasticity. Plasticity today: modelling, methods and applications, 279–310. Terada, K., Kikuchi, N., 2001. A class of general algorithms for multi-scale analyses of heterogeneous media. Computer methods in applied mechanics and engineering 190 (40), 5427–5464.

AC

CE

PT

ED

M

Van der Sluis, O., Schreurs, P., Brekelmans, W., Meijer, H., 2000. Overall behaviour of heterogeneous elastoviscoplastic materials: effect of microstructural modelling. Mechanics of Materials 32 (8), 449–462.

45