A critical evaluation of local ﬁeld statistics predicted by various linearization schemes in nonlinear mean-ﬁeld homogenization

This paper is devoted to the evaluation of various classical and more recent linearization schemes for nonlinear homogenization in terms of their eﬃciency to characterize local ﬁeld ﬂuctuations in nonlinear heterogeneous composites. It relies on an unbiased comparison between ﬁeld statistics predicted by homogenization theories and those obtained from a reference solution solved using ﬁnite element techniques, based on the same microgeometry and boundary conditions and in which local nonlinear constitutive relations are exactly veriﬁed at each point. Two categories of linearization methods have been investigated: classical approaches based on a ”stress-strain” approach (classical secant, classical and simpliﬁed aﬃne) and methods based on ”variational principles” (variational and Lahellec-Suquet procedures). For each approach, the maps and the statistical distribution functions of the local ﬁelds (strain, stress and incremental work) illustrating the intra-and inter-phase heterogeneities are provided for reinforced and porous power-law composites. This study supplements an earlier study focused on comparisons at the global level [38, 39] and provides additional information on the accuracy of some available classical and recent linearization procedures. The proposed methodology gives access to a deeper insight on nonlinear homogenization schemes and may eventually lead to improvements of these formulations.


Introduction
Recent analyses of nonlinear composites or polycrystals [16,19,24,28] suggest that in order to reproduce accurately their effective behavior with homogenization theories, it is important to capture correctly the heterogeneity of the strain (or stress) field at local scale. In addition, estimates for the statistics of the local fields within the composites could be useful to account for the evolution of the microstructure which can strongly affect the macroscopic behavior of viscoplastic composites or polycrystals undergoing finite deformation [8,10,23,24]. Moreover, information on the distribution of (residual) stresses in random composites or polycrystals may improve the prediction of the macroscopic toughness and strength of such materials [21]. In addition, estimates of local stress or strain field statistics may allow one to develop statistical theories of damage nucleation and/or growth in the context of either brittle or ductile damage [2,3,26]. In practice, estimates of local fields can be obtained by means of numerical simulations carried out using either finite element methods or the fast Fourier transform (FFT)-based algorithm developed by Moulinec and Suquet [28,29,30,31]. Such approaches induce however computational difficulties due to the large size of the numerical systems to work out, as well as issues relative to the representativeness of the computer-generated microstructures, the appropriate averaging of the results and the choice of particular boundary conditions, both in the linear case [20] and, even more critically, in the nonlinear one [13]. These approaches imply large computational expenses and may in practice often be restricted to two-dimensional problems especially when a large number of microstructural configurations or loading conditions are to be investigated. Alternatively, spatial distributions of kinematic fields can be experimentally measured using full-field strain measurement tools based for instance on digital image correlation techniques [6,14,40,41,43]. Such experiments however are very complex to carry out and usually restricted to surface observations and the statistical representativeness of the measured field fluctuations might be hard to ensure. Techniques to characterize local stress fields are still in their very early stage of development [9]. Such numerical or experimental tools are by essence restricted to the characterization of the local field statistics of a limited number of samples and loading conditions. They are very efficient tools to investigate the complex interactions governing the development of field fluctuations and to propose models to describe them. But because of their complexity they can not be used to characterize local field fluctuations in the context of structural computation, for which it is necessary to adopt some less demanding tools to describe the local field heterogeneities. To this end, nonlinear mean-field theories seem to be particularly attractive modeling tools to extract low-order statistical information about the statistical distribution of local strains in the context of multiscale modelling of various classes of materials.
In the context of linear composites, well-known exact formulas are available to express the first and second moments of the local fields in the phases, in terms of the effective elastic properties [4,32]. Such formulas are useful as they allow the extraction of estimates for the statistics of the local fields from the sole homogenization estimate for the effective thermoelastic potential of the composite. Recently, these formulas have been generalized to higher order statistics in the context of nonlinear homogenization. They can be applied to nonlinear homogenization schemes which generate an estimate for the effective nonlinear potential of the composite, such as those derived from variational formulations. For schemes based on the concept of linear comparison composite (LCC) [17,18], it is possible to compare by means of these formulas the statistical moments of various orders of the local fields in the nonlinear composite, as predicted by such nonlinear homogenization scheme, to the corresponding moments in the associated linear comparison composite subjected to the same load. It is shown that for some nonlinear schemes such as the secant variational formulation (here referred to as VAR) [34,42] or the model proposed by Lahellec and Suquet (LS) [22], the first and second moment of the fields in the nonlinear composite and the LCC do coincide. But this no longer holds in the context of the tangent second-order [35] and second-order [36,16] methods, where "corrections" terms arise [17,18] because of the lack of full stationarity of the relevant functionals with respect to the properties of the LCC. Other more classical nonlinear formulations rely on the ad hoc assumption that local fields in the LCC are good estimates of the fields in the nonlinear composite and deduce effective properties, not from an effective potential, but from appropriate averaging procedures of these local fields. This is the case of the classical secant [1] (referred to as SEC), the affine schemes (AFF-ANI) [27] and its isotropic simplification (AFF-ISOT) [12]. For such schemes, statistics of the local fields are naturally identical to those in the LCC; it should however be noted that these schemes do generally not allow to define an effective potential.
In this work, we focus our attention on the category of nonlinear homogenization theories, for which it is possible to investigate the local fields up to second order statistics in the LCC defined by these methods as consistent approximations for the corresponding quantities in the nonlinear composite. This category avoids additional analytical or numerical expenses related to the interpretation of the local field statistics in the LCC as approximations of the corresponding quantities in the nonlinear composite. In this context, it is then easy to extend the non-biased methodology presented in [39] to evaluate the accuracy of nonlinear homogenization schemes in terms of effective properties to the evaluation of their ability to predict local field statistics. We recall that this methodology relies on the analysis of a periodic composite for which both the nonlinear homogenization problem and the linear homogenization problems associated with the chosen linear comparison composite (LCC) can be solved for exactly the same microstructure, for the same loading conditions and with a similar level of numerical approximation, so that the effects of the sole linearization scheme can be evaluated without ambiguity. Both the problem of the initial nonlinear periodic unit cell and the linear one associated with the LCC are solved with high accuracy and low computational expense using classical finite element techniques. Unlike numerous previous studies, there is no approximation related to a change of microstructure between the nonlinear problem and the linear ones related to the LCC derived from various available linearization methods. This methodology allows comparisons between the local field statistics of the fully nonlinear problem and those of the LCC. For the considered linearization methods for which an effective potential can be defined, namely the VAR and LS procedures, the proposed methodology allows the computation of secondorder field statistics without any obligation to handle the effective potential of the LCC. In the literature, similar studies i.e. comparisons between field-statistics derived from full-field information and mean-field theories, have been performed but for a limited number of nonlinear linearization methods, such as Suquet [42] and Bornert [6] who made use of the VAR method to extract the isotropic trace of the second moments of the strain in two-phase, elasto-plastic composites, and compared the predictions with experimental results. Brenner et al [8] also studied field statistics in viscoplastic polycrystals by means of the "affine method". Field fluctuations for two-phase composites with "particulate" microstructures were also studied by Moulinec and Suquet [30,31], by means of the VAR homogenization method and full-field numerical simulations. More recently, Rekik [38] has provided comparisons between numerical simulations and most available linearization methods including the first version of the second-order method [35] denoted hereafter by the tangent second-order [16,17,18]. Pierard et al [33] compared the global and local responses of elasto-plastic composites reinforced with elastic and aligned ellipsoids derived both by the secant (classical and second-order methods) and incremental approaches and by Finite Element simulations performed on a random representative volume element. Idiart et al [15] proposed corrections terms for the assessment of the field statistics for the tangent second-order and a sophisticated version of this method. They also made comparisons at the local level between the numerical fields provided by the FFT method, the tangent second-order and its improved version. More recently, Idiart et al [19] made comparisons between the second-order homogenization estimates and the FFT full-field numerical simulations for the first-and second-order field statistics of the stress and strain fields.
In this paper, the proposed methodology can be relevant for a rigorous assessment of local field statistics up to second order for various linearization methods, especially those based on a "stressstrain" approach (such as classical secant and affine methods) as well as "potential-based" approaches such as the VAR and LS models, for which the "stress-strain" variant of their formulations coincide with the potential-based ones. In particular the accuracy of the LS linearization is evaluated at the local scale thus complementing its global evaluation in [38,39]. On the other hand the evaluation at the local scale of the variants of the second-order method is not addressed as similar comparisons are performed in other studies such as Idiart et al [17,18]. Finally, the proposed methodology can be extended to handle large classes of problems, whatever the choices of the local constitutive laws, the morphology of the periodic microstructure and the type of loading.
The tensor notation used herein is a fairly standard one. Products containing dots denote summation over repeated indices. For example, L : ε = L ijkl ε kl e i ⊗ e j and E :: F = E ijkl F klij where e i (i = 1, 2, 3) is a time-independent orthogonal cartesian basis and the operation ⊗ denotes the classical tensorial diadic product. Cylindrical coordinates (r, θ, z) will be used as well, with e r , e θ , e z = e 3 being the unit vectors of the orthonormal cylindrical basis and u r , u θ , u z the cylindrical coordinate components of a 3D vector u. Use will also be made of the fourth-order tensors K = I − J and J = 1 3 i ⊗ i which are the usual projectors on the subspaces of purely spherical or deviatoric second-order tensors. The tensors i and I are second and fourth-order symmetric identity tensors.

