Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy

Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550 www.elsevier.com/locate/cma

Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy V.G. Kouznetsova b

a,b,*

, M.G.D. Geers b, W.A.M. Brekelmans

b

a Netherlands Institute for Metals Research, Rotterdamseweg 137, 2628 AL Delft, The Netherlands Department of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands

Received 7 May 2003; received in revised form 17 December 2003; accepted 19 December 2003

Abstract This paper presents the detailed implementation and computational aspects of a novel second-order computational homogenization procedure, which is suitable for a multi-scale modelling of macroscopic localization and size effects. The second-order scheme is an extension of the classical (first-order) computational homogenization framework and is based on a proper incorporation of the gradient of the macroscopic deformation gradient tensor into the kinematical macro–micro scale transition. From the microstructural analysis the macroscopic stress and higher-order stress tensors are obtained, thus delivering a microstructurally based constitutive response of the macroscopic second gradient continuum. The higher-order macroscopic constitutive tangents are derived through static condensation of the microscopic global tangent matrix. For the solution of the second gradient equilibrium problem on the macrolevel a mixed finite element formulation is developed. As an example, the second-order computational homogenization approach is applied for the multi-scale analysis of simple shear of a constrained heterogeneous strip, where a pronounced boundary size effect appears.
2004 Elsevier B.V. All rights reserved. Keywords: Heterogeneous materials; Multi-scale modelling; Coarse graining; Computational homogenization; Higher-order constitutive modelling; Second gradient continuum; Size effects

*

Corresponding author. Address: Department of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands. Tel.: +31 40 247 22 45; fax: +31 40 244 73 55. E-mail address: [email protected] (V.G. Kouznetsova). 0045-7825/$ - see front matter
2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.12.073

5526

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

