Objective evaluation of linearization procedures in nonlinear homogenization: A methodology and some implications on the accuracy of micromechanical schemes

A systematic methodology for an accurate evaluation of various existing linearization procedures sustaining mean ﬁelds theories for nonlinear composites is proposed and applied to recent homogenization methods. It relies on the analysis of a periodic composite for which an exact resolution of both the original nonlinear homogenization problem and the linear homogenization problems associated with the chosen linear comparison composite (LCC) with an identical microstructure is possible. The eﬀects of the sole linearization scheme can then be evaluated without ambiguity. This methodology is applied to three diﬀerent two-phase materials in which the constitutive behavior of at least one constituent is nonlinear elastic (or viscoplastic): a reinforced composite, a material in which both phases are nonlinear and a porous material. Comparisons performed on these three materials between the considered homogenization schemes and the reference solution bear out the relevance and the performances of the modiﬁed second-order procedure introduced by Ponte Castan˜eda in terms of prediction of the eﬀective responses. However, under the assumption that the ﬁeld statistics (ﬁrst and second moments) are given by the local ﬁelds in the LCC, all the recent nonlinear homogenization procedures still fail to provide an accurate enough estimate of the strain statistics, especially for composites with high contrast.


Introduction
Most of the nonlinear homogenization procedures for heterogeneous materials implicitly or explicitly rely on two separates stages.The first one consists in approximating the actual nonlinear behavior of the individual constituents of the composite by linear constitutive equations, through a specific linearization procedure, so as to define a fictitious elastic or thermoelastic composite usually referred to as the linear comparison composite ''LCC'' (Ponte Castan ˜eda, 1991).In the second stage, the overall as well as the local responses of this fictitious LCC are approximated by a linear homogenization model appropriate for the microstructure of the LCC.If in each stage, the derivation ensures that the result is an upper bound for the effective properties, the obtained global estimate is a bound as well (Ponte Castan ˜eda, 1991) and its interpretation is clear.However in such a case, the predictions of the model might be far away from the exact response of the composite and thus be of limited practical interest.A less rigorous but more efficient linearization scheme might then be preferred and has to be selected among the numerous classical or more recently proposed formulations, such as Hill's incremental scheme (Hill, 1965), the classical (Hutchinson, 1976;Berveiller and Zaoui, 1979) or modified (Ponte Castan ˜eda, 1991;Suquet, 1995) secant approaches, the affine formulation (Masson and Zaoui, 1999;Masson et al., 2000) and its variants (Ponte Castan ˜eda, 2002;Chaboche and Kanoute ´, 2003), the second-order procedures (Ponte Castan ˜eda, 1996, 2002) or the Lahellec and Suquet procedure (Lahellec and Suquet, 2004).These models are based on more or less sophisticated derivations and generate various predictions, which deserve an objective and systematic comparison to each other and with exact solutions, both at the local and global levels to appreciate their respective merits and limitations.It is the purpose of the present study to provide a methodology for such a systematic and objective comparison, from which, in addition, guidelines for improved formulations might be deduced.
A first and usual evaluation of the various linearization procedures consists in comparing their overall predictions with the few available (mostly upper) bounds, as in some previous studies (Gilormini, 1996).However, such bounds are often not sufficiently sharp for a precise evaluation, except in some specific situations (Bornert and Ponte Castan ˜eda, 1998).For weakly inhomogeneous composites, the predictions may be compared with the expansions of the effective potential (Ponte Castan ˜eda and Suquet, 1995Suquet, , 1998) ) with respect to the contrast, at various orders.In more general situations, accurate evaluations of the various nonlinear homogenization procedures can be obtained by means of comparisons with full-field numerical solutions of the initial nonlinear problem.To this end, two different approaches may be considered.For the first one, numerical simulations are carried out on large windows of simulated microstructures (Moulinec and Suquet, 2003) supposed to depict the random microstructures addressed by the linear mean-field theories used for the homogenization of the LCC.Such an approach induces computational difficulties due to the large size of the numerical systems to work out, as well as issues relative to the representativeness of the generated microstructures, the appropriate averaging of the results and the choice of particular boundary conditions, both in the linear (Kanit et al., 2003) case and, even more critically, in the nonlinear one.This approach implies a large computational expense and is in practice often restricted to two-dimensional problems.Moreover, when three dimensional simulations are carried out, the current computer limitations require to restrict the number of comparisons and thus forbid to thoroughly explore all combinations of the parameters of the particular problem under consideration (Bilger et al., 2005).The second approach consists in performing comparisons between the nonlinear homogenization procedures and exact solutions on simple periodic microstructures.The main advantage of this approach lies in a reasonable computational expense as well as an easier control of the numerical accuracy of the results (Chaboche and Kanoute ´, 2003;Suquet, 1995).However, such an approach is questionable since the comparison between the reference numerical solution and predictions from a nonlinear mean field theory based on a linear homogenization model for a random LLC can be strongly altered by the fact that the compared procedures do not address the same microstructure.They are also distorted by the use of linear closedform estimates-which are not exact results-to evaluate the effective properties of the LCC derived by the linearization procedure of the mean field models.Such linear estimates may not be sufficiently accurate and can provide crude approximations for complex microstructures or when the mechanical contrast between constituents is high.Therefore, comparisons carried out in such a framework may lead to potentially ambiguous results because of the addition of two approximations, stemming from the linearization procedure and the linear homogenization scheme, which may be cumulative or compensative.Finally, comparisons of predictions of nonlinear homogenization schemes may of course also be compared to experimental data.However, such comparisons are very difficult since they combine approximations at various levels: the actual nonlinear constitutive law of the constituents, the unavoidable simplification of the microstructure used in the models with respect to the real one, the experimental scatters on the evaluation of microstructural parameters such as volume fractions and of the overall responses as well as the complexity of local stress or strain measurements (Bornert et al., 1994).For theses reasons, at the present stage, ''numerical experiments'' are strongly preferred for such comparisons, since their inputs are known exactly and their outputs can be controlled.The main objective of this work is to remove their aforementioned limitations in order to be left with an objective evaluation of the sole linearization methods.
To this end, a new methodology which relies on an exact treatment of both the nonlinear and the linear homogenization problems is developed.Various nonlinear homogenization formulations are compared, mainly with regard to their predictions in terms of overall responses.In addition, a preliminary evaluation of the local predictions is performed, restricted to the first and second moments of the local strain field, for the sake of brevity.At this stage, it is worth noting that, as in earlier works including very recent ones (e.g., Idiart et al., 2006), it is assumed that the local fields predicted by the nonlinear homogenization schemes are given by those in the LCC.This is consistent with the assumptions of classical formulations based on local fields averaging but can lead to some inconsistencies when overall responses are derived from evaluations of effective potentials (Masson et al., 2000;Lahellec and Suquet, 2004).As shown in a very recent paper by Idiart and Ponte Castan ˜eda (2006), not yet available at the time of submission of the present contribution, this limitation can be removed according to the general theory developed by these authors, which enables to extract in a consistent way the statistics of the local fields in the nonlinear composite from the knowledge of the effective potential of suitably perturbated composites.The corresponding correction terms, giving the discrepancy between the statistics on the LCC and those on the nonlinear composite, are not taken into account in the present study.Accordingly, the conclusions regarding the efficiency in terms of local predictions of these models should be considered as restricted to their simplified classical ''LCC-based local interpretation''.This will be discussed again with more details later on.
The structure of this paper is as follows.The methodology to evaluate the various linearization treatments sustaining nonlinear homogenization procedures is presented in Section 2. Its implementation is depicted in Section 3 where a complete description of the linear numerical elastic or thermoelastic homogenization scheme and its connections with the linearization formulations is given.The main assumptions of the various linearization schemes that will be compared are recalled at the end of this section.Namely, these models are the classical secant scheme (referred to as SEC), the variational procedure (VAR), the original affine formulation (AFF-ANI), its isotropic simplification (AFF-ISOT), the original (SOE-1) and improved (SOE-2) second-order procedures, and finally the Lahellec and Suquet (LS) formulation.The even more recent formulations such as those proposed by Idiart and Ponte Castan ˜eda (2005) and Idiart et al. (2006) have not been considered in the present work.On the other hand, Hill's original incremental procedure is not considered either since its limitations have already been established earlier (Gilormini, 1996); its isotropic simplification, recently reassessed (Doghri and Ouaar, 2003), is not evaluated either since it is very close to the classical secant formulation, both schemes coinciding exactly when the local constitutive laws are power-laws with same exponent (Hutchinson, 1976).Similarly, the popular viscoplastic self-consistent model (Molinari et al., 1987;Lebensohn and Tome ´, 1993), based on a tangent approximation of the nonlinear local constitutive law, has been omitted since it does not clearly separate the homogenization of a LCC from the linearization procedure and is restricted to random polycrystals.The presented methodology is carried out in Section 4 on three different two-phase materials: a nonlinear matrix reinforced by linear elastic particles (referred to in the sequel as the ''reinforced composite''), a composite in which both phases exhibit a nonlinear behavior (referred to as the ''two-phase material'') and a porous material.Based on the numerous comparisons performed both at the global and local scales, a comparative evaluation of the various linearization procedures is finally proposed and the conclusions are summarized in Section 5.
The tensor notation used herein is a fairly standard notation.The orders of the tensors are clear when taken in context.Products containing dots denote summation over repeated latin indices.For example, L : e = L ijkl e kl e i e j and E :: F = E ijkl F klij where (e i , i = 1, 3) is a time-independent orthogonal cartesian basis and denotes the tensorial product.Cylindrical coordinates will be used as well, with e r , e h , e z = e 3 being the unit vectors of the orthogonal cylindrical basis and u r , u h , u z the cylindrical coordinate components of a 3D vector u.