A brief outline of the considered linearization schemes
The main objective of homogenization is to predict the macroscopic behaviour of composite materials in terms of the behaviour of their constituents and prescribed statistical information about their microstructure. In this framework, we consider composite materials made of N different homogeneous constituents which are assumed to be randomly distributed in a specimen occupying a volume V , at a length scale that is much smaller than the size of the structure and the scale of variation of the loading conditions. The constitutive behaviour of each phase is characterized by a strictly convex single potential or strain energy function w r (r = 1, ..., N ), such that the stress σ and strain ε tensors are related by where ∂ ∂ε denotes differentiation with respect to ε, and the characteristic functions χ r (x) serve to describe the microstructure, being 1 if the position vector x is in phase r, and 0 otherwise. This constitutive relation corresponds to a nonlinear elastic behavior within the context of small strain. However, equation (1) can also be used either within the context of the deformation theory of plasticity or for viscoplastic materials, in which case σ and ε are the Cauchy stress and Eulerian strain rate, respectively. In the following, the constitutive behavior of the individual constituents of the composite is assumed to be nonlinear elastic. Let < . > and < . > r denote the volume averages over the composite (V ) and over phase r (V r ), respectively. The effective behaviour of the composite, which is defined as the relation between the average stressσ =< σ > and the average strainε =< ε > can be characterized by an effective strain potentialW , such that where κ(ε) = u, ε(u) = 1 2 ∇u + T ∇u and < ε(u) >= ε . The nonlinear problem (2) can be approached using nonlinear homogenization theories based on two steps: linearization and linear homogenization. One way to apprehend what these two steps consist of is first to notice that the local strain field ε(u) solution of the variational theorem (2) 2 solves the so-called nonlinear local problem consisting of the following set of equations where f r (ε) = ∂w r (ε) ∂ε denotes the constitutive law of phase r, and second to approximate the nonlinear problem (3) by the following system where L r {ε(x), x ∈ V } and τ r {ε(x), x ∈ V } are known functionals of the local strain field ε(x) of the local linear problem (4a). Problem (4a) coincides with the so-called local problem associated with a fictitious linear composite usually called the Linear Comparison Composite (LCC), the behavior of each of its individual constituents being linear thermo-elastic and defined by σ = L r : ε + τ r . The linearization stage consists in linearizing the phase constitutive relations σ = f r (ε) of the actual nonlinear composite and therefore providing analytical expressions of the functional L r {ε(x), x ∈ V )} and τ r {ε(x), x ∈ V )} while the linear homogenization stage aims at solving the local linear problem and deriving its effective behavior by means of linear homogenization techniques. The effective behavior of the LCC -which varies nonlinearly with respect to the macroscopic strain ε -is interpreted as an estimate of the nonlinear effective constitutive law of the nonlinear composite (2) 1 . Of course, for each macroscopic strain ε, the local problem associated with the LCC has to be solved several times until the nonlinear relations (4b) do not evolve anymore. System (4) can be solved by means of a fixed-point iterative procedure which is iterated on a set of reference strains which are computed from the local strain field in the LCC (ε(x), x ∈ V ). The definition of these reference strains depends on considered linearization procedure and will be specified for each of them in Sections 2.1 and 2.2. In the following the main concepts of the investigated linearization methods are recalled. Each phase constituting the LCC is therefore assumed to follow a (thermo-)elastic constitutive law σ = L r : ε + τ r , with a uniform tensor of moduli L r and a pre-stress or polarization tensor τ r .
Though the considered linearization schemes and the methodology developed here could be applied to more general situations, the presentation is restricted, for the sake of clarity, to composites made of isotropic constitutive phases, with local constitutive relations of the form: where σ m = σ : i/3 (resp. ε m = ε : i/3) and σ eq = 3 2 σ : K : σ (resp. ε eq = 2 3 ε : K : ε) are the isotropic part and the von Mises measure of the deviatoric part of the stress (resp. strain) tensor. The deviatoric unit second order tensorê is parallel to the deviatoric part of the strain. In such a constitutive relation, only the deviatoric parts of the stress and strain tensors are concerned with the nonlinearity, while the isotropic parts obey a linear relation, characterized by the bulk modulus k. The particular loading conditions, microstructure and form of the nonlinear relation σ eq = f (ε eq ) considered in the numerical simulations will be described in section 3.2.