1. Introduction Homogenization techniques are known to be excellent tools to predict the effective linearly elastic properties of heterogeneous materials. Homogenization of solids in a geometrically and physically non-linear regime is clearly more difficult [1,2]. Moreover, most of the existing homogenization methods are not suitable to deal with large deformations and complex loading paths and cannot account for an evolving microstructure. To overcome these problems a computational homogenization approach has been developed, which is essentially based on the solution of two (nested) boundary value problems (BVP), one for the macroscopic and one for the microscopic scale. Techniques of this type (i) do not require any constitutive assumption with respect to the overall material behaviour; (ii) enable the incorporation of large deformations and rotations on both micro and macrolevel; (iii) are suitable for arbitrary material behaviour, including physically non-linear and time dependent behaviour; (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 modelling technique at the microlevel for the solution of the corresponding BVP. First-order computational homogenization schemes [3–14] fit entirely into a standard local continuum mechanics framework. The macroscopic deformation (gradient) tensor is calculated for every material point of the macrostructure and is next used to formulate kinematic boundary conditions to be applied on the associated microstructural representative volume element (RVE). Only the first gradient of the macroscopic displacement field is included in the framework, which thus is called ‘‘first-order’’. After the solution of the microstructural boundary value problem, the macroscopic stress tensor is obtained by averaging the resulting microstructural stress field over the volume of the microstructural cell. As a result, the (numerical) stress–strain relationship at every macroscopic point is readily available. The first-order computational homogenization technique proves to be a versatile strategy to retrieve the macroscopic mechanical response of non-linear multi-phase materials. However, there are some severe limitations in the application of a first-order computational homogenization scheme (as also the case for other conventional homogenization methods), which in fact originate from the underlying fundamental assumption that the microstructural length scale is infinitely small compared to the characteristic macrostructural size. In a growing number of applications, however, the microstructural length scale becomes comparable with the characteristic size of the macrocomponent or with the length scale associated with the macroscopic spatial variability in the loading. This may result in different types of size effects and/or zones of highly localized deformation. These phenomena cannot be modelled using a classical (local) continuum approach and consequently neither by using a first-order homogenization method, which is based on the principles of a local continuum. To be able to deal with this class of phenomena on the macroscopic scale different authors proposed to use higher-order or non-local continua, e.g. [15–22] amongst others. The use of an enhanced continuum has the effect of introducing a material length scale into the constitutive relation providing the necessary regularization and allowing to capture size effects. These continuum macroscopic models are either purely phenomenological or only implicitly incorporate the presence of the underlying microstructure. Explicitly considering the microstructure, in the framework of closed-form homogenization, higher-order continua have been used to describe the behaviour of either the microstructural constituents or the homogenized macrostructure (or both) [23–31]. However, the above works involving higher-order homogenization were limited to geometrically linear kinematics and were mostly concerned with retrieval of homogenized material parameters for certain cases rather than with providing a general computational procedure. Recent extensions of the computational homogenization methodology to higher-order continua, resulting in a so-called second-order computational homogenization scheme [32–34], provide a systematic procedure to obtain the constitutive response of a macroscopic second gradient continuum, based on the behaviour of the underlying microstructure. The second-order scheme is based on a proper incorporation

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5527

of the gradient of the macroscopic deformation gradient tensor into the kinematical macro–micro transition. The macroscopic stress tensor and a higher-order stress tensor are retrieved using an extended version of the Hill–Mandel energy balance. A full second gradient continuum appears at the macroscale, while the microstructural boundary value problem remains classical. Thus, the second-order computational homogenization scheme directly provides a microstructurally based higher-order constitutive response without additional assumptions on the macrolevel. The general framework of the second-order computational homogenization approach has been presented in previous publications [32,33]. This paper concentrates on the detailed implementation and related computational issues of the fully coupled micro–macro second-order computational homogenization scheme. Cartesian tensors and tensor products are used throughout the paper: ~ a, A and nA denote respectively a vector, a second-order tensor and an nth-order tensor. The following notation for vector and tensor operations is employed: conjugation Acij ¼ Aji , the dyadic product ~ a~ b ¼ ai bj~ ei~ ej and the scalar products A
B ¼ Aij Bjk~ ei~ ek ;

A : B ¼ Aij Bji ;

3

. A .. 3 B ¼ Aijk Bkji ;

~ a
3A
~ b ¼ ai Aijk bk~ ej ; with ~ ei , i = 1, 2, 3 the unit vectors of a Cartesian basis. A matrix and a column are denoted by A and a,  respectively.

2. Second-order computational homogenization 2.1. General framework ~) to the For a continuum the transformation from the undeformed macroscopic state (position vector X deformed state (~ x) is considered. For an infinitesimal material line element, classical continuum mechanics leads to the linear mapping ~; d~ x ¼ FM
dX

ð1Þ

c with the deformation gradient tensor defined by FM ¼ ðr0M~ xÞ ¼ r0Mi xj~ ej~ ei , where the gradient operator is taken with respect to the macroscopic reference configuration. Here and in the following the subscript ‘‘M’’ refers to a macroscopic quantity, while the subscript ‘‘m’’ will denote a microscopic quantity. The linear relation (1) is used as the point of departure in the classical (first-order) computational homogenization. When dealing with line elements in non-homogeneously deformed volumes of a finite size, relation (1) does not apply any more, and an expression for a finite material vector D~ x in the current macroscopic configuration may be obtained using a Taylor series expansion 3

~ þ 1DX ~
3 GM
DX ~ þ OðDX ~ Þ; D~ x ¼ FM
DX 2

ð2Þ

where the third-order tensor 3GM is introduced as 3

GM ¼ r0M FM ¼ r0Mi F jk~ ei~ ej~ ek :

ð3Þ

The second-order computational homogenization approach uses both the macroscopic deformation tensor FM and its gradient 3GM = $0MFM to prescribe kinematic boundary conditions on an RVE. After the solution of the RVE boundary value problem (which will be summarized in Section 2.4), a macroscopic stress tensor PM and a higher-order stress tensor 3QM (a third-order tensor), defined as the work conjugate quantity of the gradient 3GM of the deformation gradient tensor, are derived exploiting the Hill–Mandel energy

5528

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

Fig. 1. Second-order computational homogenization scheme.

condition. Further details are given in Section 2.3. A scheme of the second-order computational homogenization is presented in Fig. 1. Contrary to the microstructural problem, which remains a classical one, on the macrolevel a full gradient second-order equilibrium problem appears. Its formulation and its finite element implementation are discussed in Section 5. In contrast to most higher-order theories, the framework presented here does not bring in complicating factors due to the direct numerical quantitative evaluation of the higher-order constitutive behaviour. Macroscopic quantities are obtained from the microstructural analysis according to a straightforward multi-scale mathematical derivation. A procedure for the extraction of macroscopic higher-order tangent operators is presented in Section 4.2. Finally, the whole framework is summarized in Section 6 where the coupled micro–macro solution scheme is outlined. 2.2. Micro–macro kinematics The boundary value problem specification for a microstructural RVE departs from the Taylor series expansion (2), where all higher-order terms are condensed into an unknown microfluctuation field D~ w. This microfluctuation field accounts for the effect of a local microscopic (fine scale) displacement field, which is superimposed on the local macroscopic (coarse scale) displacement field. Hence, ~ þ 1DX ~
3 GM
DX ~ þ D~ D~ x ¼ FM
DX w: 2

ð4Þ

This relation is applied to an undeformed volume of the microstructural representative cell with respect to ~c , located in ~ its geometrical center X xc after deformation. For simplicity the origin of the used Cartesian ~c ¼ ~ ~¼X ~X ~c ¼ X ~ vector basis is placed at geometric center of the undeformed RVE, i.e. X 0, thus DX and D~ x ¼~ x ~ xc . From a physical point of view the center of the RVE is identified as the macroscopic point, at which the macroscopic deformation gradient tensor FM and its gradient 3GM are calculated. The RVE volume represents the underlying microstructure in a finite vicinity of this point. From (4) the microscopic deformation gradient tensor Fm is determined as ~
3 GM þ ðr0m D~ xÞc ¼ FM þ DX wÞc ; Fm ¼ ðr0m D~

ð5Þ 3

where the minor symmetry GMijk = GMkji of the tensor GM, following from its definition (3), has been used. In order to proceed with the formulation of the fine scale boundary value problem some assumptions on the unknown microfluctuation field D~ w are needed. Consider, within a two-dimensional context, an initially rectangular RVE, schematically depicted in Fig. 2. In the undeformed reference state the RVE has width W, ~ L ¼ N ~ R and N ~ B ¼ N ~ T , where the subscripts L, R, B and T denote height H and boundary normals N quantities corresponding to the left, right, bottom and top boundary of the RVE. For this RVE the periodicity conditions on the microfluctuation field are written as wR ðsÞ and D~ wB ðsÞ ¼ D~ wT ðsÞ; D~ wL ðsÞ ¼ D~

ð6Þ

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5529

Fig. 2. Schematic picture of an undeformed RVE for the second-order computational homogenization.

where s is a local coordinate along the edge. It is important to notice that, although the concept of a periodic fluctuation field according to (6) is applied in the present approach, the overall RVE deformation is not periodic. In the general case of a non-zero gradient of the deformation gradient tensor, 3GM 5 30, a non-periodic RVE shape will be found. Therefore, the notion of periodicity is employed here in a ‘‘generalized’’ sense, a concept that will be clarified in the following. With account for the periodicity conditions (6) it can be shown straightforwardly that the macroscopic deformation tensor FM equals the volume average (over the undeformed volume V0 of the RVE) of the microstructural deformation tensor Fm defined by (5) Z 1 FM ¼ Fm dV 0 : ð7Þ V 0 V0 The kinematical scale transition relation (7) is one of the basic averaging theorems of the classical homogenization theory (see e.g. [35,3,36,14]), which is also preserved in the present second-order computational homogenization framework. As has been shown in [32–34], in order to fully impose the macroscale kinematics (defined by the macroscopic deformation gradient tensor FM and its gradient 3GM) on a microstructural RVE, additional constraints on the microfluctuation field are required. Again, different choices are possible. In the present framework the required constraints are obtained by establishing a relation between the macroscopic gradient of the deformation 3GM and microstructural quantities. Note, that until now only relation (7) between the macroscopic deformation tensor FM and its microstructural counterpart Fm has been considered. The following averaging theorem is introduced [33,34], relating the macroscopic gradient of the defor~ in the initial mation 3GM to the fine scale variables (represented by the microstructural position vectors X configuration and by ~ x in the deformed configuration, along with the shown Lagrangian gradients)
 Z     W2 3 1 1 3 RC RC ~Þ þ r0m ðD~ ~Þ C dV 0 ; GM þ ðII : GM Þ r0m ðD~ xX xX ð8Þ ¼ 2 2V 0 V 0 12 where I is the second-order unit tensor; conjugation for a third-order tensor 3T is defined as T Cijk ¼ T kji . The superscript RC indicates right conjugation, i.e. T RC ijk ¼ T ikj , while the superscript LC will denote left conjugation, T LC ijk ¼ T jik . Relation (8) has been obtained for a square RVE (i.e. assuming H = W). A more general relation may be found in [33,34]. Elaborating the integral in the right-hand side of (8) with the use of the divergence theorem and relation (4), leads to a new constraint on the microfluctuation field along the undeformed boundary C0 of the RVE [33,34]

5530

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

Z

~ D~ ~þX ~D~ ~ Þ dC0 ¼ 3 0: ðN wX wN

ð9Þ

C0

With account for the periodicity conditions (6), the constraint (9) reduces to two constraints along two contiguous boundaries Z Z D~ wL dC0 ¼ ~ 0 and D~ wB dC0 ¼ ~ 0: ð10Þ C0L

C0B

Comparable constraints are then automatically satisfied on the remaining opposite boundaries through the periodicity constraint (6). From (10) it is obvious that these additional constraints enforce the shape of the boundary to approximate the kinematically fully prescribed boundary in an average sense. 2.3. Stresses The determination of the RVE averaged stress measures is based on the Hill–Mandel macrohomogeneity condition [35,3]. This condition requires the microscopic volume average of the variation of work performed on an RVE to equal the local variation of the work on the macroscale. Taking into account that on the macrolevel the modelling deals with a full second gradient continuum (see e.g. [15] for the definition of the local work in this case), the extended Hill–Mandel condition becomes 1 V0

Z

. Pm : dFcm dV 0 ¼ PM : dFcM þ 3 QM .. d3 GM :

ð11Þ

V0

It is remarked that (11) is in fact the definition of the macroscopic first Piola–Kirchhoff stress tensor PM and the macroscopic higher-order stress tensor 3QM. In order to obtain the macroscopic stress measures, first the divergence theorem, with incorporation of the microstructural equilibrium equation (Eq. (14)), is applied to the left-hand side of (11), followed by the substitution of the variation of relation (4). Finally, since it can be proven that there is no contribution of the microfluctuation field to the averaged microstructural work, the following micro–macro scale transition relations for the stress quantities are obtained (see [33,34] for the detailed derivation) Z Z 1 ~ dC0 ¼ 1 ~ PM ¼ Pm dV 0 ; ð12Þ pX V 0 C0 V 0 V0 3

QM ¼

1 2V 0

Z

~pX ~ dC0 ¼ 1 X~ 2V 0 C0

Z

~þX ~Pm Þ dV 0 ; ðPcm X

ð13Þ

V0

~
Pc represents the microstructural first Piola–Kirchhoff stress vector on the RVE boundary. where ~ p¼N m The boundary integrals allow the computation of the macroscopic first Piola–Kirchhoff stress tensor PM and the macroscopic higher-order stress tensor 3QM based on microstructural variables defined on the RVE boundary only. The volume integrals in (12) and (13) provide a clear interpretation of the micro– macro averaging relations for the stress quantities in the second-order computational homogenization framework. As it is the case in the classical (first-order) homogenization, the macroscopic stress tensor PM equals the volume average of the microscopic stress tensor Pm. The macroscopic higher-order stress tensor 3QM can be interpreted as the first moment (with respect to the RVE center) of the microscopic first Piola–Kirchhoff stress tensor Pm over the initial RVE volume V0.

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5531

2.4. Microscopic boundary value problem The problem on the microstructural (RVE) level is formulated as a standard problem in quasi-static continuum solid mechanics. In the absence of body forces the equilibrium equation for the microstructural RVE in terms of the first Piola–Kirchhoff stress tensor Pm (related to the Cauchy stress tensor by Pm ¼ detðFm Þrm
ðFcm Þ1 ), has the form 0: r0m
Pcm ¼ ~

ð14Þ

The material behaviour of each microstructural constituent a (e.g. matrix, inclusion etc.) is described by its constitutive law, specifying a time and history dependent stress–strain relationship ðaÞ ðtÞ ¼ FðaÞ fFm ðsÞ; s 2 ½0; tg: Pm

ð15Þ

It is emphasized that the present framework is not specifically designed for any particular constitutive law; the microstructural material behaviour may be very complex and include a physical and/or geometrical evolution of the microstructure (e.g. phase decomposition, phase transformations, decohesion, delamination, inter or intragranular fracture etc.). The actual position vector ~ x of a point in the microstructural RVE is obtained according to (4). For a given macroscopic deformation tensor FM and its gradient 3GM, the first two (macroscopic) terms in the right-hand side of (4) are readily calculated for every microstructural material point. Thus, the only unknown contribution to the microstructural displacement field comes from the microfluctuation field D~ w. Adding boundary conditions for the microfluctuation field, i.e. the periodicity conditions (6) and zero averaged fluctuations along two adjoint boundaries (10) completes the microstructural boundary value problem (14) and (15). The solution of this boundary value problem may be obtained directly in terms of the unknown microfluctuation field. This is sometimes done in the framework of a first-order homogenization scheme [10]. However, if one wishes to use a general purpose numerical code for the solution of boundary value problems, it is more convenient to reformulate the boundary conditions (6) and (10) in terms of the total microstructural displacements, without separation of the coarse and fine scale contributions. This is done in the remainder of this section. First, relations equivalent to the periodicity conditions are derived. Applying (4) to the left and right (and similarly to the bottom and top) boundaries of the RVE and subtracting the results, with account for the periodicity conditions (6), eliminates the unknowns ~ xc and D~ w and leads to kinematic constraints between opposite boundaries according to ~L þ ~ ~ xR ¼ ~ xL þ A
X a;

~B þ ~ ~ xT ¼ ~ xB þ B
X b;

ð16Þ

with 2

~ R
3 GM
N ~R þ W N ~R ; ~ a ¼ W FM
N 2 2 ~ T
3 GM
N ~T þ H N ~T : ~ T
3 GM ; ~ b ¼ H FM
N B ¼ HN 2 xR ¼ ~ xL þ~ x2 ~ x1 and ~ xT ¼ ~ xB þ~ x4 ~ x1 (with ~ xi It is remarked that with 3G = 30, relations (16) lead to ~ denoting the current position vector of node i, Fig. 2), which are precisely the standard periodic boundary conditions on a two-dimensional periodic RVE used in the classical (first-order) homogenization procedure. Note that the microfluctuation vectors of the corner nodes D~ wi , i ¼ 1; 4 are constrained to be all identical in order to satisfy the periodicity conditions (6). Considering that the fluctuation field on the RVE is in fact determined up to a constant value, expressing a (non-relevant) rigid body displacement, which may be ~ R
3 GM ; A¼ WN

5532

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

suppressed, the fluctuation vectors of the corner nodes are set to zero, so that the positions of the four corner nodes are fully prescribed ~i þ 1DX ~i
3 GM
DX ~i ; D~ xi ¼ FM
DX 2

i ¼ 1; 4:

ð17Þ

The constraints (10), prescribing zero averaged fluctuations along two adjoint boundaries, can be rewritten in terms of the total position vectors of the boundary points in the deformed state as Z Z Z ~L dC0 þ 1 3 GLC : ~L X ~L dC0 ; X X D~ xL dC0 ¼ FM
M 2 C0L C0L C0L Z Z Z ð18Þ ~B dC0 þ 1 3 GLC : ~B X ~B dC0 : X X D~ xB dC0 ¼ FM
M 2 C0B

C0B

C0B

2.5. The two-scale problem The two-scale problem outlined in the previous sections and schematically illustrated in Fig. 1 will be tackled through a nested solution procedure. In the present work the finite element method is used for both (micro and macro) scales. However, other numerical techniques for the solution of boundary value problems may be employed, e.g. the boundary element method, meshless methods, etc. In the following, essential parts of the two-scale nested solution scheme are elaborated. First, in Section 3, the discretized microstructural boundary value problem is considered. Then, in Section 4, the micromacro scale transition relations are established. The numerical procedures for computation of the macroscopic stress and the macroscopic constitutive tangents based on microstructural quantities are presented. The next step in the two-scale scheme is the solution of the macroscopic problem, which, in the present second-order computational homogenization scheme, is formulated in the context of a full second gradient continuum theory. In Section 5 the second gradient continuum problem is outlined and its finite element implementation is presented. Finally, the nested second-order multi-scale computational scheme is summarized in Section 6.

3. Solution of the microscale problem The microstructural boundary value problem outlined in Section 2.4 is a standard non-linear quasi-static boundary value problem for a classical continuum. As it was the case in the previous section, here the attention remains focused on the two-dimensional rectangular RVE, schematically depicted in Fig. 2. Using a standard finite element procedure, the weak form of the equilibrium for the microlevel RVE (after discretization) leads to a system of non-linear algebraic equations in the nodal displacements u f int ðu Þ ¼ f ext ; 





ð19Þ

expressing the balance of internal and external nodal forces. This system has to be completed by boundary conditions. Hence, the application of the earlier introduced boundary conditions (16)–(18) have to be elaborated in more detail. First, the unknown position of the center of the RVE in the deformed state ~ xc (used in D~ xi ¼ ~ xi ~ xc ; i ¼ 1;4Þ has to be eliminated from (17). This is easily achieved by spatial fixation of the displacement of one of the corner nodes, say i = 1, (thereby suppressing rigid body translation of the RVE). ~c ¼ ~ Taking into account that X 0, ~ xc is obtained as ~1 : ~1  FM
X ~1  1 3 GLC : X ~1 X ~ xc ¼ X M 2

ð20Þ

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5533

Then the displacements of the corner nodes can be written as ~j  X ~1 Þ; ~j  X ~1 Þ þ 1 3 GLC : ðX ~j X ~1 X ~ uj ¼ ðFM  IÞ
ðX M 2 ~ u1 ¼ ~ 0:

j ¼ 2; 3; 4

ð21Þ

For a given macroscopic deformation gradient tensor FM and its gradient 3GM and a specified geometry of the RVE, the displacements of the corner nodes are now easily calculated and added to the system (19). For the initially rectangular RVE depicted in Fig. 2, the constraints (16), reflecting a linear relation between opposite edges, can be recast in terms of displacements as ~ uR ¼ ~ uL þ 12ð1  gÞð~ u2  ~ u1 Þ þ 12ð1 þ gÞð~ u3  ~ u4 Þ; ~ uB þ 12ð1  nÞð~ u4  ~ u1 Þ þ 12ð1 þ nÞð~ u3  ~ u2 Þ; uT ¼ ~

ð22Þ

with g and n the local coordinates of the corresponding points on the right–left and top–bottom boundaries, respectively (see Fig. 2), and the vectors ~ ui ; i ¼ 1; 4 are given by (21). Next, for purely non-restrictive practical reasons, it is assumed that the finite element discretization is performed such that in the undeformed state the distribution of the nodes on opposite RVE edges is equal. In a discretized format the relations (22) then easily lead to a set of homogeneous constraints (tyings) of the type C a ua ¼ 0;

ð23Þ

with Ca a matrix containing coefficients in the constraint relations and ua a column with the degrees of freedom involved in the constraints. Boundary constraints (18) prescribing zero-averaged microstructural fluctuations may be rewritten in terms of displacement vectors of the nodes on the boundary nL X

ai~ uL ðFM ; 3 GM Þ; uiL ¼ ~

i¼1

nB X

bi~ uB ðFM ; 3 GM Þ; uiB ¼ ~

ð24Þ

i¼1

where nL and nB are the numbers of nodes on the left and bottom boundary, respectively; ai and bi are coefficients following from the discretized form of the integrals in the left-hand sides of (18); ~ uL and ~ uB are vectors from the right-hand sides of (18). Eliminating ~ xc according to (20) gives for ~ uL and ~ uB  Z Z ~L  X ~1 Þ dC0 ; ~L  X ~1 Þ dC0 þ 1 3 GLC : ~L X ~1 X ~ ðX ðX uL ¼ ðFM  IÞ
M 2 C0L C0L Z Z ð25Þ ~B  X ~1 Þ dC0 : ~B  X ~1 Þ dC0 þ 1 3 GLC : ~B X ~1 X ~ ðX ð X uB ¼ ðFM  IÞ
M 2 C0B

C0B

Consequently, the vectors ~ uL and ~ uB are known for any given FM and 3GM and RVE geometry. Contrary to the constraints (22), which gave rise to the homogeneous tying relations (23), the constraints (24) result in a set of non-homogeneous constraints of the type 

C b ub ¼ qb : 



ð26Þ

Procedures for imposing the constraints (23) and (26) include the direct elimination of the dependent degrees of freedom from the system of equations, or the use of Lagrange multipliers or penalty functions. In the following, the constraints (23) and (26) are enforced by elimination of the dependent degrees of freedom. Although such a procedure may be found in many textbooks on finite elements (e.g. [37]), here it is summarized for the sake of completeness but also in the context of the derivation of the macroscopic tangent stiffness, to be presented in Section 4.2.

5534

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

First, the constraints (23) and (26) are combined into Cu ¼ q; 

ð27Þ



where the column q has non-zero values only at the positions corresponding to the non-homogeneous con straints (26). In the following derivations no distinction is made between constraints with zero and with non-zero values in q and the right-hand side column q in (27) is treated as a column containing given values.   Next, (27) is partitioned according to " # u i ½ Ci Cd  ð28Þ ¼ q; u  d where ui are the independent degrees of freedom (to be retained in the system) and ud are the dependent degrees of freedom (to be eliminated from the system). Because there are as many degrees of freedom ud as there are independent constraint equations in (28), matrix Cd is square and non-singular. Solving for u yields d " # u   i 1 1 ¼ ð29Þ u C C d di d  q ; with C di ¼ C d C i : 

This relation may be further rewritten as " # " #
 u u I O i i ; ¼ T q ; with T ¼ u C di C 1 d d 

ð30Þ

where I is a unit matrix of size [Ni · Ni] and O is a zero matrix of size [Ni · Nd], with Ni and Nd the number of the independent and dependent degrees of freedom, respectively. With the transformation matrix T a linear system of equations K d ¼ r is transformed into a new system K 0 d0 ¼ r 0 using d ¼ T d0 , r 0 ¼ T Tr and K 0 = TTK T. Standard linearization of the non-linear system of Eq. (19) leads to a linear system in the iterative corrections du to the current estimate u. This system may be partitioned as  2 3 2 3
 dui dri K ii K id  4 5¼4 5; ð31Þ du dr K di K dd d d   with the residual nodal forces at the right-hand side. Noting that all the constraint equations considered above are linear, and thus their linearization is straightforward, application of the transformation (30) to the system (31) gives 3 " #2 du 3 2 dri þ C Tdi drd i K ii þ K id C di þ C Tdi K di þ C Tdi K dd C di ðK id þ C Tdi K dd ÞC 1 d 4 5¼4 5: ð32Þ T T T dq ðC 1 Þ dr d ðC 1 Þ ðK þ K C Þ ðC 1 Þ K C 1 d

di

dd

di

d

dd

d



d



The prescribed degrees of freedom (among which are the prescribed displacements of the corner nodes and the prescribed values of dq) can be further eliminated from the system (32), leaving the displacements of the  independent nodes as the only unknowns in the system.

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5535

4. The scale transition 4.1. Calculation of the macroscopic stress and higher-order stress After the analysis of a microstructural RVE is completed, the RVE averaged stress tensors are extracted and transported to the macroscopic level as the macroscopic stress tensors in the integration point. Of course, the averaged stress and higher-order stress tensors can be calculated by numerically evaluating the volume or surface integrals in (12) and (13). However, for the particular implementation proposed here, computationally more efficient formulas may be obtained. It may be verified that all the forces involved in the homogeneous kinematic constraints (22) cancel out from the surface integrals in (12) and (13), while the forces involved in the non-homogeneous constraint relations (24) have to be properly accounted for. After somewhat lengthy but straightforward mathematical manipulations the surface integrals in (12) and (13) can be transformed into 1 X ~ ~H PM ¼ ð33Þ f X ; i ¼ 1; 2; 3; 4; L ; B V 0 i ðiÞ ðiÞ 3

QM ¼

with ~H X ðiÞ

YH ðiÞ

1 X ~ H LC ; f ðiÞ YðiÞ 2V 0 i

i ¼ 1; 2; 3; 4; L ; B ;

8 ~ð1Þ ; ~ X for X > > < R ðiÞ ~L  X ~ð1Þ Þ dC0 ; for ðX ¼ C0L > R > : ~B  X ~ð1Þ Þ dC0 ; for ðX C0B

ð34Þ

i ¼ 1; 2; 3; 4; i ¼ L ;

ð35Þ



i¼B :

8 ~ X ~ð1Þ ; ~ð1Þ X ~ X for X > > < R ðiÞ ðiÞ ~ ~ ~ ~ ðX L X L  X ð1Þ X ð1Þ Þ dC0 ; for ¼ C0L > > : R ðX ~B  X ~ð1Þ Þ dC0 ; for ~B X ~ð1Þ X C0B

i ¼ 1; 2; 3; 4; i ¼ L ;

ð36Þ



i¼B :

In (33) and (34) ~ f ðiÞ , i ¼ 1; 4 are the external forces in the four prescribed corner nodes and ~ f ðL Þ and ~ f ðB Þ are the resultant forces necessary to enforce the non-homogeneous constraints (24). If a Lagrange multiplier technique would have been used to impose the constraints (24), the forces ~ f ðL Þ and ~ f ðB Þ might be identified as the Lagrange multipliers. If elimination of the dependent degrees of freedom, as outlined in the previous section, is used, these forces are extracted from the column in the right-hand side of (32). In a finite element program the forces ~ f ðiÞ , i = 1, 2, 3, 4, L*, B* are readily available for the converged solution, from which the macroscopic stress tensor PM and the higher-order stress tensor 3QM can be easily calculated using the formulas (33) and (34). 4.2. Macroscopic constitutive tangents For the finite element solution of the macroscopic problem a stiffness matrix at every macroscopic integration point is required. In computational homogenization schemes there is no explicit form of the macroscopic constitutive behaviour assumed a priori. Therefore the tangent operator has to be determined numerically. This has been realized by some authors by numerical differentiation based on a perturbation technique [38]. An alternative approach is the condensation of the total microstructural stiffness to the local macroscopic stiffness. In the framework of a first-order computational homogenization scheme such a

5536

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

procedure, in combination with the Lagrange multiplier method to impose boundary constraints, has been recently elaborated in [9,39]. An alternative and more efficient scheme for the first-order homogenization, which employs the direct condensation of the constrained degrees of freedom has been presented in [14]. In the present paper, the latter procedure is extended to obtain constitutive tangents for a second gradient continuum in the framework of the presented second-order computational homogenization scheme. In order to derive the macroscopic stiffness matrix, the system of Eq. (32) is written as K H du ¼ df H and   rearranged to the form "

KH pp KH fp

3 2 H3 2 H3 #2 df p df p du  p 4 5¼4 5  4  5; du df fH KH 0 ff f

KH pf

ð37Þ





where the subscript p refers to ‘‘prescribed’’ degrees of freedom (degrees of freedom through which the macroscopic tensors FM and 3GM are imposed on the RVE). The subscript f refers to all remaining ‘‘free’’ nodes. System (37) is taken at the converged end of a microstructural increment, thus the residual forces in the free nodes can be neglected df fH  0. Elimination of du from the system (37) then leads to the ref  duced stiffness matrix K H that relates the variations of the prescribed degrees of freedom to the variations M of the associated forces KH ¼ df pH ; M du p 

1

H H H H with K H M ¼ K pp  K pf ðK ff Þ K fp :

ð38Þ

Next, relation (38) needs to be transformed to obtain an expression relating variations of the macroscopic stress tensors to variations of the associated kinematical quantities. Since in the second-order computational homogenization framework the macrostructure is modelled as a full second gradient continuum, the linearized constitutive relations can always be written in the form ð1Þ ð2Þ . dPM ¼ 4 CM : dFcM þ 5 CM .. d3 GRC M ;

ð39Þ

ð3Þ ð4Þ . d3 QM ¼ 5 CM : dFcM þ 6 CM .. d3 GRC M ; 4 ð1Þ CM ,

ð40Þ 5 ð2Þ CM

5 ð3Þ CM

ð4Þ

the fifth-order tensors and and the sixth-order tensor 6 CM where the fourth-order tensor are the macroscopic consistent constitutive tangents. In order to obtain these constitutive tangents departing from the reduced matrix K H M , relation (38) is first rewritten in a vector/tensor format X ðijÞ KM
d~ uðjÞ ¼ d~ f ðiÞ ; i; j ¼ 1; 2; 3; 4; L ; B ; ð41Þ j ðijÞ

where the components of tensors KM are simply found in the tangent matrix K H M at the rows and columns of the degrees of freedom corresponding to i and j. Next, the expression for the variation of the forces (41) is substituted into the relations for the variations of the macroscopic stress and higher-order stress obtained by varying (33) and (34), which leads to 1 X X ðijÞ ~H ; dPM ¼ ðKM
d~ uðjÞ ÞX ð42Þ ðiÞ V0 i j d3 QLC M ¼

1 X X ðijÞ ðKM
d~ uðjÞ ÞYH ðiÞ : 2V 0 i j

ð43Þ

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5537

The vectors d~ uðjÞ are now obtained from (21) and (25) as H

~
dFc þ 1YH : d3 GRC : d~ uðjÞ ¼ X M M ðjÞ 2 ðjÞ Substitution of (44) into (42) and (43) gives

1 X X ~H ðijÞ ~H LC 1 ~H ðijÞ H LC .. 3 RC . d X ðiÞ KM X ðjÞ dPM ¼ : dFcM þ X K Y G ; ðiÞ M ðjÞ M V0 i 2 j d3 Q M ¼



C23 . 1 X X  H ðijÞ ~H C23 1 ðijÞ H .. d3 GRC ; YðiÞ KM X ðjÞ : dFcM þ YH K Y M 2V 0 i 2 ðiÞ M ðjÞ j

ð44Þ

ð45Þ

ð46Þ

where the superscript C23 indicates conjugation on the second and third indices. Comparing (45) and (46) with (39) and (40), the consistent tangents are finally identified as 1 X X ~H ðijÞ ~H LC 4 ð1Þ X ðiÞ KM X ðjÞ CM ¼ ; V0 i j 1 X X ~H ðijÞ H LC 5 ð2Þ X ðiÞ KM YðjÞ CM ¼ ; 2V 0 i j ð47Þ 1 X X  H ðijÞ ~H C23 5 ð3Þ CM ¼ YðiÞ KM X ðjÞ ; 2V 0 i j C23 X X 1 ð4Þ ðijÞ H 6 CM ¼ YH : ðiÞ KM YðjÞ 4V 0 i j Thus, the constitutive tangents of the macroscopic second gradient continuum are directly obtained through static condensation of the microscopic global stiffness matrix. Notice that consistency is preserved through this micro–macro transition (which is only important if a consistent microstructural tangent operator is used). In a non-linear case the obtained macroscopic tangents automatically include both geometrical and material contributions, evidently assuming that these parts have been properly dealt with on the microscale. Use of a multi-scale procedure as described above overcomes one of the major difficulties in higher-order constitutive modelling: identification and quantification of many constitutive parameters. When a linear elastic microstructural material behaviour is considered, this scheme leads to a direct quantification of macroscopic elastic stiffnesses (e.g. of MindlinÕs type) based on a microstructural analysis. In geometrically and physically non-linear regimes such an approach is particularly valuable, since in these cases it is especially difficult to make a well-motivated assumption on the closed-form of the (evolving) underlying higher-order constitutive relation.

5. Solution of the macroscale problem 5.1. Second gradient boundary value problem Contrary to the microstructural problem, which remains a classical one, on the macroscale an equilibrium problem for a full second gradient continuum appears. Such a continuum has been studied in the context of elasticity in the sixties [40–43] and used in several strain gradient plasticity theories later [15,44]. In a second gradient continuum theory the internal work is assumed to be determined by the deformac tion gradient tensor F ¼ ðr0~ xÞ and its (Lagrangian) gradient 3G = $0F combined with their respective work conjugated quantities: the first Piola–Kirchhoff stress tensor P and the higher-order stress tensor

5538

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

3

Q. Starting from this assumption after some mathematical manipulations the strong form of local equilibrium for a second gradient continuum can be derived [43,15,34]   r0
Pc  ðr0
3 QÞc þ ~ b ¼~ 0; ð48Þ where ~ b is the body force per unit undeformed volume. The natural boundary conditions associated with this system of partial differential equations are expressed in (i) the surface traction ~ t  c  c c c s s 3 3 3 ~ ÞN ~
ðN ~
QÞ  r
ðN ~
QÞ ; ~
P  ðr0
QÞ þ ðr
N ~ ð49Þ t¼N 0 0 ~N ~ Þ
r0 , and (ii) the double stress where the surface gradient operator is defined as rs0 ¼ ðI  N traction ~ r ~: ~
3Q
N ~ r¼N

ð50Þ

In the case that the surface of the body is non-smooth (has edges) also an additional line load appears. The kinematic boundary conditions for the second gradient continuum include prescribed displacements ~
r0 . Mixed boundary conditions are also ~ u and normal gradients of displacements D0~ u with D0 ¼ N possible. In the above relations and throughout this section the subscript ‘‘M’’ is omitted, since all quantities in this section pertain to the macroscopic level. Finally, constitutive equations relating the first Piola–Kirchhoff stress tensor P and the higher-order stress tensor 3Q to the history of the deformation tensor F and its gradient 3G, should be added. The constitutive relations in a general format may be written as ~; tÞ ¼ FP fFðX ~; sÞ; 3 GðX ~; sÞ; s 2 ½0; tg; PðX 3

~; tÞ ¼ FQ fFðX ~; sÞ; 3 GðX ~; sÞ; s 2 ½0; tg: QðX

ð51Þ ð52Þ

The main advantage of the micro-macro framework presented is that the required macroscopic constitutive response is obtained directly from the numerical analysis of the underlying microstructure. No assumptions on the form of the higher-order constitutive relations (51) and (52) need to be made. Actually this is one of the major difficulties in the closed-form constitutive modelling of higher-order continua, where many constitutive parameters need to be identified and quantified. If additionally, the non-linear microstructural behaviour is strongly evolving, this task becomes nearly impossible. Using the computational homogenization strategy this problem is easily overcome. Moreover, for the nested solution of the resulting non-linear macroscopic boundary value problem macroscopic tangent operators are required. In the computational homogenization approach the tangents are obtained directly through static condensation of the microscopic global stiffness matrix. This procedure was presented in detail in Section 4.2. 5.2. Finite element implementation The direct application of the conventional displacement-based finite element method to solve the second gradient boundary value problem (48)–(52) requires at least C1 continuous displacement interpolation functions, i.e. both the displacements and their first derivatives are required to be continuous across interelement boundaries. Several C1 continuous elements have been proposed for higher-order continua (e.g. [45–47]). However, development of C1 continuous elements suitable for general multi-dimensional nonlinear problems is not straightforward, since these elements usually require additional assumptions (e.g. geometrical) and hence suffer from various restrictions. Therefore, the use of C0 continuous elements is here preferred. The successful use of C0 continuous elements for higher-order continua is based on a mixed (multi-field) formulation [48,45,49–51]. The strategy is to introduce another unknown field (in this case a

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5539

b in addition to the unknown displacement field ~ second-order tensor field, denoted by F) u. The kinematic constraint between these two fields is weakly enforced by means of Lagrange multipliers, also forming a second-order tensor field. The weak form, which accounts for the appropriate field equations, (natural) boundary conditions and c b ¼ F in the volume, enforced by the Lagrange multiplier k, and F b s ¼ ðrs~ the relationships F 0 xÞ on the surbs ¼ F bF b
N b enforced by another Lagrange multiplier l, ~N ~ the surface part of the tensor F), face (with F may be written as Z Z   

 . b : dFc þ 3 QðF; 3 GÞ b .. d3 G b þ d k : ðF b c  rs~ b c  Fc Þ dV 0 þ d l : ð F PðF; 3 GÞ 0 xÞ dS 0 s V0

¼

Z

~ b
d~ u dV 0 þ

V0

Z

~ t
d~ u dS 0 þ

S0

Z

S0 c

b dS 0 ; R : dF

ð53Þ

S0

Here for simplicity it is assumed that the body has a smooth surface, so that no contributions related to line forces on edges are present. It is a mathematically straightforward exercise to show the equivalence between (53) and the strong form (48)–(52). In (53) the deformation gradient tensor F is obtained from the displacec c b is a third-order ment field ~ u via the standard kinematics relation F ¼ ðr0~ xÞ ¼ I þ ðr0~ uÞ . The variable 3 G b tensor defined as the gradient of the field F, symmetrized on the first and the last indices b ijk ¼ 1 ðr0i Fb jk þ r0 Fb ji Þ. In the right-hand side of (53) R is a double traction tensor, which is related G k 2 ~ . In the weak form (53) the dependence of the stress tensor to the double traction vector ~ r by R ¼ ~ rN and the higher-order stress tensor on kinematic quantities has been explicitly indicated. This choice has been made in the present implementation, however, other alternatives are possible. b ¼ F is satisfied in the volume in an averaged sense, the surface integral in the leftAssuming that F hand side of (53) may be neglected. This allows to rewrite the weak form as Z Z Z Z  

.. 3 c c c 3 b b c dS 0 ; ~ b P : dF þ Q . d G þ d k : ð F  F Þ dV 0 ¼ t
d~ u dS 0 þ R : d F b
d~ u dV 0 þ ~ V0

V0

S0

S0

ð54Þ which for a sufficiently fine mesh does not introduce a significant error, while it simplifies the finite element implementation. The same assumption has also been made by other authors ([49,51]). However, in [50] the surface integral term has been retained. In order to correctly formulate a boundary value problem for a smooth second gradient continuum body, six independent boundary conditions are required for a three-dimensional problem (four independent boundary conditions for a two-dimensional problem). A way to prescribe the boundary conditions associb has been addressed in [49]. Let a superscript * indicate a ated with the primary unknown fields ~ u and F  prescribed quantity on the surface. If the stress traction ~ t and the double stress traction ~ r  are given, the boundary conditions to be prescribed are simply  ~ t ¼~ t

and

~: R ¼~ r N

ð55Þ 

If the displacement ~ u  and the displacement gradient normal to the surface ðD0~ uÞ are given, the boundary conditions are prescribed as ~ u ¼~ u

and

 bc ¼ I þ N ~ ðD0~ F uÞ þ rs0~ u ;

ð56Þ

where the surface gradient rs0~ u  is given on the surface. Using these u  is known because the displacement ~ kinematic boundary conditions no errorc is introduced in (54) upon neglecting the boundary integral term in b ¼ rs~ (53), since in this case the constraint F 0 x is automatically satisfied by (56). Mixed boundary condis tions are obtained by combination of the dynamic and kinematic boundary conditions (55) and (56).

5540

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

b and their corresponding weighting Only the first-order derivatives of the primary unknown fields ~ u and F b appear in the weak form (53) (or (54)), so they are only required to be piece-wise diffunctions d~ u and d F ferentiable, i.e. C0 continuous. No derivatives of the Lagrange multipliers k are present in this weak form, ~, i.e. C1 continuous so the Lagrange multipliers should only be piece-wise continuous functions of X functions. The finite element approximations of the unknown fields are written as b ¼ NF F b F A A;

~ x ¼ N xA~ xA ; x

F

k ¼ N kA kA ;

ð57Þ

0

where N and N are C continuous interpolants for the position vector field and the relaxed deformation gradient field, respectively, and Nk are (discontinuous across elements) interpolants for the Lagrange multiplier field. Here and in the following capital Latin subscripts indicate nodal quantities, while small Latin subscripts will refer to components of vectors and tensors. The Einstein summation convention over repeated indices is used in both cases. The range of the capital Latin letters covers all the nodes (of the mesh or the element, depending on the context); the range of the small Latin letters equals the number of dimensions of the considered configuration. Following the standard Galerkin approach the interpolants Nx, NF and Nk are also used for the respecx BA and tive weighting functions. The discrete approximations to the gradient relations read F ¼ ~ xA ~ F 3b b A , with G ¼~ BA F  F F x oN xA F b ijk ¼ 1 oN A Fb Ajk þ oN A Fb Aji : ~ BA ¼ and ~ BA defined such that G ~ 2 oX i oX k oX Discretizing the weak form (54) and taking into account that it must be satisfied for arbitrary weighting parameters results in a system of equations constituting the discrete force balances and additionally the discretized kinematic constraint Z Z Z x x ~ ~ ~ P
BA dV 0 þ kB
C AB ¼ ð58Þ bN A dV 0 þ ~ tN xA dS 0 ; V0

Z

V0 3

F Q
~ BA dV 0 þ kcB EAB ¼

V0

Z

S0

Rc N FA dS 0 ;

ð59Þ

S0

b c ¼ 0; ~ C AB~ xA þ EAB F A

ð60Þ

where the following notation has been introduced Z Z x k ~ ~ C AB ¼  BA N B dV 0 ; EAB ¼ N FA N kB dV 0 : V0

ð61Þ

V0

Due to the presence of geometrical and/or material non-linearities the system of Eqs. (58)–(60) is generally non-linear. The solution of this set of equations is determined iteratively using the Newton–Raphson procedure. The non-linear system of equations is linearized in a standard manner. Consistent linearization of the constitutive relations (51) and (52) leads to . b RC ; DP ¼ 4 Cð1Þ : DFc þ 5 Cð2Þ .. D3 G

ð62Þ

. b RC ; D3 Q ¼ 5 Cð3Þ : DFc þ 6 Cð4Þ .. D3 G

ð63Þ

3

where DP and D Q are the iterative corrections for the stress tensor and the higher-order stress tensor, b are calculated from the iterative corrections D~ b respecrespectively. In (62) and (63) DF and D3 G x and D F,

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5541

tively. It is important to note that in the computational homogenization framework presented here, the consistent tangent stiffnesses 4C(1), 5C(2), 5C(3), 6C(4) are obtained numerically, based on the microstructural behaviour, according to (47). The total number of scalar unknowns in the linearized system of equations in the 2D case is equal to 2Kx + 4KF + 4Kk, with Kx, KF and Kk the number of nodes in the mesh for the respective unknowns. This is significantly larger compared to a standard displacement-based finite element formulation for a classical continuum. However, taking advantage of the weaker continuity requirements on the Lagrange multiplier k, allows to eliminate k at the element level prior to the assembly. This reduces the total number of unknowns in the assembled system of equations. For this, however, a penalized version of the weak form has to be used

Z   . b þ d k : ðF b c  Fc Þ þ 1 k : dkc dV 0 P : dFc þ 3 Q .. d3 G a V0 Z Z Z b c dS 0 ; ~ ¼ ð64Þ t
d~ u dS 0 þ R : d F b
d~ u dV 0 þ ~ V0

S0

S0

where a is a penalty parameter. As usual for large deformation problems, the external load is applied in a number of loading steps (increments). Besides, it is worth mentioning that in the previous derivations the total Lagrangian formulation has been adopted, i.e. all the derivatives and integrals are taken with respect to the Lagrangian (material) ~. The same theory can also be formulated in terms of an updated Lagrangian procedure, as coordinates X has been presented in [52,51]. The updated Lagrangian formulation has an advantage over the total Lagrangian formulation when implementing rate-type constitutive models. However, the actual implementation has been developed for use in the context of a computational homogenization scheme, where the stress updates and the consistent tangents are obtained numerically from the microstructural analysis. As demonstrated in the previous sections this can be done in a convenient way with a total Lagrangian framework. Several plane strain quadrilateral elements for the analysis of a second gradient continua have been developed and implemented. The performance of the elements has been evaluated on a patch test and by comparing the numerical results with the analytical solution for the boundary shear layer assuming MindlinÕs second-order elastic material model. Two elements that successfully pass the tests are shown in Fig. 3. The element with nine displacement nodes, nine deformation gradient nodes and four Lagrange multiplier nodes located in the 2 · 2 Gauss quadrature points (Fig. 3a) shows very good convergence of the numerical solution towards the analytical one (upon mesh refinement). However the obtained system of equations for this case is about two times larger than for the other element considered. Moreover, in order to avoid zero energy modes, this element should be used only in combination with a 3 · 3 Gauss integration

Fig. 3. Finite elements for the second gradient continuum formulation. (d) displacements u1, u2; (h) deformation gradients F^ 11 , F^ 22 , F^ 12 , F^ 21 ; (·) Lagrange multipliers k11, k22, k12, k21.

5542

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

scheme. This latter observation is rather crucial for the selection of an element in a computational homogenization analysis, where generally a microstructural analysis has to be performed at every macroscopic integration point. The use of four integration points instead of nine integration points would decrease the number of required microstructural analyses more than twice. Therefore, in the subsequent coupled second-order homogenization computations the element shown in Fig. 3b will be used. This element has eight displacement nodes and four deformation gradient nodes, while the Lagrange multipliers are evaluated in the center of the element. In most problems this element performs well with a 2 · 2 Gauss integration scheme. More details on the finite element implementation of the second gradient equilibrium and the element selection for the computational homogenization framework are beyond the scope of this paper and may be found in [34].

6. The nested multi-scale solution scheme The direct application of the presented second-order computational homogenization framework, is next performed through a nested solution scheme for the coupled multi-scale numerical analysis. The macroscopic structure to be analyzed is discretized by finite elements. To each macroscopic integration point a unique microstructural RVE is assigned. Note that in general, it suffices to make this scale jump only in those zones that are highly non-linear or evolving, i.e. where closed-form constitutive equations become clearly inaccurate. In order to initiate the macroscopic finite element analysis, constitutive tangents at every integration point are required. To obtain these tangents from the microstructural properties a preparing microstructural analysis is performed. During this initialization the stiffness matrix of an undeformed RVE is assembled and used to derive the initial macroscopic constitutive tangents at a macroscopic integration point. During the actual analysis the external macroscopic load is applied in increments. For every step of the macroscopic incremental-iterative procedure, and in each macroscopic integration point, the macroscopic deformation gradient tensor FM and its gradient 3GM are calculated based on the current (iterative) macroscopic displacement field. These deformation tensors and gradients of deformation tensors are transferred to the microlevel, where they are used to construct the boundary value problem for the RVE, corresponding to the respective macroscopic integration point (see Section 2.4). Upon the solution of every RVE problem, the averaged stress tensor PM and the higher-order stress tensor 3QM are obtained using (33) and (34). Additionally, the constitutive tangents are extracted according to (47) and returned to the macroscopic FE program. When the analysis of all RVEs is finished, the stress tensor, the higher-order stress tensor and the consistent constitutive tangents are available at every macroscopic integration point. Hence, macroscopic internal nodal forces can be calculated, higher-order equilibrium can be evaluated and, if required, the next macroscopic iteration can be performed. If equilibrium is achieved the calculations can be continued for the next increment. This solution scheme is summarized in Table 1. Obviously, the multi-scale algorithm described above is parallel by its nature. All RVE calculations for one macroscopic iteration can be performed at the same time without any exchange of data between them. So the use of parallel processors for the RVE analyses significantly reduces the total micro–macro calculation time, a strategy which is adopted in the current implementation.

7. Computational homogenization modelling of boundary shear layers In this section, as an example, the presented second-order computational homogenization approach is applied to the modelling of simple shear of a heterogeneous strip with boundary constraints. Simple shear

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5543

Table 1 Incremental-iterative nested multi-scale solution scheme for the second-order computational homogenization

of a constrained strip is often used in literature as a benchmark problem for higher-order (plasticity) constitutive models, e.g. in [18,53,54]. The particular microstructure considered in this example is of purely academic nature. The boundary constraints, however, can often be found in practice and may be due to, for example, rigid substrates, oxide layers or a hard coating. The problem of a constrained boundary may arise for example in the mechanical response of micro-systems and thin films (e.g. in micro-electronics, MEMS, etc.), for which a sensitivity to the characteristics (and constraints) of boundary or interface layers should be explored. On the macroscopic level, plane strain shear of a two-dimensional strip with a height H in the X2 direction and an ‘‘infinite’’ length in the X1 direction, is considered (see Fig. 4a). In order to model constraints at the (upper and lower) edges, higher-order boundary conditions prescribing vanishing shear at the surface are applied. The boundary conditions in the present analysis are

5544

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550 ^

d

^

d ^

^

(a)

(b)

Fig. 4. The macroscopic model (a) and the microstructural unit cell (b) for the computational homogenization analysis of the boundary shear layer problem.

u1 ¼ 0; u1 ¼ U  ; U*

Fb 12 ¼ 0;

u2 ¼ 0;

c*H,

u2 ¼ 0; c*

Fb 12 ¼ 0;

Fb 22 ¼ 1 on X 2 ¼ 0; Fb 22 ¼ 1

on X 2 ¼ H ;

= with the prescribed shear. Since for this simple shear problem all quantities are indewhere pendent of X1, one single row of two-dimensional plane strain elements, as indicated in Fig. 4a, suffices. Periodicity conditions are applied between the left and right sides of the mesh, tying the associated degrees of freedom. The material considered consists of an elasto-plastic matrix with randomly distributed voids (12% volume fraction). The material parameters of the matrix material correspond to a commercial steel, T67CA, YoungÕs modulus E = 210 GPa, PoissonÕs ratio m = 0.3, initial yield stress ry0 = 507 MPa and (constant) hardening modulus h = 200 MPa. For the microscopic analysis the unit cell plotted in Fig. 4b is used. The average diameter of the voids is 32 lm. The size of the square microstructural unit cell is taken as d = 0.2 mm. In the present calculations the height of the macroscopic layer H has been varied between 1 mm and 10 mm, while keeping the size of the microstructural cell d unchanged. This allows to investigate the influence of the ratio H/d, between the height of the strip and the characteristic size of the microstructure, on the development of the boundary layer. Macroscopic shear loading has been applied up to an average shear of c* = 0.01. When the boundary value problem stated above is modelled in the context of a classical (local) continuum (with the higher-order boundary conditions left out of consideration), the solution is typically given by a homogeneous shear across the layer. The first-order computational homogenization is fully consistent with a local continuum formulation, and thus also predicts a macroscopically homogeneous deformation. The second-order computational homogenization, on the other hand, allows the complete incorporation of the boundary constraints via the higher-order boundary conditions. This leads to a non-uniform ‘‘macroscopic’’ deformation distribution over the strip. Fig. 5a shows the distribution of the shear (shear component F12 of the deformation gradient tensor normalized by the prescribed average shear c*) along the height of the strip for different ratios of H/d at the final shear value c* = 0.01. Boundary layers with a vanishing shear can be clearly observed. The width of the boundary layers reduces with an increasing ratio of H/d. For the smaller values of H/d (H/d = 5 and 10) the boundary layer extends practically over the full thickness of the strip, while for larger values of H/d the zones of reduced shear and of homogeneous deformation are clearly distinguishable.

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5545

dimensionless position, x2 /H

200

shear stress, MPa

150

100 H/d=5 H/d=10 H/d=20 H/d=50

50

0

(a)

12

0

0.002

(b)

0.004 0.006 0.008 averaged shear

0.01

Fig. 5. (a) Shear strain distribution along the height of the strip for several values of H/d at the prescribed average shear c* = 0.01. (b) Overall shear stress–strain response of the strip obtained from the presented second-order computational homogenization model for several values of H/d.

dimensionless position, x2/H

The overall shear stress–strain response obtained from the presented second-order computational homogenization model for different ratios H/d is presented in Fig. 5b. When H/d is reduced from 50 to 5 the elastic stiffness, the yield strength and the hardening rate slightly increase. Thus, the presence of constraint boundaries gives rise to a size effect, which in this case and from a macroscopic point of view can be formulated as: ‘‘smaller is stronger’’. The appearance of a size effect in the elastic regime is inherent to the second-order approach which applies to both the elastic and the elasto-plastic regimes. This clearly differs from some of the gradient plasticity theories e.g. [18], in which the higher-order boundary conditions are related to the plastic deformation only and thus no size effect is taken into account in the elastic region. The development of boundary shear layers is further illustrated in Fig. 6, which shows the shear distributions over the strip (normalized by the averaged shear c*) at four selected values of the prescribed average strain c* in the range of 0.0025–0.01 for a ratio H/d = 20. Clearly, the width of the boundary layer does not

12

Fig. 6. Shear strain distribution along the height of the strip at various values of the average shear c* for a ratio H/d = 20.

5546

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

only depend on the ratio H/d (as was illustrated in Fig. 5a), but also on the level of the prescribed average shear. At small values of the prescribed shear the boundary layer is relatively thin, while it widens as the loading progresses. However, after having reached a certain level of loading the thickness of the boundary layer does not evolve further. The evolution of the shear profiles with ongoing loading may be related to the formation of microstructural shear bands. Fig. 7 illustrates the development of a shear band in the microstructural unit cell located near X2 = H/2 for the case H/d = 20. As can be seen, at small values of the average shear (e.g. c* = 0.0025) the microstructural cell is deformed only slightly and the overall deformation of the cell can be characterized as almost uniform shear (Fig. 7a). This corresponds to the macroscopic shear profile in Fig. 6 with a relatively thin boundary layer and a substantial region of uniform deformation across the macroscopic strip. With increasing macroscopic shear a microstructural shear band starts to develop between the voids. At the macrolevel this corresponds to an extension of the boundary layer. Finally, as the shear band across the unit cell is fully developed, further deformation of the microstructural cell is almost completely realized through plastic yielding within this band. As a result the macroscopic shear profile does not alter its geometry anymore upon further loading (at least not in the range of the macroscopic deformation considered). Fig. 8 shows the deformed macroscopic finite element meshes (displacements are 10 times magnified) and the distribution of the equivalent plastic strain in three microstructural unit cells corresponding to three

Fig. 7. Development of a shear band in the microstructural unit cell located near X2 = H/2 for a ratio H/d = 20 at increasing values of the prescribed macroscopic shear c*; distribution of the equivalent plastic strain. (a) c* = 0.0025; (b) c* = 0.005; (c) c* = 0.0075; (d) c* = 0.01.

Fig. 8. The deformed macroscopic mesh (displacements are 10· magnified) and the distribution of the equivalent plastic strain in the deformed unit cells corresponding to three different points of the macroscopic layer at a prescribed shear c* = 0.01 for different ratios H/d: (a) 5; (b) 50.

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5547

different points of the macroscopic layer. Fig. 8a applies for the case H/d = 5 and Fig. 8b for H/d = 50, both at a prescribed shear c* = 0.01. The constraining conditions along the boundary of the macroscopic strip result in a very small amount of deformation in the microstructural cells close to the boundary. The macroscopic shear is mostly accommodated through the formation of microstructural shear bands in the bulk of the strip. As has been pointed out already, for H/d = 5 the shear profile is non-uniform across the whole height of the strip. At the microscale this is reflected by a different amount of deformation of the unit cells at different macroscopic locations. The microscopic equivalent plastic strain concentration is considerably higher in the unit cell near the center line of the macroscopic strip, than in a unit cell further away from the center (see Fig. 8a). On the other hand, for a ratio H/d equal to 50, the deformation profile shows a region of uniform shear (as can be seen also in Fig. 5a). Consequently, at the microstructural level it can also be observed that the microstructural unit cells within this region are in a similar state of loading. Finally, it is remarked that overall (macroscale) results obtained here (i.e. a size effect and a dependence of the shear strain profiles on a microstructural length scale parameter, as well as the evolution of the deformation profile with increasing prescribed shear) are qualitatively similar to the results presented in literature for this benchmark case obtained using various higher-order plasticity formulations [18,53,54] and discrete dislocation simulations [53]. Whereas it is often claimed that a particular constitutive description is required to capture these qualitative phenomena, the present analysis shows that it is essentially the second-order continuum approximation along with its enriched boundary conditions that governs the obtained results. Despite its simplicity, the benchmark example discussed in this section clearly demonstrates the ability of the second-order computational homogenization scheme to capture higher-order boundary effects (e.g. the effect of a constrained boundary). The added value of the second-order computational homogenization as a natural way to retrieve the higher-order continuum response is also worth mentioning. Moreover, as has been illustrated by this example, a computational homogenization technique allows a direct investigation of the relation between a macroscopically observed response and microstructural phenomena (up to the scale at which they are modelled).

8. Conclusions This paper dealt with a novel second-order computational homogenization procedure with a special focus on implementation issues. The second-order scheme is based on a proper incorporation of the macroscopic gradient of the deformation tensor into the kinematical micro-macro framework. The higher-order constitutive behaviour, which is difficult to capture in a closed-format with its many constitutive parameters unknown, is (numerically) obtained in a natural way through the solution of a microstructural boundary value problem. Moreover, it has been shown that the macroscopic consistent tangent operators, necessary for the finite element solution of the macroscopic problem, can be extracted from the microstructural stiffness matrix. Thus, the behaviour of a non-linear evolving multi-phase microstructure is directly taken into account in the macrostructural analysis. An important feature of the second-order computational homogenization technique is that it incorporates a length scale (being representative for the macroscopic second gradient continuum) via the size of the microstructural representative cell, which is directly related to the size of the microstructural domain in which the basic microstructural deformation mechanisms occur. This allows the analysis of certain phenomena, not addressed within the first-order scheme, such as macroscopic localization [34,55] and size effects. A typical example of the multi-scale analysis of the size effect originating from the presence of a constrained boundary has been given in this paper. It should be mentioned, that, compared to the first-order framework, the only additional computational effort in a second-order computational homogenization is the solution of the higher-order equilibrium problem on the macrolevel, since the microstructural boundary value problem remains classical. Moreover,

5548

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

the calculation time required for a coupled numerical analysis can be substantially reduced through the use of parallel computations. Another option is a selective usage of the multi-scale framework, i.e. where noncritical regions are modelled by continuum closed-form homogenized constitutive relations (or by the constitutive tangents derived from the microstructure in case an incremental update of these tangents can be omitted, e.g. when the material hardly evolves or unloads) while in critical regions a multi-scale analysis of the microstructure is fully performed. Besides the coupled multi-scale analysis of micro–macro structure-property relations of multi-phase materials, the computational homogenization schemes may also be used for material parameter identification based on the reduced microstructural tangent stiffnesses (e.g. for a higher-order elastic material) or verification of microstructurally based higher-order constitutive models.

Acknowledgment This research was carried out under project number ME97020 ‘‘Multi-scale computational homogenization’’ in the framework of the Strategic Research Programme of the Netherlands Institute for Metals Research in the Netherlands (www.nimr.nl).

References [1] P. Ponte Castan˜eda, P. Suquet, Nonlinear composites, Adv. Appl. Mech. 34 (1998) 171–302. [2] J. Fish, K. Shek, M. Pandheeradi, M.S. Shephard, Computational plasticity for composite structures based on mathematical homogenization: theory and practice, Comput. Methods Appl. Mech. Engrg. 148 (1997) 53–73. [3] P.M. Suquet, Local and global aspects in the mathematical theory of plasticity, in: A. Sawczuk, G. Bianchi (Eds.), Plasticity Today: Modelling, Methods and Applications, Elsevier Applied Science Publishers, London, 1985, pp. 279–310. [4] J.M. Guedes, N. Kikuchi, Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods, Comput. Methods Appl. Mech. Engrg. 83 (1990) 143–198. [5] K. Terada, N. Kikuchi, Nonlinear homogenization method for practical applications, in: S. Ghosh, M. Ostoja-Starzewski (Eds.), Computational Methods in Micromechanics, vol. AMD-212/MD-62, ASME, 1995, pp. 1–16. [6] S. Ghosh, K. Lee, S. Moorthy, Two scale analysis of heterogeneous elastic–plastic materials with asymptotic homogenisation and Voronoi cell finite element model, Comput. Methods Appl. Mech. Engrg. 132 (1996) 63–116. [7] R.J.M. Smit, W.A.M. Brekelmans, H.E.H. Meijer, Prediction of the mechanical behaviour of non-linear heterogeneous systems by multi-level finite element modeling, Comput. Methods Appl. Mech. Engrg. 155 (1998) 181–192. [8] C. Miehe, J. Schro¨der, J. Schotte, Computational homogenization analysis in finite plasticity. Simulation of texture development in polycrystalline materials, Comput. Methods Appl. Mech. Engrg. 171 (1999) 387–418. [9] C. Miehe, A. Koch, Computational micro-to-macro transition of discretized microstructures undergoing small strain, Arch. Appl. Mech. 72 (2002) 300–317. [10] J.C. Michel, H. Moulinec, P. Suquet, Effective properties of composite materials with periodic microstructure: a computational approach, Comput. Methods Appl. Mech. Engrg. 172 (1999) 109–143. [11] F. Feyel, J.-L. Chaboche, FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fiber SiC/Ti composite materials, Comput. Methods Appl. Mech. Engrg. 183 (2000) 309–330. [12] K. Terada, N. Kikuchi, A class of general algorithms for multi-scale analysis of heterogeneous media, Comput. Methods Appl. Mech. Engrg. 190 (2001) 5427–5464. [13] S. Ghosh, K. Lee, P. Raghavan, A multi-level computational model for multi-scale damage analysis in composite and porous materials, Int. J. Solids Struct. 38 (2001) 2335–2385. [14] V. Kouznetsova, W.A.M. Brekelmans, F.P.T. Baaijens, An approach to micro–macro modeling of heterogeneous materials, Comput. Mech. 27 (2001) 37–48. [15] N.A. Fleck, J.W. Hutchinson, Strain gradient plasticity, Adv. Appl. Mech. 33 (1997) 295–361. [16] E.C. Aifantis, On the microstructural origin of certain inelastic models, Trans. ASME J. Engrg. Mater. Tech. 106 (1984) 326–330. [17] R. de Borst, H.-B. Mu¨hlhaus, Gradient-dependent plasticity: formulation and algorithmic aspects, Int. J. Numer. Meth. Engrg. 35 (1992) 521–539. [18] N.A. Fleck, J.W. Hutchinson, A reformulation of strain gradient plasticity, J. Mech. Phys. Solids 49 (2001) 2245–2271.

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

5549

[19] G. Pijaudier-Cabot, Z.P. Bazˇant, Nonlocal damage theory, J. Eng. Mech. 113 (1987) 1512–1533. [20] R.H.J. Peerlings, R. de Borst, W.A.M. Brekelmans, J.H.P. de Vree, Gradient-enhanced damage for quasi-brittle materials, Int. J. Numer. Meth. Engrg. 39 (1996) 3391–3403. [21] R.A.B. Engelen, M.G.D. Geers, F.P.T. Baaijens, Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour, Int. J. Plasticity 19 (2003) 403–433. [22] M.G.D. Geers, R.L.J.M. Ubachs, R.A.B. Engelen, Strongly nonlocal gradient-enhanced finite strain elastoplasticity, Int. J. Numer. Meth. Engrg. 56 (2003) 2039–2068. [23] V.P. Smyshlyaev, N.A. Fleck, Bounds and estimates for linear composites with strain gradient effects, J. Mech. Phys. Solids 42 (12) (1994) 1851–1882. [24] W.J. Drugan, J.R. Willis, A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites, J. Mech. Phys. Solids 44 (4) (1996) 497–524. [25] N. Triantafyllidis, S. Bardenhagen, The influence of scale size on the stability of periodic solids and the role of associated higher order gradient continuum models, J. Mech. Phys. Solids 44 (11) (1996) 1891–1928. [26] H.T. Zhu, H.M. Zbib, E.C. Aifantis, Strain gradients and continuum modeling of size effect in metal matrix composites, Acta Mech. 121 (1997) 165–176. [27] O. van der Sluis, P.H.J. Vosbeek, P.J.G. Schreurs, H.E.H. Meijer, Homogenization of heterogeneous polymers, Int. J. Solids Struct. 36 (1999) 3193–3214. [28] M. Ostoja-Starzewski, S.D. Boccara, I. Jasiuk, Couple-stress moduli and characteristic length of a two-phase composite, Mech. Res. Comm. 26 (4) (1999) 387–396. [29] V.P. Smyshlyaev, K.D. Cherednichenko, On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media, J. Mech. Phys. Solids 48 (2000) 1325–1357. [30] R.H.J. Peerlings, N.A. Fleck, Numerical analysis of strain gradient effects in periodic media, J. Phys. IV 11 (2001) 153–160. [31] S. Forest, F. Pradel, K. Sab, Asymptotic analysis of heterogeneous Cosserat media, Int. J. Solids Struct. 38 (2001) 4585–4608. [32] M.G.D. Geers, V. Kouznetsova, W.A.M. Brekelmans, Gradient-enhanced computational homogenization for the micro–macro scale transition, J. Phys. IV 11 (2001) 145–152. [33] V. Kouznetsova, M.G.D. Geers, W.A.M. Brekelmans, Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme, Int. J. Numer. Meth. Engrg. 54 (2002) 1235–1260. [34] V. Kouznetsova, Computational homogenization for the multi-scale analysis of multi-phase materials, Ph.D. thesis, Eindhoven University of Technology, Eindhoven, The Netherlands, 2002. [35] R. Hill, Elastic properties of reinforced solids: some theoretical principles, J. Mech. Phys. Solids 11 (1963) 357–372. [36] S. Nemat-Nasser, Averaging theorems in finite deformation plasticity, Mech. Mater. 31 (1999) 493–523. [37] R.D. Cook, D.S. Malkus, M.E. Plesha, Concepts and Applications of Finite Element Analysis, Wiley, Chichester, 1989. [38] C. Miehe, Numerical computation of algorithmic (consistent) tangent moduli in large-strain computational inelasticity, Comput. Methods Appl. Mech. Engrg. 134 (1996) 223–240. [39] C. Miehe, Computational micro-to-macro transitions for discretized micro-structures of heterogeneous materials at finite strains based on the minimization of averaged incremental energy, Comput. Methods Appl. Mech. Engrg. 192 (2003) 559–591. [40] R.A. Toupin, Elastic materials with couple-stress, Arch. Ration. Mech. Anal. 11 (1962) 385–414. [41] R.A. Toupin, Theories of elasticity with couple-stress, Arch. Ration. Mech. Anal. 17 (1964) 85–112. [42] R.D. Mindlin, Micro-structure in linear elasticity, Arch. Ration. Mech. Anal. 16 (1964) 51–78. [43] R.D. Mindlin, N.N. Eshel, On first strain-gradient theories in linear elasticity, Int. J. Solids Struct. 4 (1968) 109–124. [44] H. Gao, Y. Huang, W.D. Nix, J.W. Hutchinson, Mechanism-based strain gradient plasticity––I. Theory, J. Mech. Phys. Solids 47 (1999) 1239–1263. [45] Z.C. Xia, J.W. Hutchinson, Crack tip fields in strain gradient plasticity, J. Mech. Phys. Solids 44 (10) (1996) 1621–1648. [46] M.R. Begley, J.W. Hutchinson, The mechanics of size-dependent indentation, J. Mech. Phys. Solids 46 (10) (1998) 2049–2068. [47] A. Zervos, P. Papanastasiou, I. Vardoulakis, A finite element displacement formulation for gradient elastoplasticity, Int. J. Numer. Meth. Engrg. 50 (2001) 1369–1388. [48] L.R. Herrmann, Mixed finite elements for couple-stress analysis, in: S.N. Atluri, R.H. Gallagher, O.C. Zienkiewicz (Eds.), Hybrid and Mixed Finite Element Methods, John Wiley & Sons Ltd., 1983, pp. 1–17. [49] J.Y. Shu, W.E. King, N.A. Fleck, Finite elements for materials with strain gradient effects, Int. J. Numer. Meth. Engrg. 44 (1999) 373–391. [50] E. Amanatidou, A. Aravas, Mixed finite element formulations of strain-gradient elasticity problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 1723–1751. [51] T. Matsushima, R. Chambon, D. Caillerie, Large strain finite element analysis of a local second gradient model: application to localization, Int. J. Numer. Meth. Engrg. 54 (2002) 499–521. [52] J.Y. Shu, C.Y. Barlow, Strain gradient effects on microscopic strain field in a metal matrix composite, Int. J. Plasticity 16 (2000) 563–591.

5550

V.G. Kouznetsova et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 5525–5550

[53] J.Y. Shu, N.A. Fleck, E. van der Giessen, N. Needleman, Boundary layers in constrained plastic flow: comparison of nonlocal and discrete dislocation plasticity, J. Mech. Phys. Solids 49 (2001) 1361–1395. [54] E. Bittencourt, A. Needleman, M.E. Gurtin, E. van der Giessen, A comparison of nonlocal continuum and discrete dislocation plasticity predictions, J. Mech. Phys. Solids 51 (2003) 281–310. [55] M.G.D. Geers, V.G. Kouznetsova, W.A.M. Brekelmans, Multi-scale first-order and second-order computational homogenization of microstructures towards continua, Int. J. Multiscale Comput. Engrg. in press.