Principle
The methodology and the numerical tool which are proposed in this paper enable a systematic and accurate evaluation of the linearization procedures without suffering from the limitations mentioned in introduction.We address a problem where the exact nonlinear solution, regarded as the reference solution, may be computed with high accuracy at a low computational expense.Moreover, the homogenization of the LCC, which retains the microstructure of the nonlinear composite, is accurately computed by the same numerical tool as the one handling the nonlinear composite.Thereby, the difficulties related to the numerical approximations, to the change of microstructure as well as to the approximations induced by the linear estimates of the LCC are avoided.The effect of the sole linearization procedure may be assessed without any ambiguity.
The chosen microstructure is periodic and classical finite element techniques associated with periodic homogenization are used to solve both the initial nonlinear problem and the linear problem associated with the LCC.This methodology can be adopted for any type of local nonlinear constitutive relations, any number of phases and any microstructure, since large unit cells, able to mimic a random microstructure, may in principle be used.However, in order to restrict the computational expense, only small unit cells which describe simple periodic composites will be used as a first application.Because of this practical restriction, the proposed methodology cannot evaluate exhaustively all situations.It however allows to handle large classes of problems, since various choices of local constitutive laws, contrasts between constituents, volume fractions, local geometries or loadings are possible.It thereby extends considerably the earlier comparisons of this type and, most importantly, provides unambiguous comparisons.In addition, since the linearization procedures are often supposed to apply to any type of microstructure, they should also work for periodic ones and can therefore be tested in this case.Let us also stress the fact that even if the effective properties and the heterogeneity of the local fields are quantitatively different in a random microstructure and a periodic one with similar constituents, they are qualitatively close so that comparisons performed between models for periodic microstructures should be sufficiently representative of what might be observed on random ones.Note finally that a similar approach has independently been used by Lahellec and Suquet (2004) to evaluate the specific model proposed by these authors, as well as by Moulinec and Suquet (2004) for a comparison between the classical and the modified secant linearization procedures.In both these contributions, only a limited number of models have been compared and comparisons were restricted to incompressible two-dimensional microstructures made of fiber reinforced matrices loaded in their transverse plane, the response of which is characterized by a one-dimensional relation.In the present work, three dimensional microstructures with a wider range of phase contrasts are considered, including porous materials exhibiting a two-dimensional response.An additional fundamental difference between the proposed methodologies can be noticed: in these earlier works, the numerical simulations used for the evaluations of both the exact nonlinear behavior and the linear properties of the LCC aimed at reproducing complex random microstructures made of randomly distributed fibers.Expensive computational tools had to be used and statistical averages over results obtained with various realizations of complex unit cell had to be taken, giving rise to aforementioned corresponding questions relative to numerical and statistical convergence of the predictions.In the present proposed methodology, the explicit choice of simple periodic microstructures intrinsically avoids such difficulties and allows one to explore more exhaustively the performances of many formulations, under various conditions, at a very low numerical cost.This is at the price of the anisotropy of the properties of both the nonlinear medium and the LCC, which can however easily be dealt with.

The nonlinear periodic composite
As a first illustration of the proposed methodology, we consider composites with isotropic nonlinear elastic constituents, the behavior of which is governed by a single potential or strain energy function w(e), according to r ¼ f ðeÞ ¼ ow oe ðeÞ where r and e are the local stress and strain fields, respectively, in the framework of small strains.Note that this case is similar to the rate problem for a nonlinear viscoplastic composite.For the sake of simplicity, the studied material reduces to a periodic two-phase composite, made of aligned spherical inclusions (particles or pores) embedded in a matrix.The inclusions are distributed on a hexagonal network in the transverse plane and are aligned along the third direction, so that a cylinder with a hexagonal basis and a single spherical inclusion can be used as unit cell.All the loadings considered in this study are axially symmetric along the third direction.In a first and classical approximation (Michel et al., 2001), this hexagonal cell can be replaced by a cylindrical cell of height 2H with a circular basis of diameter 2R.Because of the invariance with respect to any rotation along the third direction of the geometry of the cell, of the constitutive relations of the phases and of the applied load, the 3D initial problem reduces to a 2D axially symmetric problem which can be solved at a very low numerical cost.
For numerical applications, the constitutive law r = f r (e) followed by phase r is a Ramberg-Osgood type relation with or without a threshold, an initial isotropic elastic behavior and a nonzero compressible part, which reads with e m ¼ r m 3k r and e eq ¼ g r ðr eq Þ ¼ In these expressions, k r ; l r e ; r r 0 and r r y are the elastic bulk modulus, the elastic shear modulus, the flow stress and the threshold stress of the phase r, respectively.The parameter m = 1/n is the work-hardening parameter (n is the nonlinearity exponent) satisfying 0 6 m 6 1, 0 is a reference strain and P(a) is the positive part of a.
The von Mises equivalents of the strain and the stress tensors are defined as usually by e eq ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 3 e : K : e q and r eq ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 3 2 r : K : r q , while r m ¼ 1 3 i : r and e m ¼ 1 3 i : e are their respective spherical parts.The tensors i and I are second and fourth-order identity tensors, and J ¼ 1 3 i i and K = I À J are the usual projectors on the subspaces of purely spherical or deviatoric second-order tensors.The nonlinearity exponent n is the same in both phases, even if more general situations could be considered.Unlike usually assumed, l r e and k r are finite in order to avoid practical numerical difficulties related to infinite tangent and secant moduli at e = 0 and to material incompressibility.We assume a perfect bonding between the inclusion and the matrix so that the displacement and surface traction are continuous at the interface.

Overall and local response under axial load
In the ensuing calculations, the unit cell is submitted to an axisymmetric macroscopic deformation where ê ¼ ðK : eÞ=e eq ¼ e 3 e 3 À 1 2 ðe 1 e 1 þ e 2 e 2 Þ.In the considered small strain formalism, e is the average of the local strain, say e ¼ hei, where the notation hCi ¼ 1 jV j