Elastic LCC
Some linearization schemes define an elastic LCC (τ r = 0). The linearity of the local problem ensures a linear dependence of the local strain field in the LCC with the macroscopic strainε which reads where A(x) is the fourth-order localization tensor. The per-phase average of relation (6) reads ε r =< ε > r = A r :ε where the fourth-order tensor A r =< A > r is not necessarily symmetric and only depends on the microstructure and the local properties, but not on the overall load.

Classical secant (SEC)
For this formulation, the constitutive behaviour per phase inside the LCC reads: where L r is defined as the isotropic tensor of secant moduli evaluated at the reference strain ε r , 3ε r is the secant shear modulus of the nonlinear constitutive relation of phase r. According to the SEC method, the reference strain is defined as the equivalent per phase average of the strain in the LCC: ε r =ε r eq =< ε > r eq = (A r :ε) eq .

Modified secant extension (VAR)
This formulation coincides with the variational approach of Ponte Castañeda [34,42]. In the LCC, the constitutive behaviour of the phases is identical to that defined by the SEC method but the reference strain ε r is set equal to the scalar second-order moment of the strain in phase r in the LCC: ε r =ε r eq = ε 2 eq r . It is worth recalling that the variational procedure provides an upper bound for the effective potential.

Thermoelastic LCC
In a thermoelastic LCC (τ = 0), the local strain field and its per-phase averages read The fourth-order tensor A(x) coincides with the localization tensor field of the elastic LCC. The second-order tensor a(x) is the local strain under vanishing overall strain. In the framework of a two-phase composite and in particular for the matrix-inclusion microstructure considered hereafter, we apply Levin's relation [25] which allows one to compute a(x) from A(x) as follows with B = (△L) −1 : △τ and C =ε + B.
The superscripts "m" and "p" refer to the matrix and the inclusion phases considered hereafter, respectively.

Classical affine (AFF-ANI)
This formulation coincides with the original affine approach proposed by Masson et al [27]. The linear constitutive behaviour of the individual constituents of the LCC follows the thermoelastic law where L r and τ r are the tensor of elastic moduli and the polarization tensor defined by L r = L r tgt (ε r ) = ∂ 2 w r ∂ε 2 (ε r ) and τ r = ∂w r ∂ε (ε r ) − L r :ε r , respectively. Again ε r =ε r in the LCC. We note that the thermoelastic law of each individual constituent of the LCC is tangent to the nonlinear constitutive law of each phase of the real nonlinear composite σ = ∂w r ∂ε (ε) at ε =ε r . Unlike the secant methods, the tensor L r is anisotropic and can be written as where µ r tgt (ε eq ) = dσ eq 3dε eq (ε eq ) is the tangent shear modulus. The fourth order tensor E r = 2 3ê r ⊗ê r withê r = (K :ε r ) /ε r eq is the projector onto the direction of the strainε r and F r = K − E r . When ε r is uniaxial, L r is transversely isotropic with respect to the same axis.

Simplified affine (AFF-ISOT)
Due to the anisotropy of L r in the original affine formulation, the application of this scheme requires more analytical developments and numerical expenses than the linearization schemes defining an isotropic LCC. To simplify this formulation, Chaboche and Kanouté [12] have proposed a new variant referred to as AFF-ISOT with an isotropic tensor of elastic moduli defined by: L r = 3k r J +2µ r tgt (ε r )K, where ε r =ε r eq . Note that both relations coincide for local strains ε proportional toε r but not in general.