R
V CðxÞ dx denotes the mean value of C over the cylindrical unit cell V, whose measure is jVj = 2pR 2 H. Likewise and for a later use, hCi r ¼ 1 jV r j R V r CðxÞ dx denotes the average of C over phase r, which occupies the domain V r in the unit cell.Note that in the sequel e m and e eq are always positive quantities.The ratio s e ¼ e m =e eq characterizes the macroscopic strain triaxiality ratio.
The local stress r and strain e fields in the unit cell are solution of the following set of equations (Michel et al., 2001) uðxÞ where v r (x) is the characteristic function of phase r and where the general notations # and À# mean that the fluctuating part of the displacement vector u * and the surface traction r AE n (n being the outer unit normal) are periodic and anti periodic on the cell boundary oV, respectively.The superscripts m and p denote the matrix and the particles, respectively.The macroscopic stress r, given by r ¼ hri, is axially symmetric because of the rotational invariance of these equations along the axis e 3 .It reads r ¼ r m i þ 2 3 ir eq ê where r m and r eq are the overall hydrostatic and von Mises equivalent stresses, given by r m ¼ r 33 þ 2r 11 3 and r eq ¼ jr 33 À r 11 j; ð5Þ and where i is equal to +1 or À1.In general, the overall response of the nonlinear composite is thus characterized by the two-dimensional relation If the local constitutive relations were positively homogeneous functions of the same degree m, the local strain, solution of the system of equations ( 4), would be a positively homogeneous function of degree one of the imposed strain e, the local stress a positively homogeneous function of degree m and so would be the overall stress (Ponte Castan ˜eda and Suquet, 1998).The general relation (6) could then be cast either in the form or Since the constitutive phases are compressible and their elastic shear moduli are finite, the local constitutive relations never exhibit exactly such homogeneity relations, even when the threshold stresses r r y vanish.However, such homogeneity relation are asymptotically true for large strains, when the linear parts of the deformation ( rm 3k r and req 3l r e ) tend to be negligible with respect to the nonlinear one (e 0 ð req r r 0 Þ n ).For a reinforced composite or a two-phase material, the macroscopic behavior may be considered as almost incompressible, since the bulk moduli k r will be chosen sufficiently large.As a consequence, the sole overall quantity of interest is the deviatoric response r eq ¼ f ðe eq Þ, which is characterized by the single scalar r0 when the threshold stresses vanish.On the other hand, for the porous case, the effective response is compressible, even if the matrix is almost incompressible, and is characterized by the two-dimensional relation (6), which may asymptotically restrict to the pair of functions r0 ðs e Þ; r0 ðs e Þ as illustrated in Section 4. The comparison between the exact nonlinear solution and the estimates provided by various nonlinear homogenization schemes can be performed on these macroscopic quantities, seen as functions of the parameters of the nonlinear problem.In Section 4, the dependance of the overall response with respect to the volume fraction f p of the inclusions and the work-hardening parameter m will in particular be thoroughly investigated.
In addition to such macroscopic quantities, the comparisons can also be performed on predictions of local fields.Of course, one should not expect a same level of accuracy of the evaluated homogenization theories in terms of local field predictions as in terms of effective properties.However, as such models might be used for micromechanical applications, in which the local field predictions might serve to predict the activation of local physical mechanisms, such as damage or other microstructure evolutions, it is useful to assess the predictive capabilities of these models even at the local scale, mostly in terms of statistics.More precisely, one can investigate the so-called inter-phase heterogeneity of the stress and strain fields, characterized by the first-order moments e r ¼ hei r and r r ¼ hri r , as well as the intra-phase heterogeneities, given by the second-order moments he ei r and hr ri r , or the variances C r e ¼ hðe À e r Þ ðe À e r Þi r ¼ he ei r À hei r hei r and C r r ¼ hðr À r r Þ ðr À r r Þi r ¼ hr ri r À hri r hri r , which quantify the fluctuations of the local fields in the phases with respect to their average.Because of the rotational invariance of the nonlinear homogenization problem, the first-order moments can be decomposed into their spherical and deviatoric parts, the latter being proportional to e ˆ; they are thus characterized by two scalars.Similarly, the second-order moments and variances are transversely isotropic fourth-order symmetric tensors, which are in general characterized by 5 scalars (see for instance Walpole (1981)).Among these, the following quantities are of special interest since, as will be recalled later, several linearization schemes are based on them where E ¼ 2 3 ê ê and F = K À E are the classical fourth-order projectors (Ponte Castan ˜eda, 1996) relative to the directions parallel and perpendicular to the overall load.A thorough evaluation of the linearization schemes at the local scale is out of the scope of the present paper which is focussed on the comparisons at the macroscopic scale.A few results relative to these quantities will however be presented at the end of Section 4.More detailed comparisons could also be performed in terms of statistical distribution functions (as in Bornert et al., 1994) or even in terms of the spatial distribution of local fields in the form of strain maps on the unit cell.The reader is invited to refer to Rekik (2006) for such additional results.It is noted that such comparisons are possible within the proposed methodology, since the local fields can be computed exactly both in the nonlinear composite and in the LCC; in more traditional approaches in which the LCC is homogenized by means of a mean-field linear model, only the overall properties and the per phase first and secondorder moments of the local fields are available.At this stage it is useful to recall that the local predictions of the homogenization schemes are here assumed to be given by the local fields in the LCC.This is consistent with most evaluated formulations; exceptions will be indicated later.

Implementation
We address here more technical issues relative to the computation of the reference exact nonlinear solution (Section 3.1) and the homogenization of the linear elastic (Section 3.2.1)or thermoelastic (Section 3.2.2) LCC.The main assumptions of the evaluated linearization schemes are recalled in Section 3.3.and Section 3.4 is a short discussion on the accuracy of the numerical results.Let us before recall that, because of the rotational invariance along the third axis of all the data of the nonlinear problem (cell geometry, isotropic nonlinear local constitutive laws, macroscopic imposed strain), all quantities derived from these data will exhibit the same invariance.More specifically, • Second-order tensors, such as the overall stress and per phase averages of stress and strain, admit the decomposition where c m and c d are (positive or negative) scalars.If c is a strain, jc d j = c eq while jc d j ¼ 2 3 c eq if c is a stress.• Symmetric second-order tensor fields (such as local stress and strain fields) are axisymmetric and so are first-order tensor fields like displacements • Fourth-order tensors, such as the aforementioned covariance tensors, are transversely isotropic and admit the decomposition The adopted notation is that of Bornert and Suquet (2001a), which is derived from the decomposition of transversely isotropic tensors introduced by Walpole (1981).The definition of the tensors E L , J T , F 0 , K T and K L , is recalled in Appendix A.
Such decompositions, which are quite obvious for the nonlinear problem, hold true for the quantities associated with the LCC as well.Indeed the linearization procedures of all the considered schemes do not introduce any additional preferential direction so that the LCC retains the invariance properties of the nonlinear problem.In particular, the tensors of moduli of the phases admit decomposition (13), with c = c 0 .In addition, as shown later, the determination of the relevant effective properties of the LCC can be performed with the application of overall axisymmetric strain tensors similar to Eq. (3).

Reference solution
The reference nonlinear solution is obtained by solving system (4) with standard finite element techniques.Because of the rotational invariance, the displacement field admits decomposition ( 12) and the stress and strain fields can be written as in Eq. ( 11), so that two-dimensional axisymmetric elements can be used.Because of the additional symmetry with respect to the transverse plane, only the domain (r, z) has to be meshed, the origin of the coordinates being centered on the inclusion.The symmetry and periodicity conditions reduce then to the simple set of equations (see Tvergaard (1982); Koplik and Needleman (1988) for details) u r ðR; zÞ ¼ Re 11 and u r ð0; zÞ ¼ 0; 0 6 z 6 H ; with e 11 ¼ e m À e eq =2 and e 33 ¼ e m þ e eq : ð15Þ The FE code CAST3M has been used for the determination of the local fields.An appropriate iterative procedure is used to address the nonlinearity of the system.After convergence at an imposed macroscopic strain, the global stress components are obtained by averaging the other components being zero, as already noted.The overall mean and deviatoric stresses are then obtained from Eq. ( 5).Similar relations hold true for the computation of the per phase average strains and stresses, with an integration restricted to the corresponding phase and jVj replaced by its volume.The scalar second-order moment e r eq is obtained by averaging the square of the local equivalent strain The parallel and perpendicular components of the (tensorial) second-order moments are obtained from To get all the desired quantities to characterize the inter-and intra-phase heterogeneity of the strain field, one has thus to integrate over each phase domain the four scalars e zz , (e rr + e hh ), e 2 eq and ðe : êÞ 2 .This can be performed with standard operators of the FE code CAST3M.Similar relations hold for the stress field.

Linear homogenization
The equations governing the global and local response of the LCC are similar to those relative to the nonlinear problem (system 4), the only difference being that the local constitutive relation is now linear, taking the form The fourth-order tensor L r is the homogeneous tensor of moduli and the second-order tensor s r is the uniform polarization of phase r.Their values depend on the linearization scheme under consideration, as detailed in the next section.However, because of the global symmetry of the problem, these quantities are invariant with respect to any rotation along the symmetry axis.The tensor L r then reads The following results are given under the assumption that L r takes this general form.Indeed, as recalled later, the actual tensors defined by the evaluated linearization schemes admit the decomposition (Ponte Castan ˜eda, 1996) where E r ¼ 2 3 êr êr , F r = K À E r and êr ¼ ðK : e r Þ=e eq .Since e r is axially symmetric, we have K : e r d ¼ e r eq ê and therefore The tensor L r is therefore transversely isotropic and can be written in the more general form ( 23) with the following relations between the constants (k r , k r , l r ) and For some linearization schemes, L r is isotropic.For these, the previous relations still hold, with the simplification k r = l r .The polarization s r vanishes for all linearization schemes that define a linear elastic LCC.For all other schemes, s r is axially symmetric as a consequence of the general symmetry properties of the considered problem.In addition, it turns out that the polarization s r is trace free, because the considered local nonlinear constitutive equations are linear on their spherical part.The tensor s r then takes the form where js r d j is the von Mises equivalent stress of s r .The tensors L r and s r are ''inputs'' defining the LCC.The other inputs are its microgeometry, which, consistently with the presented methodology, is exactly the same as the microstructure of the nonlinear composite and can be represented by exactly the same cylindrical unit cell under periodic boundary conditions.The third input is the overall load applied to the LCC.A full homogenization of the LCC would require to investigate the global and local responses of the LCC under any load.However, all required ''outputs'' turn out to be obtained with axially symmetric overall strains.As a consequence, the same 2D axisymmetric cell, with the boundary conditions ( 14), can be used again.More precisely, the expected ''outputs'' are the overall stress induced in the LCC by the overall strain imposed on the nonlinear composite, the first-order moment of the strain in each phase of the LCC, which, as previously, admits the decomposition (10), and the traces (9) of the second-order moments or the variances of the strain fields.Similar stress-like quantities might also be useful.Note that these local quantities are of interest, first for comparison purposes as was already the case for the nonlinear homogenization problem, but also in order to define the constitutive relations of the LCC (i.e., L r and s r ), according to the principles of the linearization scheme under consideration, in an iterative and (hopefully) converging procedure.
In principle, since linear elastic constitutive relations are particular cases of nonlinear ones, all the relations given in Section 3.1 apply to the LCC.Because of the linearity, additional properties apply.They are detailed in the following two sections, devoted first to the case s r = 0 and then to the general thermoelastic situation.It is in particular shown how the simpler linear elastic results can be efficiently extended to the second case, without additional computations.

Elastic case
The linearity of the local problem ensures the existence of the localization tensor field A(x) and its per phase averages The tensors A r are fourth-order tensors, not necessarily symmetric, which depend on the microstructure and on the local properties, but not on the overall load.They are transversely isotropic and read Note that the coefficient of the localization tensor of one phase are deduced from those of the other by the relation where the f r = hv r i denote the volume fractions of the phases.The overall stress is obtained from The tensor of effective moduli L ¼ hL : Ai admits the decomposition (13) with coefficients deduced from Eqs. ( 33), ( 31) and ( 23) From Eqs. ( 31), ( 32) and ( 34), it is clear that the derivation of the overall stress and the per phase first-order moments of the strain induced by an axisymmetric overall deformation only requires the determination of the four scalars ða A r ; b A r ; c A r ; c 0 A r Þ relative to one of the constituents.This can easily be achieved by solving the localization problem with the FE tool for two independent axisymmetric loads.For instance e ¼ e 33 e 3 e 3 gives a A r ¼ he zz i r and c A r ¼ he rr þ e hh i r = ffiffi ffi 2 p while e ¼ e 1 e 1 þ e 2 e 2 leads to b A r ¼ he rr þ e hh i r =2 and . Alternatively, one may also choose the elementary loads e ¼ i and e ¼ ê.In the sequel, the corresponding local fields will be called e e¼i ðxÞ ¼ AðxÞ : i and e e¼ê ðxÞ ¼ AðxÞ : ê, respectively.The relations giving the constants ða A r ; b A r ; c A r ; c 0 A r Þ from the averages of e e¼i and e e¼ê over phase r are easily obtained but are not given here for brevity.Note that the computation of these quantities requires the construction of one stiffness matrix, its inversion, two matrix products and four spatial integrations of scalar quantities.
To obtain the traces (9), additional computations are required as described hereafter.The local field under the general axisymmetric load ðe m ; e eq Þ is given by superposition as eðxÞ ¼ e m e e¼i ðxÞ þ e eq e e¼ê ðxÞ: ð35Þ Its second-order moment in phase r is given by he ei r ¼ e 2 m he e¼i e e¼i i r þ e 2 eq he e¼ê e e¼ê i r þ e m e eq ðhe e¼ê e e¼i i r þ he e¼i e e¼ê i r Þ; ð36Þ so that e r eq reads e m e eq he e¼ê : K : Similarly, we have from which the quantities e r k and e r ?can be deduced as was done at the end of Section 3.1.This shows that six additional integrations are then required, among which two involve quantities that couple the local fields induced by two different loads.Again this can be performed by the object-oriented operators of CAST3M.

Thermoelastic extension
The previous relations can be extended to situations with s r 5 0. The local strain field and their averages read now where A(x) coincides with the localization tensor field of the elastic LCC introduced in Section 3.2.1.The second-order tensor field a(x) is the local strain under vanishing overall strain.Since we are dealing with twophase composites, the so-called Levin's relations apply (Levin, 1967) and allow to compute a(x) from A(x) In addition, the overall stress-strain relation reads We refer to Willis (1983) for details on such derivations.The computation of L : e can be performed exactly as in the purely elastic case and requires only the four constants ða A r ; b A r ; c A r ; c 0 A r Þ relative to one of the constituents.Since s r is proportional to e ˆ, the computation of s r : A r requires the same four constants.The overall stress induced by a given macroscopic strain is thus obtained from the same quantities as those required in the elastic case.
For the per phase first-order moment of the strain field, we have Since Ds is proportional to e ˆand L r takes the form L r = 3k r J + 2k r E + 2l r F, the product (DL) À1 Ds will also be proportional to ê, so that B and C read B = B eq e ˆ(B eq being positive or negative) and C = C m i + C eq e ˆwith C m ¼ e m and C eq ¼ e eq þ B eq .The per phase first-order moments of the strain are then obtained from Eq. ( 45), which requires again the four constants ða A r ; b A r ; c A r ; c 0 A r Þ obtained as in the elastic case.Traces of the second-order moments are obtained with relations similar to Eqs. ( 37) and ( 39) and require the same additional integrations.More specifically The computation of the traces along K, E and F requires the same six quadratic integrals already required in the elastic case: hðe e¼i eq Þ 2 i r , hðe e¼ê eq Þ 2 i r , hðe e¼i : êÞ 2 i r , hðe e¼ê : êÞ 2 i r , he e¼i : K : e e¼ê i r and he e¼i : E : e e¼ê i r , as well as the two averages he e¼i : êi r and he e¼ê : êi r , which derive from the four constants ða A r ; b A r ; c A r ; c 0 A r Þ.This shows that, even if the final expressions are slightly more complex, the numerical complexity of a linearization scheme based on a thermoelastic LCC is identical to that of a scheme using a linear LCC.

Nonlinear formulations
In this section, we describe the different nonlinear formulations which provide the elastic and thermoelastic LCC and spell out how to solve the whole nonlinear system consisting of both the linear (thermo)elastic homogenization problem and the linearization procedure.
The various linearization procedures can be split into two categories.The first consists of the nonlinear extensions which define a LCC where the phase behavior is described by a local (thermo)elastic stress-strain relation r(x) = L r (e r ): e(x) (+s r (e r )).For some linearization procedures such as the secant or affine approaches, the reference strain e r is set to the average of the strain field over the phase r, namely e r ¼ e r eq .For other formulations such as the secant modified extension-i.e., the variational procedure -, the strain field fluctuations are taken into account in such a way that the proposed prescription for the reference strain becomes e r ¼ e r eq .For those linearization procedures which we will refer to as the ''stress-strain approaches'', the macroscopic stress r is defined as the mean value of the local stress field, say r ¼ hri, computed according to Eq. ( 44).
For the second class of linearization procedures which we will refer to as the ''potential based approaches'', the effective stress is defined according to Hill's theorem as the derivative of the effective potential with respect to the effective strain r ¼ ow oe ðeÞ.The mathematical form of the overall potential depends on the chosen linearization procedures.Generally, the main idea is to define a thermoelastic LCC whose effective potential can be used to estimate the effective potential of the nonlinear composite.The thermoelastic strain potential of the thermoelastic LCC reads where the thermoelastic strain potentials w r T ðeÞ are second-order Taylor-type expansions (or third-order expansion for the LS approach) defined as In Eq. ( 48), w r (e) is the nonlinear strain-energy function of the phase r of the nonlinear composite.The reference strains e r and the moduli L r are constant inside each phase r.Those tensors are chosen such that they generate the best possible estimates of the nonlinear effective potential w from known estimates of the thermoelastic effective potential wT .In general, the optimal values of e r and L r are derived from stationarity conditions.For the initial second-order method and the LS procedure, L r ¼ ðo 2 w=oe 2 Þðe r Þ.For the second-order modified procedure, the expression of the tensors L r is more complicated (Eq.( 52)).As mentioned above, the macroscopic constitutive behavior of the potential based approaches is evaluated by the derivative of the effective potential with respect to the effective strain.Since both the effective stress and strain are axially symmetric, the effective behavior can be shown to be defined, as in the purely isotropic case, by two scalar equations r m ¼ ow oe m ðeÞ and r eq ¼ ow oe eq ðeÞ; ð49Þ where the derivatives are found numerically.

The nonlinear problem
Before presenting in more details the various nonlinear formulations, we first focus our attention on the whole nonlinear system to be solved.For all the linearization procedures but the modified second-order approach (Ponte Castan ˜eda, 2002), the nonlinear problem can be written in the following form nonlinear relations where L r (e) and s r (e) are known functions whose exact expressions depend on the chosen linearization procedure.To solve this system of equations, a fixed-point iterative procedure is used.An initial value denoted p r 0 is ascribed to the reference strain e r for each phase.From Eq. (50e), we obtain the starting values of the moduli tensors L r and the polarization s r .Then, the numerical procedures described in Sections 3.2.1 and 3.2.2 are applied to provide a first evaluation of the first e r eq and second-order e r eq moments of the strain field in the LCC.According to the chosen linearization procedure, Eq. (50f) provides a new evaluation of the reference strain denoted p r 1 .This iterative procedure goes on until the discrepancy p r 1 À p r 0 meets the convergence criterion.Generally, three or four iterative steps are sufficient to reach convergence.For the modified second-order procedure, the nonlinear problem is solved by the same numerical scheme except that the moduli L r , as shown in Eq. ( 52), depend on two reference strains e r and êr .