Lahellec and Suquet method (LS)
This formulation retains the energetic framework of the initial second-order formulation and modifies it such that the field formulation (σ =< σ >) is in exact agreement with the energetic formulation (σ = ∂w ∂ε ). The strain-energy function is approximated by a third-order Taylor expansion around a reference strain ε 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 ε r andε r such that the following expressions are generated for the LCC moduli: It is worth noting that the main difference between the LS and AFF-ANI approaches lies in the polarization tensor τ r , which in the case of the LS formulation, explicitly takes into account the strain field fluctuations and generates a "softer" LCC material (at least in the case of power-law materials). For both formulations, the tensor of elastic moduli L r is defined in the same way.

General definitions of the local field statistics
In this part we give general expressions of quantities useful for the evaluation of the per-phase averages of the strain field and its fluctuations both in the reference medium (i.e. the nonlinear composite) and in the various LCCs (elastic or thermoelastic). Similar expressions related to the stress field are reported in Appendix A. The so-called inter-phase heterogeneity of the stress and strain fields are characterized by their first-order momentsε r andσ r . In order to analyze the intra-phase fluctuations of the local field, it is appropriate to identify two "components" of the deviatoric strain (respectively stress) which represent its projection "parallel" ε (respectively σ ) and "perpendicular" ε ⊥ (respectively σ ⊥ ) to the deviatoric per phase average of the strainε r d , which in the present situation, as it will be shown in section 3.3, is proportional to the macroscopic deviatoric strainε d . These components can be derived through the following relations The fourth order tensor E and F are the classical projectors defined as previously, based on the overall load [35] and related to the directions parallel and perpendicular to the overall load, respectively. Note that the parallel and perpendicular components of the local fields satisfy the relations: ε 2 eq = ε 2 + ε 2 ⊥ and σ 2 eq = σ 2 + σ 2 ⊥ . Furthermore, one has ε (x) = |2ê : ε(x)/3|. Based on the parallel and perpendicular components of the fields introduced in (13), we define three different measures of the fluctuations of the local strain field over the phase by where C r ε = ε ⊗ ε r − ε r ⊗ ε r is the phase covariance strain tensor. In the sequel, δ r (ε), δ r ⊥ (ε) and δ r eq (ε) will be referred to as the parallel, perpendicular and isotropic (or equivalent) Measures of the Intraphase Strain Fluctuations (MISF) over the phase r, respectively. Their stress counterparts are defined in Appendix A.1.
To fully characterize the fluctuations of the strain field, it is necessary not only to compute its phase covariance tensors C r ε but also the spatial distributions of its fluctuations. We restrain the present study to the evaluation of the parallel δ r (ε(x)), perpendicular δ r ⊥ (ε(x)) and isotropic δ r eq (ε(x)) fluctuations of the strain fields defined by means of the trace of the strain local covariance tensor C r ε (x) = (ε(x) −ε r ) ⊗ (ε(x) −ε r ) with respect to the fourth-order projector E, F and K, respectively, as follows The measures of the strain fluctuations defined in (14) are the second-order moments of their local . However, it is worth noting that they are not the standard deviations of the strain components ε (x), ε ⊥ (x), ε eq (x) introduced in (13) except for the parallel measure δ r (ε) for which we have δ r (ε) = ε 2 r − ε r 2 .
In addition to the computation of the spatial distribution of the local strain field fluctuations δ r (ε(x)), δ r ⊥ (ε(x)) and δ r eq (ε(x)), we also determine their probability density functions. The computations are carried out by means of a procedure described in Appendix B which allows one to smooth a probability density function when the numbers of available values in the classes are not so high.
To supplement these investigations of the fluctuations of these local kinematic and static quantities, classically considered in the context of nonlinear homogenization, we focus also our attention on an energetic quantity, namely the local density of the incremental work σ : dε due to a change of overall load fromε toε + dε. By Hill-Mandel's macrohomogeneity principle, one hasσ : dε =< σ : dε >. The spatial distribution of this quantity is the signature of the local contributions of the medium to the overall work, for the considered loading history, and the statistical characterization of its fluctuations is another micromechanical signature of the way a material reacts to an overall load. The ability of a nonlinear homogenization scheme to accurately predict such an information is an important factor of its overall performance. Note that in the present case of nonlinear elasticity, without dissipation, this quantity also reads σ : dε = ∂w ∂ε : dε = dw(ε). It is however numerically easier to compute it from the local stress and strain fields. In practice, it is evaluated for finite but small loading increments and is averaged over each element e of the finite element mesh according to where subscripts (i) and (i+1) refer to two successive loading steps -the averages over the elements are carried out through the procedure INTG of the FE software Cast3M used in this study [11]. As loading increments are small and identical in the nonlinear solution and the solutions based on the tested nonlinear homogenization procedures, this approximation is sufficient.

The nonlinear periodic composite
In order to avoid large numerical costs and the computational difficulties mentioned in introduction when solving computational homogenization problems associated with complex microstructures, we consider as a first illustration of the proposed methodology a two-phase composite with a simple periodic microstructure. The periodic composite is made of identical isotropic linear elastic spherical inclusions embedded in a nonlinear isotropic matrix. The matrix is assumed to obey a Ramberg-Osgood equation defined by where k m and µ m e are respectively the bulk and shear modulus of the elastic part of the matrix constitutive law. Both moduli are chosen to be finite in order to avoid numerical difficulties related to infinite tangent or secant modulus at ε = 0 and to material incompressibility. Besides, m = 1/n is the strain-hardening parameter (n is the nonlinearity exponent) such that 0 ≤ m ≤ 1, ε 0 is an auxiliary strain, σ 0 is the flow stress. The constitutive relation σ eq = f r (ε eq ) is obtained by means of a numerical inversion procedure. Matrix and inclusion are assumed to be perfectly bonded, thus ensuring the continuity of the displacement and the stress vectors at the inclusion/matrix interface. Note that such type of composite enables to study the case of a rigidly reinforced or porous material by assigning to the bulk k p and shear µ p moduli of the inclusion numerical values close to infinity or zero, respectively. The spherical inclusions are distributed according to a hexagonal network in the transverse plane and are aligned along the third direction, such that a cylinder with a hexagonal basis with a single spherical inclusion located at the cylinder's center can be used as unit cell. The macroscopic loading conditions that will be considered are axial symmetric along the third direction. In a first and classical approximation [37], this hexagonal cell can be replaced by a cylinder with a circular basis, and the periodic boundary conditions (BC) replaced by symmetry conditions, making the whole unit cell problem fully invariant with respect to any rotation along this direction. This 3D homogenization problem then reduces to a 2D axial symmetrical problem which can be solved at a particularly low numerical cost. The unit cell being in addition symmetric with respect to the transverse plane, only one fourth of the cross section of the cylinder needs to be meshed. Symmetry conditions are prescribed on the transverse plane and along the symmetry axis. Homogeneous longitudinal (resp. radial) displacements are imposed on the upper (resp. lateral) face to ensure (quasi-)periodic BC -see Eq. (14) in [39] for a detailed expression of the BC. With such a simplification of the boundary conditions, the considered nonlinear FE solution will no longer provide the exact fields in the periodic composite. However the corresponding approximations are likely to be small. Most importantly, it should be emphasized that exactly the same approximation is used in the nonlinear solution and in the LCC, so that the only difference between both calculations is the substitution of the exact local nonlinear constitutive relation by the linearized, per-phase uniform, approximate relation defined by the considered linearisation scheme. For this reason, the nonlinear solutions will be sometimes referred to in the following as the "exact" result (with quotation marks). Indeed it does not provide the exact fields in the nonlinear periodic composite, but gives the local fields in the axisymmetric unit cell associated with the exact nonlinear local constitutive relation, which serve as references to the approximate fields in the same cell under the same conditions but associated with linearized local constitutive relations. In the ensuing calculations, the unit cell is submitted to a monotonic uniaxial purely deviatoric extension along the symmetry axis, such that the overall strain isε =ε eqê whereê is reduced here to the axisymmetric tensor e 3 ⊗ e 3 − 1 2 (e 1 ⊗ e 1 + e 2 ⊗ e 2 ).

Computations of the field statistics
Hereafter, we address more technical issues relative to the computation of the reference exact field statistics in the nonlinear composite ( §3.3.1) and the predicted ones in the elastic or thermoelastic LCC ( §3.3.2) derived by the various considered linearization schemes. 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: • The deviatoric strain average over the phase r readsε r d = K :ε r =ε r eqê ; thereforê • 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, |c d | = c eq while |c d | = 2 3 c eq if c is a stress.
• Symmetric second-order tensor fields (such as local stress and strain fields) are axisymmetric c(x) = c rr (r, z)e r ⊗ e r + c rz (r, z)(e r ⊗ e z + e z ⊗ e r ) + c zz (r, z)e z ⊗ e z + c θθ (r, z)e θ ⊗ e θ .
• 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 [7], which is derived from the decomposition of transversely isotropic tensors introduced by Walpole [44]. The definition of the tensors E L , J T , F 0 , K T and K L , is recalled in Appendix C.
Such decompositions and the associated invariance properties are valid for both the nonlinear problem and the quantities associated with the LCC. In particular, the tensors of moduli of the phases admit decomposition (18), with γ = γ ′ .

"Exact" results
The FE code Cast3M has been used for the determination of the local fields. An appropriate iterative procedure is used to resolve the nonlinear problem. After convergence at the imposed macroscopic strainε, the hydrostaticε r m and equivalent strainsε r eq are given bȳ Moreover, parallel, perpendicular and equivalent Measures of Intraphase Strain Fluctuations (MISF) read 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 functions of the four scalars ε zz , (ε rr + ε θθ ), ε eq and (ε :ê). Similar relations hold for the stress field. Lastly, the incremental work density δw e per mesh element e reads In the following numerical simulations, note that the loading strain increment between steps i and i + 1 is always equal to 2.5% of the final macroscopic strain which was set equal toε eq = 4%.