Stress-strain approaches
Classical and modified secant formulations.For both formulations, the LCC is elastic (s r = 0).The tensor L r is defined as the isotropic tensor of secant moduli evaluated at the reference strain e r , namely L r ¼ 3k r J þ 2l r sct ðe r ÞK where l r sct ðe eq Þ ¼ r eq ðe eq Þ=3e eq is the secant shear modulus.For the Ramberg-Osgood constitutive law, e eq is related to r eq through Eq. (2b) and r eq = f r (e eq ) is derived by means of a numerical inversion procedure.For the classical secant formulation SEC, the reference strain is prescribed to be the per phase average of the strain field in the LCC-that is e r ¼ e r eq -whereas for the modified secant extension VAR, which coincides with the variational approach of Ponte Castan ˜eda (1991) (Suquet, 1995), the reference strain becomes e r ¼ e r eq .It is worth recalling that the variational procedure provides an upper bound for the effective energy.
Affine formulations.For the original affine approach referred to as AFF-ANI (Masson et al., 2000), the linear constitutive behavior of the individual constituents of the LCC follows a thermoelastic constitutive law where the elastic moduli and the polarization tensors are defined by L r ¼ L r tgt ðe r Þ ¼ o 2 w r oe 2 ðe r Þ and s r ¼ ow r oe ðe r Þ À L r : e r , respectively.Note that the thermoelastic law of each individual constituent of the LCC is tangent to the constitutive law of each phase of the real nonlinear composite r ¼ ow r oe ðeÞ at e ¼ e r .Unlike the secant methods, the tensor L r is transversely isotropic and can be written as L r ¼ 3k r J þ 2l r sct ðe r eq ÞF þ 2l r tgt ðe r eq ÞE where l r tgt ðe eq Þ ¼ dreq 3deeq ðe eq Þ is the tangent shear modulus.Due to the anisotropy of L r , applying the affine formulation is less straightforward than the isotropic secant procedures.To get round this drawback, two simplified isotropic versions have been proposed.In the first variant referred to as AFF-ISOT (Chaboche and Kanoute ´, 2003), L r is defined by L r ¼ 3k r J þ 2l r tgt ðe r eq ÞK.In the second variant AFF-ISOI, it is defined as the projection of the actual transversely isotropic tangent tensor on the subspace of the isotropic tensors, namely L r ¼ 3k r J þ 2l r inv ðe r eq ÞK with l r inv ðe r eq Þ ¼ 4l r sct ðe r eq Þþl r tgt ðe r eq Þ 5 .

Potential based approaches
Second-order formulations.For both second-order formulations SOE-1 and SOE-2, the thermoelastic strain potential of the LCC follows Eqs. ( 47) and ( 48).Regarding the initial second-order formulation referred to as SOE-1 (Ponte Castan ˜eda, 1996), the effective potential W is estimated by the effective potential W T ¼ hw T i of the thermoelastic LCC, i.e., W ¼ W T .The mean value over the phase of the strain field e r is prescribed as the reference strain e r .The moduli are defined as L r ¼ L r tgt ðe r Þ. Making use of those prescriptions inside the definition of the thermoelastic strain-energy function, the following expression is generated for the effective potential We refer to Ponte Castan ˜eda (1996) for the details on this calculation.Note that the estimated local strain fields in the LCC are the same as for the AFF-ANI formulation.It should be noted that the discrepancies between the affine and initial second-order formulations are due to the intra-phase fluctuations of the local fields which are implicitly taken into account in the SOE-1 approach as shown in Eq. ( 37) of Masson et al. (2000).
For the modified second-order formulation SOE-2 (Ponte Castan ˜eda, 2002), based on stationarity conditions related to the evaluated effective energy, e r is still prescribed as the reference strain e r while the moduli tensors associated with strain-energy function ( 48) are defined by the generalized secant relation where the additional strains êr follow from the introduction of suitable error measures (Ponte Castan ˜eda, 2002).Those additional reference strains, used to characterize intraphase strain fluctuations, are assessed using the relations êr ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi .In this procedure, the tensors L r are still written in the general anisotropic form as given by Eq. ( 24).Finally, the effective strain-energy function is estimated by the following expression Unlike the aforementioned formulations AFF-ANI and SOE-1, the SOE-2 procedure explicitly accounts for the intra-phase strain fluctuation in the definition of the LCC for all the constitutive phases.Like the SOE-1 formulation, a dual version of the SOE-2 based on a stress (instead of strain) energy function is available in Ponte Castan ˜eda (2002).In the following, for the sake of simplicity, we only consider the second-order procedures based on the strain-energy function w(e).Lahellec Suquet formulation (LS).This formulation (Lahellec and Suquet, 2004) retains the energetic framework of the initial second-order formulation and modifies it such that the field formulation (r ¼ hri) is in exact agreement with the energetic formulation (r ¼ ow oe ).The strain-energy function is approximated by a third-order Taylor expansion around a reference strain e r and the cubic term is linearized around an additional reference strain êr .The approximate effective potential is then rendered stationary with respect to both variables e r and êr such that the following expressions are generated for the LCC moduli It results the following expression for the effective potential The main difference between this new approach and the AFF-ANI and SOE-1 approaches lies in the fact that the polarizations s r , which explicitly allow for strain field fluctuations, generate a ''softer'' LCC in the LS formulation.However, the moduli tensors do not explicitly depend on these fluctuations.

Numerical accuracy
In order to quantify the numerical errors generated by the proposed numerical scheme, we studied the influence of the following parameters on the macroscopic response: (i) the type of finite element, (ii) the number of finite elements, (iii) the value of the convergence criterion associated with the fixed point iterative procedure described in Section 3.3.1 and (iv) the value of the step used to compute the derivative of the effective potential.Only the main results of this study are reported here.To minimize the numerical errors, we selected the following values.All the ensuing computations have been carried out by using a quadratic element with 8 nodes and a mesh constituted of 1200 elements.The convergence criterion associated with the fixed point iterative procedure described in Section 3.3.1 is defined by jðp r 1 À p r 0 Þ=ðp r 1 þ p r 0 Þj 6 g and the value of the parameter g is set to 0.01.For ''potential based approaches'', the effective hydrostatic and equivalent stresses defined by Eq. ( 49) are computed from the following approximation for the derivative of a function f(x) In order to avoid oscillations which can occur especially when the contrast between the phases is infinite, h is set to 0.1.Once the values of the aforementioned parameters have been set, the numerical procedure described above generally provides very accurate results for which the numerical errors are very small and far lesser than the discrepancies observed between the various nonlinear formulations.However, for some specific cases, difficulties of numerical order have been encountered.For composites with soft inclusions-0 < c ¼ r p 0 =r m 0 6 1:8-or for highly nonlinear composites-m close to 0-the fixed-point iterative procedure does not always converge.For this reason, the case of composites with inclusions that are weaker than the matrix is not addressed.The lowest values of the work-hardening parameter m range from 0.01 to 0.04, depending on whether the convergence of the fixed-point iterative scheme is attained or not.
As mentioned in Section 2.3, although the constitutive behavior of the individual constituents are Ramberg-Osgood laws, they asymptotically behave as power-law with the same exponent m for sufficiently large strains.Accordingly, the effective flow stresses are computed for a sufficiently large value of the macroscopic strain in order to guarantee that Eqs. ( 7) or (8) are satisfied.For instance, to check Eq. (7b), we compute lnðr eq Þ and choose one of the values of the macroscopic strain for which an affine evolution of lnðr eq Þ with respect to lnðe eq Þ is observed.

Results and discussion
The present methodology is carried out for three different types of heterogeneous materials: (i) ''reinforced composites'' composed of a nonlinear elastic matrix reinforced by linear isotropic elastic particles, (ii) ''two-phase materials'' made up of a nonlinear elastic matrix and a nonlinear elastic particulate second phase, (iii) ''porous media'' with a nonlinear elastic matrix.
All the considered nonlinear elastic individual constituents are assumed to obey Ramberg-Osgood constitutive equations as defined in Eq. ( 2).First, we focus our attention on the influence of the linearization procedure on the macroscopic responses.For the reinforced composites and the two-phase materials, only the response r eq ¼ f ðe eq Þ is estimated since these materials may be considered as incompressible.For the investigated porous media, the macroscopic response is compressible and two-dimensional and is characterized by Eq. ( 6).Then, the study of the influence of the linearization formulation on the macroscopic response is supplemented by an analysis at the local scale where the evolutions of some first or second moments of the strain field as functions of m or f p are computed.For all the F.E. computations, the number of finite elements inside the matrix N m and inside the inclusion N p has been set to N m = 900 and N p = 300, respectively.In the hereafter reported illustrative results, the material parameters have been chosen as indicated in Table 1.Additional parameters such as f p , m and the two-phase contrast c ¼ r p 0 =r m 0 are specified when necessary.