Local fields evaluated in the elastic and thermoelastic LCC
For the considered problem, the polarization tensor τ r is axisymmetric and exclusively deviatoric (τ r m = 0). Therefore it can be given the form τ r = 2 3 τ r eqê . Accordingly the tensor B is axisymmetric and purely deviatoric B = B eqê . Moreover, as the applied macroscopic strain is axisymmetric: ε =ε m i +ε eqê , the tensor C is axisymmetric and reads C = C m i + C eqê where C m =ε m and C eq =ε eq + B eq . Note that C eq , B eq andε eq are the algebraic scalars describing the deviatoric parts of the second order tensors C, B andε.
Substituting such a decomposition in equation (9), the local strain field induced by this loading in a thermoelastic LCC can be written in the form It is then a superposition of two independent local strains ε i (x) = A(x) : i and εê(x) = A(x) :ê induced in an elastic LCC by the specific axisymmetric loadingsε = i andε =ê, respectively. Accordingly, it is necessary to solve independently two elastic problems related to the elastic LCC submitted to the loadings i andê, respectively, in order to get access to their response in terms of the strain fields ε i (x) and εê(x). This will allow to assess field statistics (first-order moment and measures of the strain fluctuations) in the elastic and thermoelastic LCCs. It is noted that strain field statistics in the elastic LCC are deduced from that in the thermoelastic LCC by replacing C eq byε eq as the tensor B is zero in the elastic LCC. In the following we give expressions of J :: C r ε , E :: C r ε and I :: C r ε useful to compute the field statistics in the thermoelastic LCC (Eq. (14)), the traces K :: C r ε and F :: C r ε being obtained by means of K = I − J and F = K − E. To this end we define a fourth-order tensor ξ r and a second-order tensor A r so that I :: where the vector V is given by V = t C m , C eq and the components of the second-order tensor A r by A r aα =< e a : ε α > r = e a :< ε α > r where The second-order tensor A r aα defines the localization of the strain field in the space generated by the orthonormal basis {e a , a ∈ {i, e}} and is not symmetric. The components of the fourth-order tensor ξ r read ξ r aαβb = (e a : ε α )(ε β : e b ) r = (e a ⊗ e b ) ::< ε α ⊗ ε β > r where a, b ∈ {i, e} and α, β ∈ {I, E}. This tensor admits the following symmetry ξ r aαβb = ξ r bβαa and characterizes the second-order moments of strain fields induced by the two different loadings i andê. The quantities ξ r .αβ. are defined by ξ r .αβ. = (ε α : ε β ) r . The point in ξ r .αβ. can be interpreted as an implicit summation over the subscript a -i.e., ξ r .αβ. = a (e a : ε α )(ε β : e a ) r -when the orthonormal basis {e a , a ∈ {i, e}} = {e i , e e } is enlarged to {e i , e e , e f , e g } with e f = (e r ⊗ e r − e θ ⊗ e θ )/ √ 2 and e g = (e r ⊗ e z + e z ⊗ e r )/ √ 2. The enlarged orthonormal basis spans the space of symmetrical second-order tensor fields satisfying rotational invariance along e z . It is worth noting that B eq does not appear in the expressions of J :: C r ε , E :: C r ε and I :: C r ε . Again these integrations can be performed by the object-oriented operators of CAST3M.
Lastly, the incremental work density ∆w e in the thermoelastic LCC, from load incrementε toε+△ε associated with the states (i) and (i+1), respectively, is evaluated numerically by means of Eq. (16) where ε (i+1) (x) and σ (i) (x) are given by respectively (see Eqs. (22) and (A.4)). The incremental work density in the elastic LCC is deduced from that in the corresponding thermoelastic LCC by substituting the terms B eq and τ r eq by 0.

Illustrative example using the classical affine linearization scheme: case of a porous medium
To bear out the relevance of the proposed methodology in terms of prediction of the local field statistics stemming from linearization methods we propose in this section to focus on the assessment of a particular linearization method: the AFF-ANI scheme. The predictions of the affine linearization for the local equivalent strain ε eq (x), the parallel ε (x)−ε m eq and perpendicular ε ⊥ (x) strain field fluctuations normalized by the macroscopic equivalent strain ε eq are compared with the "exact" results derived from the numerical simulation of the periodic cell with a pore volume fraction of 25% and a low work-hardening exponent m = 0.15 implying a high nonlinearity for the matrix of the considered porous media subjected to a pure deviatoric extension.
At this point, it should be noted that numerical investigations carried out for different values of ε eq and m have shown that the applied strainε eq = 0.04 is representative of all range of possible macroscopic strains for the nonlinear part of the composite behavior, even for low work-hardening exponent. The maps illustrating the spatial distributions of the local considered fields for the exact and AFF-ANI estimations are reported in Fig. 1. The probability density functions quantifying the distributions of these fields in the matrix phase are plotted in the same figure. It is worth noting the complementary between the maps and density functions to interpret the local field statistics. In particular colormaps in the strain maps have been chosen to emphasize the spatial distributions and may not be the same for the exact and AFF-ANI results. Distribution functions plotted on the same graphs help to compare both field fluctuations. The evolutions of the equivalent average strainε m eq , the measures of parallel δ m (ε) and perpendicular δ m ⊥ (ε) intraphase strain fluctuations (MISF) over the matrix for different work-hardening exponents m and pore volume fractions f p are plotted in Fig. 2a and 2b, respectively, in order to study the effect of these parameters on the accuracy of the affine homogenization model for power-law porous media. In the sequel, the exact solution of the nonlinear composite is referred to as NL. Similar maps, probability density functions and evolutions are reported in Fig. 3 and 4 for the stress field in the same porous medium.
These investigations demonstrate that the affine linearization scheme predicts well the zones of strong localization of the local equivalent strain located both at the corners of the unit cell and at the boundary of the pore. The distributions of the parallel and perpendicular fluctuations of the AFF-ANI strain field are in satisfactory agreement with the exact (NL) solution. However, these quantities are overestimated by this model at the localization regions, namely a conical region tangent to the pore for the parallel strain fluctuations and the matrix/pore interface for the perpendicular strain fluctuations.
The evolutions of the normalized matrix parallel and perpendicular MISFs with respect to the work-hardening exponent show that the AFF-ANI scheme overestimates the strain fluctuations especially for strong and moderate non-linearities (0 ≤ m ≤ 0.4). This is in agreement with the fact that the AFF-ANI model overestimates the fluctuations of the strain field in zones of strong localization (Fig. 2). However, it is noted that for the same range of work-hardening exponents, this model slightly underestimates the equivalent average strain in the matrix. The evolutions ofε m eq , δ m (ε), δ m ⊥ (ε) and δ m eq (ε) as functions of the pore volume fraction show similar trends as those observed for strong and moderate nonlinearities. However, it is observed on Fig.  2b that the parallel MISF curves cross each others at f p ≃ 0.57 thus showing that the AFF-ANI procedure fails to reproduce qualitatively their evolution near the percolation limit. This is probably due to the fact that the AFF-ANI approach does not incorporate the strain field fluctuations which become significant near the percolation limit. For the stress field statistics, it is noted that the AFF-ANI procedure significantly overestimates the local equivalent stress in both the diffuse and localization regions in the matrix of the porous medium. On the other hand, it generally underestimates the parallel stress field fluctuations in both the diffuse area and the conical localization region. This result is confirmed by the probability density function which is globally shifted towards the left, and exhibits a much shorter tail. Moreover, it can be observed that AFF-ANI provides reasonably good estimates for the perpendicular stress fluctuations in the diffuse zone, with a slight overvaluation, while it overestimates this field more significantly in the localization region situated at the boundary of the pore.
These tendencies observed for a particular choice of m and f p turn out to be rather general. Indeed the results related to the evolution of the equivalent average stress with respect to the workhardening exponent show that AFF-ANI overestimates the reference equivalent average stress in the matrix for any nonlinearity. For the Measures of the Intraphase StresS Fluctuations (MISSF) over the matrix, it is also observed that AFF-ANI yields notably larger estimates than the reference δ m ⊥ (σ) for high nonlinearities while it slightly underestimates the reference δ m (σ). Accordingly, AFF-ANI overestimates δ m eq (σ) for large nonlinearity. Note in addition that these trends hold for any pore volume fraction and that the differences between the AFF-ANI estimates and the "exact" solution for δ m ⊥ (σ) and the variance of the stress field become important for f p ≥ 50% (i.e. near the percolation threshold) as the AFF-ANI model fails to reproduce qualitatively the trend of the reference evolutions for high range of pore volume fractions. At last, the AFF-ANI scheme provides a fair estimation of the incremental work density which, as mentioned in section 3.1, presents the advantage of combining the whole local tensorial stress and strain components without favoring any direction.