Reinforced composites
The macroscopic responses r eq ¼ f ðe eq Þ derived from the various linearization formulations for f p = 0.3 and m = 0.1 are reported in Fig. 1.First, it is observed that the SEC procedure yields the stiffest response.As expected from Bornert and Ponte Castan ˜eda (1998) and Masson et al. (2000), the AFF-ANI approach is also too stiff.It is worth noting that a more usual but too hasty comparison for two different microstructures would have led to a more optimistic evaluation: using for the AFF-ANI approach the generalized self-consistent linear scheme (Christensen and Lo, 1979;Chabert et al., 2004)-known to provide a good estimate for a microstructure described by Hashin's composite spheres assemblage instead of the actual periodic microstructure-leads to a softer and apparently better macroscopic response (AFF-ANI-GSC estimate in Fig. 1).The variational approach provides a slightly softer response.Again, the accurate comparison of the VAR estimate with the exact solution shows a larger discrepancy than that reported by Suquet (1995) where the microstructure is different.This illustrates clearly the bias introduced in the comparisons, either by the change of microstructure, or by the approximation induced by the linear homogenization scheme.The simplified AFF-ISOT formulation, based on an isotropic approximation of the tensor of tangent moduli, excessively softens the macroscopic response such that it lies significantly below the exact nonlinear solution.This result qualifies the conclusions of Chaboche and Kanoute ´(2003) where again the comparisons rely on different microstructures: a periodic nonlinear solution and the Mori Tanaka linear model for the LCC.The other isotropic version (AFF-ISOI) of the affine approach provides a far too stiff response.From a pragmatic point of view, both  former results suggest the possibility to define an intermediate linearization procedure between the AFF-ISOT and the AFF-ISOI formulations which would be more efficient and still easy to implement.The two variants of the second-order formulation provide the closest estimates of the exact solution, SOE-1 being still too stiff and SOE-2 too soft.This confirms the great improvement due to the integration of local field fluctuations in these theories, with the price of a higher complexity.However, these results suggest also that the modifications introduced in the newer formulation might still not be optimal, since the induced softening is somewhat too strong.To that respect, the very recent new version of the second-order formulation proposed by Idiart and Ponte Castan ˜eda (2005) and aimed at improving the two previous ones, definitely deserves an evaluation with the present methodology.This is left for further investigations.The LS formulation yields a macroscopic response which is in quasi perfect agreement with the one derived by the initial second-order approach SOE-1.This probably results from the similarity of the additional term for the polarization introduced in the LS approach and the expression of the difference between the affine and the SOE-1 estimates for the overall stress (see Masson et al., 2000)).Note that these results still hold for a nonzero threshold r m y as illustrated in Rekik et al. (2005).In view of deriving more general trends, we also computed the evolution of the normalized effective flow stresses r0 r m 0 (see Eq. 7) as a function of either the work-hardening exponent m or the inclusions volume fraction f p .The results are reported in Fig. 2. We find roughly the same trends as in Fig. 1 with some slight modifications.In Fig. 2a, it is noted that the AFF-ISOT approach provides very accurate results for m P 0.2.For smaller m values, as already noticed in Fig. 1, the AFF-ISOT approach excessively softens the macroscopic response.It is also observed that the modified second-order procedure yields the best estimate for m close to zero.This shows its ability to relevantly account for the strain field fluctuations which are still more significant in the ideally plastic limit.When the inclusion concentration increases (Fig. 2b), it can be seen that the SOE-1 and LS procedures excessively stiffen the macroscopic response for large values of the volume fraction and even violate the variational upper bound for f p P 0.5.The modified second-order procedure, as expected (Ponte Castan ˜eda, 2002), satisfies the bound irrespective of the values of m or f p .Further, it provides the best estimate for the reinforced composite at the macroscopic scale.The fact that the SOE-1 procedure violates the variational bound near percolation is well known (Leroy and Ponte Castan ˜eda, 2001;Ponte Castan ˜eda, 2002) and is due to significant strains fluctuations which are not relevantly taken into account in the SOE-1 procedure.This suggests that the LS procedure does not take the strain field fluctuations into account in an optimal way since it still violates the bound near percolation.

Two-phase incompressible materials
The evolution of the effective equivalent stress with respect to the macroscopic strain is depicted in Fig. 3 whereas Fig. 4 reports the variation of the effective flow stress r0 as a function of either the work-hardening exponent or the volume fraction.
The overall trends are similar to those observed for the reinforced composites.However, it can be seen in Fig. 4b that all linearization formulations qualitatively reproduce the evolution of the effective flow stress even for high concentrations unlike the reinforced composite where the SEC, AFF-ANI, LS and SOE-1 estimates blow up for high volume fractions.This result is probably due to the fact that the elastic particles of the composite are almost rigid so that their stress field is more difficult to assess.It should be noted that our implementation of the AFF-ISOT did not yield acceptable values of the effective flow stress for m 6 0.2 or f p < 0.3 (these values are not reported in Fig. 4).Further, some oscillations can be observed in Fig. 4 for the secondorder procedures, resulting from the computation of the derivative of the effective potential as mentioned in Section 3.4.Results for a softer second phase which are not reported here are qualitatively similar; nevertheless, they have been only determined for c > 1.8 since, most of the time, the fixed-point iterative procedure of our numerical scheme does not converge when applied to the second-order and LS approaches for c 6 1.8.

Porous media
The two-dimensional macroscopic response described by Eq. ( 6) is investigated for three different types of loading according to the value of the macroscopic strain triaxiality ratio s e ¼ em eeq : an isochoric extension s e = 0 (that is e ¼ e eq ê), a purely hydrostatic extension (s e = 1) and a mixed extension (s e finite and nonzero).
Isochoric extension.(s e = 0) The macroscopic tensile curves r eq ¼ f ðe eq Þ and r m ¼ f ðe eq Þ are depicted in Fig. 5 for f p = 0.3 and m = 0.4.In Fig. 5a, the same trends as in the case of the reinforced composite can be observed except that this time the AFF-ISOT approach slightly overestimates the reference solution.Moreover, the initial second-order as well as the LS estimates yield better results than the modified second-order procedure.Note that the same results are found for a nonzero threshold r m y .These trends hold for any pore volume fraction (see Fig. 6b).In Fig. 6a, which depicts the variation of the effective flow stress r0 with respect to m, the main results shown by Figs.5a and 6b are confirmed.However, although the SOE-1 and LS approaches provide the best results for m > 0.15, the SOE-2 procedure is the only linearization formulation which reproduces the evolution of the exact solution for highly nonlinear behaviors (m close to zero).Thus, this latter case turns out to be a very discriminant evaluation test.In Fig. 5b, it is noted that the linearization formulations based on an isotropic approximation of the local constitutive behavior-that is the SEC, VAR and AFF-ISOT approaches -yield very poor estimates of the reference effective hydrostatic stress unlike the anisotropic linearization procedures such as the AFF-ANI procedure.The results related to the ''potential based approaches'' such as the second-order procedures have not been reported in Fig. 5b for purely technical reasons.
Pure hydrostatic loading.(s e = 1) The macroscopic response is then no more characterized by r0 but by the effective flow stress r0 defined in Eq. (8a).The evolution of the effective flow stresses r0 with respect to either m or f p is depicted in Fig. 7.As can be seen, all linearization formulations but the VAR and SOE-2 approaches fail to capture the evolution of the reference effective flow stress.It is well known (Bornert and Suquet, 2001b) that the nonlinear formulations which only make use of the per phase average of the strain field without accounting for strain field fluctuations are not able to provide relevant estimates for isotropic porous materials subjected to hydrostatic loading.Indeed, under the classical assumption that the deviatoric part of the matrix behavior is nonlinear and the hydrostatic part is linear, as described by Eq. ( 2), the LCC moduli coincide with the initial linear (thermo)elastic moduli of the matrix-the deviatoric part of the average strain in the matrix is zero-and therefore lead to a linear and very stiff effective response.In our present case, the microstructure is not isotropic because of the periodic spacial distribution of the pores.For that reason, the deviatoric average strain in the matrix is not exactly zero and can become significantly nonzero when the pore concentration increases.On the other hand, for low concentrations, pore interactions can be neglected so that the effective behavior is almost isotropic and therefore the deviatoric part of the average strain in the matrix tends to zero.In this limit, the predictions of all the linearization schemes which do not account for strain field fluctuations in the definition of the moduli of the LCC are excessively stiff (Fig. 7b) as in the case of purely isotropic materials.For larger volume fractions, the presence of a deviatoric part in the matrix average strain induces softer but still too stiff responses (Fig. 7).In addition, one can note that the LS procedure significantly improves on the SOE-1 estimate as shown in Fig. 7b, which was not the case in previous situations.Mixed loading.The variation of the effective flow stresses r0 and r0 as a function of the work-hardening exponent m and the pore concentration f p is depicted in Figs. 8 and 9, respectively, for a fixed value of the macroscopic strain triaxiality ratio s e = 1/3.The effective response is now characterized by the evolution of the effective flow stress with respect to both the deviatoric and the hydrostatic parts of the macroscopic strain.This is illustrated by the fact that, although the VAR procedure yields an upper bound for the effective potential, it does not do so anymore for the effective flow stresses r0 and r0 as it did for a purely deviatoric or hydrostatic extension: this can be clearly seen in Fig. 8 where the reference solution lies above the VAR prediction.
Again, due to their relevant account for strain field fluctuations, the only formulations which reproduce the evolution of the reference effective flow stress r0 in Fig. 8 are the VAR and SOE-2 procedures, the latter leading to more accurate predictions.In Fig. 9 which is concerned with the effective flow stress r0 r m 0 , the trends are similar to those observed for the reinforced composite except that now the AFF-ISOT, SOE-1 and LS procedures provide fairly good estimates even if they are still less accurate than the VAR and SOE-2 estimates.In Fig. 10, the variations of the effective flow stresses r0 and r0 , respectively, are reported as a function of the strain triaxiality ratio s e .Note first that the iterative fixed point procedure associated with the SOE-1, SOE-2 and LS estimates converges only for small values of s e .It is likely that more sophisticated numerical schemes could have been successful for other triaxialities.This is an open technical question which has not been investigated in this study.As can be seen in Fig. 10a, none of the linearization formulations manages to capture the evolution of the reference flow stress r0 .On the other hand, the VAR, SOE-2 and AFF-ISOT approaches provide accurate estimates for the effective flow stress r0 (Fig. 10b), although the AFF-ISOT approach fails to reproduce the plateau depicted by the exact solution.