Specificities arising from the comparison of linearization scheme predictions at the local level
As shown above for the particular case of the AFF-ANI scheme, it is possible to extend the proposed methodology to all linearization schemes based on the conjecture that local fields and their statistics up to second-order in the defined LCC constitute reasonable approximations for the corresponding nonlinear quantities. Hence, by means of the proposed systematic methodology, the evaluation and the comparison at the local scale of the various linearization procedures considered in this study (SEC, VAR, AFF-ANI, AFF-ISOT and LS) were performed in an exhaustive way, thus allowing to explore some common points and some specificities related to their performances or limits in terms of assessing the local field fluctuations of porous media and reinforced composites under isochoric extension. As there are numerous results relative to each of the considered linearization methods [38], we report below only on the most salient features.

Porous media
• Although the AFF-ANI scheme reproduces qualitatively well the probability densities of the local fields in the matrix, it over-estimates the average of the stress over this phase (see Fig. 4) and hence the global response of the porous media as shown on Fig. 6a in [39]. In contrast, the LS procedure estimates well both the local fields (see Fig. 7b-i, 7b-ii, 8b-i, 8b-ii) and the global response [39]. This improved performance is due to the additional term in the polarization field, which does not significantly modify local strains but strongly shifts local and global stresses.
• In the case of porous media submitted to isochoric extension (i.e. a purely deviatoric loading), it is noted that linearization methods defining anisotropic phases for the LCC estimate the heterogeneity of the local fields better than those leading to an isotropic LCC. Indeed, for the local parallel stress fluctuations, isotropic linearizations (SEC, VAR and AFF-ISOT) predict additional localization regions -the green zones on Fig. 7c-i, 7d-i and 7e-i at the corners of the periodic cell -by comparison to the reference map. This explains why the probability density functions provided by isotropic linearizations are flat in opposition to the reference function and also the LS (Fig. 7b-ii) and AFF-ANI (Fig. 3c-ii) anisotropic linearization schemes which all exhibit a main narrow peak and a long tail. Moreover, the isotropic models highly underestimate the perpendicular stress fluctuations (see Fig. 8) at the diffuse zone in the matrix and hence yield to distribution functions with a peak shifted to the left in contrast to the "exact" results and the anisotropic schemes. The evolutions of the parallel δ m (σ) and perpendicular δ m ⊥ (σ) measures of the intraphase stress fluctuations as functions of the parameters m and f p (Fig. 5 and 6) prove again the relevance of anisotropic linearizations compared to isotropic ones. Indeed, the latter overestimate (resp. underestimate) the parallel (resp. perpendicular) measures of the intraphase stress fluctuations for all tested values of porosity and nonlinearity exponent. Moreover, we note some convergence problems of the fixed-point iterative scheme for the AFF-ISOT approach. • Although it has been known for a long time that the VAR procedure provides better estimates than SEC at the global level [7,39,42] as it incorporates the intraphase strain fluctuations through a reference strain equal to the second-order moment of the strain, it is noted here (Fig. 5 and 6) that both the VAR and SEC estimates for field fluctuations are almost superimposed for porous media submitted to a purely deviatoric extension. Appendix D provides an analytical proof of this result.

Rigidly reinforced composite
For the sake of brevity, only general features are reported here, without illustrations. More details can be found in [38].
• All models show a good ability to reproduce qualitatively (but often not quantitatively) the localization regions of the spacial distribution of the local fields, especially at the particle/matrix interface. Generally, the distribution functions only show one pronounced peak centered on the mean value of the estimated variable. Some times, they show a second less pronounced peak associated to a mean value in some diffuse zone. Finally, they exhibit a tail spread out towards the right demonstrating a strong localization of the studied variable in a few number of finite elements constituting the matrix mesh and located in the localization cone [38].
• Local strain fields resulting from linearization methods defining isotropic phases for the LCC as SEC, VAR and AFF-ISOT are identical. A proof of this result is given in Appendix D. Fig. 11a in [39] and Fig. 3.46 in [38] provide illustrations of this result.
• Similarly, the local strain fields resulting from tangent linearizations such as LS and AFF-ANI are identical in the matrix. Appendix E provides an analytical proof of this result and again illustrations can be found in [39] (see Fig. 11a) and [38] (Fig. 3.46).
• For the matrix, it is observed that the AFF-ISOT scheme highly underestimates local fluctu- ations of the stress field in the perpendicular direction σ ⊥ (x) and the measures of the perpendicular intraphase stress fluctuations σ m ⊥ over the matrix (Fig. 3.35b and 3.45c in [38]). This can be explained by the fact that this model is based on an empirical isotropic simplification of the original version of this model as it substitutes the secant shear modulus µ sct (ε m eq ) by the tangent one µ tgt (ε m eq ) in the perpendicular direction and that µ tgt (ε m eq ) << µ sct (ε m eq ). Accordingly it does not rest on rigorous physical arguments. On the other hand, it is worth noting that VAR scheme which also relies on an isotropic linearization, gives good estimates for the evolutions of the measures of the perpendicular intraphase stress fluctuations in the matrix as function of the work-hardening parameter m (again Fig. 3.45c in [38]) but overestimates the parallel fluctuations (Fig. 3.45b in [38]).