Local fields
In order to gain a deeper insight in the evaluation of the various linearization schemes sustaining nonlinear homogenization procedures, we also carried out some comparisons at the local scale, focused on the first and second moments of the strain in the nonlinear composite, evaluated from their counterparts in the LCC.Let us stress the fact that, for all the stress-strain approaches and for the LS procedure which can also be interpreted as a stress-strain approach, it is fully consistent to approximate the local fields in the nonlinear composite by those in the LCC.However, for the potential-based second-order procedures (SOE-1 and SOE-2), the strain statistics are not simply given by their corresponding quantities in the LCC, but there are additional ''correction'' terms which could be evaluated with the general relations very recently derived by Idiart and Ponte Castan ˜eda (2006).As mentioned in introduction, these corrections terms are not taken into account in the present study.Their evaluation would require additional algorithmic developments in order to compute the corresponding partial derivatives.Such developments are left for further investigations.Therefore, the results derived in this section relative to the efficiency of the second-order procedures in terms of local predictions should be considered as restricted to the classical ''LCC-based local interpretation'' for which, once again, it is assumed that the strain statistics in the nonlinear composite are given by their LCC counterparts.
First, we considered the reinforced composite defined in Section 4.1 and computed the evolution of the equivalent e r eq strain average over the phases (r) and the parallel e r k À e r eq and perpendicular e r ?averages of the strain fluctuations over the phases (r) as functions of the work-hardening exponent m.For this reinforced composite, those quantities are very weak inside the particles and only the strain field average and fluctuations inside the matrix are reported in Fig. 11a.Although all the linearization schemes provide an accurate estimate of the matrix average strain of the exact solution as expected from the fact that the particles are almost rigid, none succeeds in reproducing the evolution of the matrix average of the reference strain field fluctuations which become significant for large nonlinearity.The reasons for that can be understood better from the analysis of the snapshots (not reported here) of both the parallel and perpendicular reference strain field fluctuations for large nonlinearity (m = 0.1).We noticed that the matrix strain field fluctuations of the reference solution strongly localize into a very thin band tilted at an angle of 45°relative to the cylinder axis and located near the matrix-inclusion interface.Thus,Figs. 11a.ii and 11a.iii show that all the linearization schemes fail to capture this very strong localization of the strain field fluctuations.This conclusion has to be tempered especially when applied to the SOE-2 procedure and, to minor extent, to the AFF-ANI and LS approaches.Indeed, it is observed in Fig. 11a.ii that these procedures partially reproduce the evolution of the exact solution thus suggesting a weak localization of the parallel strain field fluctuations in the matrix.These results are not fully consistent with the macroscopic results depicted in Fig. 2a since the SOE-2 approach which relies strongly on the strain field fluctuations through the definition of the generalized secant modulus provides a very accurate estimate of the macroscopic response for large nonlinearity despite the fact that the strain field fluctuations are poorly assessed.This apparent inconsistency may probably be explained by the fact that for large non linearity, the work-hardening is weak and therefore the generalized secant modulus does not evolve significantly for large enough strains.
For the two-phase material defined in Section 4.2, the evolutions of e r eq , e r k À e r eq , e r ?as functions of the second phase volume fraction f p are reported in Fig. 11b.Again, it is observed that no linearization procedure accurately reproduces the evolution of the matrix average strain even if they all capture the appearance of the  plateau near the ''percolation'' limit.The best results in this case are provided by the SOE-2, AFF-ANI and LS procedures.Regarding the evolution of the matrix average of the strain field fluctuations, it should be noted that only the SOE-2 approach provides very accurate results.However, all linearization schemes reproduce qualitatively the evolution of the exact solution near the percolation limit.Note that all results are fully consistent with those observed at the macroscopic scale in Fig. 4b.For the porous material defined in Section 4.3 and subjected to an isochoric extension, Fig. 12 depicts both the evolutions of the strain averages over the phase-the matrix and the pore-and the evolutions of the matrix average of the strain field fluctuations as a function of the pore concentration f p .As can be seen, the AFF-ANI and LS procedures are in very good agreement with the exact solution while the SOE-2 predictions deviate significantly from the reference solution, especially for large porosity.Thus, although the SOE-2 approach leads to accurate results at the macroscopic scale (Fig. 6b), it apparently fails in this case to provide an accurate estimate of the first and second moments of the strain field.This conclusion should however be tempered by the fact that the actual local fields predicted by the SOE-2 procedure differ from those in the LCC.Unlike the SOE-2 approach, the LS procedure yields accurate results both at the global and local scale.The trends observed for the SOE-2 procedure still hold for the AFF-ISOT approach.Further, this latter leads to the worst predictions at the local scale.

Synthesis
The main general conclusions which can be derived from the comparisons performed in the previous sections are threefold.(a) and pore e p eq (b) strain averages, of the parallel e m k À e m eq (c) and perpendicular e m ?(d) matrix average of the strain fluctuations as functions of the pore concentration f p .The results are normalized by the equivalent macroscopic strain e eq .Work-hardening parameter: m = 0.2.
• First, the approaches which take the strain fluctuations into account provide overall better predictions of the global and local responses than those which only account for the average quantities.This assertion is especially well illustrated, even for weak strain intraphase heterogeneities, by the case of a porous material under a purely hydrostatic loading or under mixed extension with a large triaxiality ratio.As already mentioned in Section 4.3, the average strain in the matrix is then close to zero for weak enough pore volume fractions and the account for strain fluctuations becomes essential to describe the nonlinearity of the material behavior; when use is only made of the per phase average of the strain field, the LCC moduli are close to the initial linear (thermo) elastic moduli of the matrix and therefore lead to a much too stiff effective response.• Second, the nonlinear formulations making use of an anisotropic thermoelastic LCC provide better results than those relying on an elastic LCC with isotropic moduli.Of course, this assertion makes sense if compared approaches incorporate the same degree of information regarding the local fields (only their first, or first and higher-order moments).• Third, an explicit dependence of the moduli on the strain fluctuations seems to be the most appropriate way to incorporate field fluctuations in the linearization procedure.Indeed, to that respect, it is noted an overall superiority of the SOE-2 approach and, to a more limited extend but clear in the porous case, of the VAR approach which are the only procedures which take explicitly into account the intra-phase field fluctuations in the definition of the moduli while the LS and SOE-1 formulations take them only into account either in the additional polarizations (LS procedure) or in an implicit way through the derivation of the macroscopic stress from the effective potential (SOE-1 approach).
Further more specific comments can be made from other results.For the three types of materials under study, the SEC approach always leads to the stiffest predictions at the macroscopic scale.Moreover, as soon as the strain heterogeneities become significant-namely for high nonlinearity or near the percolation limit-it fails to reproduce qualitatively the evolution of the macroscopic response.At the local scale, the SEC approach also provides a poor estimate of the strain field.As shown by these results, the SEC approach is much less efficient than all the other approaches and therefore its use is only recommended when the simplicity of the implementation is the primary concern.
Due to its bound status, the variational approach yields too stiff macroscopic predictions.However, it is the sole approach together with the SOE-2 procedure which almost every time qualitatively reproduces the evolution of the macroscopic response, even when the strain heterogeneity becomes very significant.At the local scale, its results are generally not in good agreement with the reference solution although they qualitatively reproduce the mean features of the reference solution.Nevertheless, it should be noted that all the results obtained at the local scale have been derived for situations where the strain heterogeneities are weak with respect to the mean values of the strain field over the phase and one could expect better predictions when the material develops significant strain fluctuations in a given phase with respect to the phase average of the strain, as for the porous material subjected to a purely hydrostatic or mixed extension.For practical use, the implementation of the modified secant approach is much simpler than that of the SOE-2 procedure and even if the VAR predictions are less accurate than the SOE ones, the VAR approach can still be used as a convenient substitute for isotropic enough responses (porous materials under isotropic or quasi-isotropic load).
For the AFF-ANI approach, as expected, the associated macroscopic response is too stiff and it even often violates the VAR bound.Although the AFF-ANI procedure improves on the SEC approach a lot since it uses a tangent linearization of the local behavior instead of a secant one, it provides poor estimates as soon as strain fluctuations are large, due to the fact that neither the tangent moduli nor the polarization tensors account for the strain fluctuations.Nevertheless, better results are obtained at the local scale at least for the first moments of the strain field and sometimes also for the second moments as shown in Fig. 12 where the AFF-ANI and LS estimates are in very good agreement with the reference solution; but for other cases such as the reinforced composite, the AFF-ANI approach fails to capture the evolution of the second moments of the strain field fluctuations as shown in Fig. 11a.As a whole, one can think that, compared with the other more elaborate procedures presented in this paper, the main advantage of the AFF-ANI approach lies in the fact that it can also handle materials for which the constitutive behavior is governed by two poten-tials such as nonlinear viscoelastic, elastoplastic or elastoviscoplastic materials.For practical use, since the AFF-ANI approach provides too stiff responses and is rather difficult to implement due to the anisotropy of the tangent moduli, it is advised to confine its use to such two-potential materials for which other procedures such as the potential based approaches cannot be applied.
The AFF-ISOT approach leads to ambiguous results and for this reason should be considered with special care.On the one hand, it can yield macroscopic predictions which are in very good agreement with the reference solutions such as for reinforced composites and two-phase materials or for porous media under isochoric extension, though the isotropic approximation generally induces an excessive softening; nevertheless, like the AFF-ANI approach, the agreement is poorer as soon as the strain field fluctuations become significant (especially, once again, for porous materials under pure hydrostatic or mixed extension).On the other hand, it can generate local strain and stress fields which are very far from the reference solution, e.g., for porous media subjected to an isochoric extension (Fig. 12).This apparent contradiction can probably be explained by the addition of compensating approximations at the local scale which finally yield a good estimate of the macroscopic response.For instance, Fig. 12a shows that the matrix strain average is underestimated thus leading to a matrix tangent modulus which is too stiff.But the isotropic simplification which consists in replacing the secant modulus l r ¼ l r sct ðe r Þ associated with the perpendicular component of the strain field by the tangent modulus k r ¼ l r tgt ðe r Þ with k r < l r , softens the matrix tangent modulus.The superposition of these two opposite effects can yield a satisfactory macroscopic estimate.Accordingly, the AFF-ISOT is a convenient treatment because of its simple and easy to implement formulation and its apparent good performance in terms of overall predictions, but only in one dimensional situations, as for reinforced composites or two-phase incompressible materials for which earlier evaluations had led to similar observations (Chaboche and Kanoute ´, 2003).However, in more general situations with a multidimensional response-as for the porous material-this scheme leads to poor predictions and cannot be recommended.
For the reinforced composites, the two-phase materials and the porous materials under purely deviatoric extension, the SOE-1 procedure gives very accurate results when the strain heterogeneities are weak enough, with a slight overestimation with respect to the reference solution.Despite the fact that the SOE-1 estimate takes the strain fields fluctuations implicitly into account, it no longer captures the evolution of the reference solution when the strain heterogeneities become significant.This is illustrated by most of the curves reported for the porous material.It can also be observed for the reinforced composite near the percolation limit as shown in Fig. 2b.This shortcoming of the SOE-1 approach led to the development of the SOE-2 procedure which incorporates a direct dependence on the second moments of the field fluctuations in the phases.However, it should be noted that for some particular cases where the strain field fluctuations are weak, the SOE-1 approach leads to more accurate estimates of the reference solution than the SOE-2 approach, as shown in Figs.2a and 6a for m P 0.2.At the local scale, it is known that the SOE-1, within the context of the ''LCC-based local interpretation'' of the models, and AFF-ANI estimates lead to the same local strain fields predictions and thus the comments given for the AFF-ANI approach are still valid.
For the LS procedure, we find similar conclusions as for the SOE-1 method.Both approaches yield estimates which are often in excellent agreement although the LS slightly improves on the SOE-1 predictions when the fluctuations are essential (Fig. 7).Note that the LS procedure, unlike the SOE-1 approach, takes the strain fields fluctuations explicitly into account and still does not capture the trends of the exact solution for significant fluctuations (Fig. 2b).This suggests that the account for field fluctuations by means of the additional polarizations is not optimal and can still be improved.At the local scale, it is observed that the estimates of the first and second moments of the local strain field derived from both the LS procedure and the AFF-ANI approach are in very close agreement.Therefore, the conclusions regarding the AFF-ANI approach could also be applied to the LS procedure.
The SOE-2 procedure provides the best results at the macroscopic scale.It is the only approach which accurately reproduces the evolution of the exact solution whatever the material considered even for large nonlinearities, near to the percolation limit or for hydrostatic or combined extensions, thus showing that the field fluctuations are incorporated in a proper way.Although the SOE-2 formulation improves at lot on the initial second-order approach when the strain field fluctuations are significant, the effective response is a little bit too soft.Moreover, for some situations where the strain field fluctuations are weak, the initial second-order method provides better results than the current SOE-2 version (Figs. 2a and 6a for m P 0.2).This suggests that the SOE-2 approach could still be improved.Note that new results recently obtained by Idiart and Ponte Castan ˜eda (2005) and Idiart et al. (2006) have shown that the choice of the reference strains e r , which are the only variables not optimized in the construction of the method, is crucial.These authors proposed the use of a macroscopic strain as the reference instead of the phase average of the strain in the computation of the generalized secant moduli (Eq.( 52)).Moreover, the macroscopic predictions associated to this new prescription are in general stiffer than those of the SOE-2 version and softer than those of the SOE-1, at least in the case considered in Idiart and Ponte Castan ˜eda (2005) and Idiart et al. (2006), thus suggesting an improvement of the SOE-2 procedure since, as shown in the previous sections, the SOE-2 predictions are a bit too soft.The results derived from the LCC-based local interpretation of the SOE-2 approach are more qualified at the local scale.Indeed, when applied to the two-phase material, this procedure provides a very accurate estimate of the strain field fluctuations as shown in Fig. 11b, which confirms the results derived by Idiart et al. (2006).However, when applied to composites with infinite contrast such as the reinforced composite or the porous material under deviatoric extension, it fails to capture the evolution of the first-moment of the strain field in the matrix (Fig. 12a) and also that of the second moments of the matrix strain field .On the whole, the SOE-2 procedure yields the more accurate predictions among all the presented linearization procedures, as expected from the fact that it is the only one that combines all the general main features for an accurate scheme suggested at the beginning of this section: use of an anisotropic thermoelastic LCC, dependence on field fluctuations, explicit dependence of the moduli on the latter.The observed inadequacies at the local scale are likely to originate from the classical interpretation of the local fields in the LCC.These local predictions might be improved by the use of the newer and more consistent evaluation of field statistics proposed by Idiart and Ponte Castan ˜eda (2006).

Conclusion
The main objective of this paper was to carry out an unambiguous and systematic evaluation of various linearization treatments sustaining nonlinear homogenization procedures.To perform this evaluation, we presented an objective methodology which eliminates different bias which are often included in such evaluations.Those bias are mainly induced by the use of linear classical closed form estimates to evaluate the behavior of the LCC which may not be adequate enough for the considered materials, by the implicit or explicit consideration of different microstructures or by numerical or statistical convergence issues relative to the computation of the reference nonlinear solution.To remove the aforementioned limitations, we chose a composite material with a periodic microstructure so as to make possible a numerical resolution of the nonlinear problem and, after application of various linearization approaches, to be left with different linear homogenization problems with the same periodic microstructure, which can also be solved exactly.Then we proceeded to the evaluation of various classical or more recent linearization procedures, mostly in terms of predicted overall properties.A few comparisons have also been performed at the local scale, under the classical assumption that the local fields in the nonlinear composite can be approximated by those in the LCC.The application of such a methodology on three different types of materials with various local phase contrasts has led to provisional conclusions which are summed up in Section 4.5.
The evaluation has intentionally been limited to the simplest three-dimensional microstructure for which the proposed procedure applies.It can be extended to more general microstructural morphologies, by simply changing the geometry of the unit cell, or to other overall loads, such as for instance pure shear.In this case, three dimensional meshes might be necessary, but still with computational costs compatible with nowadays computers.Other constitutive behavior including elastoplasticity or elastoviscoplasticity could also be addressed.A more detailed evaluation of the local performances of recent formulations based on the prediction of effective potentials would require the implementation of the general relations derived by Idiart and Ponte Castan ˜eda (2006).This is left for further investigations.
The conclusions drawn from the presented comparisons are not definitive.However they allow us to establish the general superiority of recent formulations, taking into account field fluctuations and making use of anisotropic thermoelastic comparison composites, with respect to older treatments using elastic isotropic LCCs and ignoring intraphase heterogeneities.They suggest also the existence of room for even more efficient procedures, still to be developed.

InclusionFig. 4 .
Fig. 4. Two-phase composite: variation of the normalized effective flow stresses r0 r m 0 associated with the various linearization procedures with respect to (a) the work hardening exponent m and (b) the second phase volume fraction f p .Material parameters: c ¼ r p 0 =r m 0 ¼ 5.In (a) f p = 0.25; in (b) m = 0.2.

Fig. 8 .Fig. 9 .
Fig. 8. Porous material under mixed loading s e ¼ 1 3 : variation of the normalized effective flow stress r0 r m 0 with respect to (a) the work hardening exponent m and (b) the pore concentration f p .In (a) f p = 0.3; in (b) m = 0.2.

Fig. 10 .
Fig. 10.Porous material under mixed loading: variation of the normalized effective flow stresses (a) r0 r m 0

Fig. 11 .
Fig. 11.Variation of the matrix strain average e m eq (i), of the parallel e m k À e m eq (ii) and perpendicular e m ?(iii) matrix average of the strain fluctuations for a reinforced composite (a) and for a two-phase material (b) as a function of either the work-hardening exponent m or the inclusion concentration f p .The results are normalized by the equivalent macroscopic strain e eq .Material parameters: for the reinforced composite (a): f p = 0.3; for the two-phase material (b): m = 0.25, c ¼ r p 0 =r m 0 ¼ 5.

Table 1
Parameters used for the finite element calculations