Conclusion and perspectives
In this work, we have developed a methodology for a systematic and accurate evaluation of the ability of a large class of nonlinear homogenization schemes to characterize local field fluctuations up to second order. In order to illustrate the relevance of this methodology, the local solution resulting from the special case of the AFF-ANI model was treated in an exhaustive way and was compared with the reference local solution. Further, field statistics in a rigidly reinforced composite and a porous material were also evaluated exhaustively for all models. However, due the significant number of results, only main features of the comparative analysis are reported in this paper. Though the conclusions resulting from this comparative analysis are not general since they depend on the type of material -here reinforced composite or porous material -on the type of the loading -axially symmetric loading without spherical part -and of the microstructure, the methodology is rather general and could be applied to other situations. For a porous material, it is shown that LS and AFF-ANI schemes are the models providing the best estimates for the local fields. Other above-mentioned conclusions are specific to the case of a reinforced composite. The proposed evaluation is limited by the fact that the stress and strain fields used to evaluate the effective behaviour have to correspond to those inside the LCC which is not always the case as for instance with the second-order methods. For this reason, the second-order procedures, which are known to provide the best estimates of both the effective behavior and the statistics of the local fields in most available results, e.g. [16,17,18,38], have not been evaluated at the local scale in this work. However the presented methodology can be partially extended (at least for the first and second-order moments of the local fields) to homogenization schemes not satisfying the aforementioned conditions by adding appropriate corrective terms to the first and second-order moment of the local fields as given by Idiart et al. [15,17,18]. Accordingly, new insight about the efficiency of the nonlinear homogenization techniques can be expected by applying this methodology to other more sophisticated and more recent models such as for the third variant of the second-order method and to other types of periodic materials with a unit cell composed of a larger number of inclusions and/or submitted to other types of loadings. Corresponding analysis are left for further investigations.
in the matrix of an isotropic LCC are solution of the following set of equations u kinematically admissible with prescribed displacement boundary conditions, ε = 1 2 (∇u + t ∇u), tr(ε) = 0, σ = −p(x)i + 2µ m K : ε, −p(x)i + 2µ m K : ε + τ m , isotropic secant linearizations (SEC and VAR), isotropic tangent linearizations (AFF-ISOT), div(σ) = 0, u = 0, or σ.n(x) = 0, at the matrix/particle interface for a rigidly reinforced composite, at the boundary of the pore for a porous material. (D.1) Theorem. For two-phase rigidly reinforced incompressible composites with constant shear modulus µ m and polarization tensor τ m in the matrix, the strain field solution to problem (D.1) associated to isotropic secant or tangent linearization schemes does not depend on µ m and τ m . For porous materials with incompressible matrix, the same result holds true for isotropic secant linearization schemes only. Accordingly, for two-phase rigidly reinforced incompressible composites (resp. porous materials with incompressible matrix), the isotropic linearization formulations SEC, VAR and AFF-ISOT (resp. SEC and VAR) lead to the same strain fields under a macroscopic strainε.
Proof. Let us first consider the case of isotropic secant linearization procedures for both two-phase rigidly reinforced incompressible composites and porous materials with incompressible matrix. For such a constitutive behavior, the LCC, which only depends on one scalar µ m , will be referred to as LCC(µ m ) in the following. Let us denote (u(x), ε(x), σ(x), p(x)) the solution of the local problem (D.1) associated with LCC(µ m ). Obviously, ∀ α > 0, (u(x), ε(x), ασ(x), αp(x)) is the solution of the local problem associated with LCC(αµ m ). Therefore, ε(x) does not depend on µ m .
Let us now consider the case of isotropic tangent linearizations for two-phase rigidly reinforced incompressible composites only. The corresponding LCC depends on two quantities µ m and τ m and is denoted LCC(µ m , τ m ). Let (u(x), ε(x), σ(x), p(x)) be the solution of the local problem (D.1) associated with LCC(µ m , τ m ). Since τ m is constant in the matrix, one has div(τ m ) = 0, so that (u(x), ε(x), σ(x), p(x)) is also solution of the problem LCC(µ m , 0), the kinematic field of which does not depend on µ m as seen above. Note that this proof associated with isotropic tangent linearizations can not be extended to porous materials because the BC σ.n(x) = 0 at the border of the pore is no longer satisfied. Note also that this property does neither hold for a compressible matrix, as the local strain depends in such a case on the ratio µ m /k m , i.e. on the Poisson's ratio, of the matrix.

Appendix E. Strain local field for tangent linearizations in rigidly reinforced composites
For tangent linearization methods defining anisotropic thermoelastic LCC as the AFF-ANI and LS approaches, the local stress and strain fields in a two-phase rigidly reinforced incompressible composite subjected to kinematic boundary conditions are solution to the following set of equations u kinematically admissible with prescribed boundary displacement conditions, u = 0 at the matrix/particle interface, is also the same for both procedures. Accordingly, the local strain fields solutions to problem (E.1) associated with the LCC derived from either the LS or the AFF-ANI approach are equal.