Franc ¸ois Willot
email: francois.willot@ensmp.fr

for computing the mechanical response of composites with accurate local fields

Keywords: FFT methods, Homogenization, Heterogeneous media, Linear elasticity, Computational mechanics, Spectral methods

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Fourier-based algorithms, or "FFT" methods for short, are an efficient approach for computing the mechanical response of composites. Initially restricted to linear-elastic media, FFT tools are nowadays employed to treat more involved problems, ranging from viscoplasticity [1] to crack propagation [START_REF] Li | A non-local fracture model for composite laminates and numerical simulations by using the fft method[END_REF]. In FFT methods, the microstructure is defined by 2D or 3D images and the local stress and strain tensors are computed along each pixel or "voxel" in the image. Coupled with automatic or semi-automatic image segmentation techniques [START_REF] Faessel | Segmentation of 3D microtomographic images of granular materials with the stochastic watershed[END_REF], this allows for the computation of the mechanical response of a material directly from experimental acquisitions, like focused ion beam or 3D microtomography techniques [START_REF] Willot | Microstructure-induced hotspots in the thermal and elastic responses of granular media[END_REF]. The latter often deliver images containing billions of voxels, for which FFT methods are efficient [START_REF] Dunant | A critical comparison of several numerical methods for computing effective properties of highly heterogeneous materials[END_REF][START_REF] Willot | Elastic and electrical behavior of some random multiscale highly-contrasted composites[END_REF]. This allows one to take into account representative volume elements of materials which are multiscale by nature such as concrete or mortar [START_REF] Escoda | Estimation of local stresses and elastic properties of a mortar sample by FFT computation of fields on a 3D image[END_REF]. From a practical viewpoint, the simplicity of FFT methods is attractive to researchers and engineers who need not be experts in the underlying numerical methods to use them. Nowadays, FFT tools are available not only as academic or free softwares [START_REF]Morphhom software[END_REF][START_REF]CraFT software[END_REF] but also as commercial ones [START_REF]GeoDict software[END_REF].

In the past years, progress has been made in the understanding of FFT algorithms. Vondřejc and co-workers have recently shown that the original method of Moulinec and Suquet [START_REF] Moulinec | A fast numerical method for computing the linear and non linear mechanical properties of the composites[END_REF] corresponds, under one technical assumption, to a particular choice of approximation space and optimization method [START_REF] Vondřejc | An FFT-based Galerkin method for homogenization of periodic media[END_REF] (see also [START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites[END_REF]). This property allows one to derive other FFT schemes that use standard optimization algorithms, such as the conjugate gradient method. In this regard, making use of variational formulations, efficient numerical methods have been proposed that combine FFTs with an underlying gradient descent algorithm [START_REF] Zeman | Accelerating a FFT-based solver for numerical homogenization of a periodic media by conjugate gradients[END_REF][START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF].

Different approximation spaces or discretization methods have also been proposed, where, contrary to the original scheme, the fields are not trigonometric polynomials anymore. Brisard and Dormieux introduced "energy-based" FFT schemes that rely on Galerkin approximations of Hashin and Shtrikman's variational principle [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites[END_REF] and derived modified Green operators consistent with the new formulation. They obtained improved convergence properties and local fields devoid of the spurious oscillations observed in the original FFT scheme [START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites[END_REF][START_REF] Brisard | Periodic homogenization using the lippmannschwinger formalism[END_REF]. In the context of conductivity, accurate local fields and improved convergence rates have also been obtained from modified Green operators based on finitedifferences [START_REF] Willot | Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF]. These results follow earlier works where continuum mechanics are expressed by centered [START_REF] Müller | Mathematical vs. experimental stress analysis of inhomogeneities in solids[END_REF][START_REF] Brown | Discrete fourier transforms and their application to stress-strain problems in composite mechanics: a convergence study[END_REF] or "forward and backward" finite differences [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwisedisordered porous media[END_REF].

This work focuses on the effect of discretization in FFT methods. It is organized as follows. We first recall the equations of elasticity in the continuum (Sec. 2). We give the Lippmann-Schwinger equations and the "direct" and "accelerated" FFT schemes in Sec. [START_REF] Faessel | Segmentation of 3D microtomographic images of granular materials with the stochastic watershed[END_REF]. In Sec. (4), a general formulation of the Green operator is derived that incorporates methods in [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwisedisordered porous media[END_REF], and a new discretization scheme is proposed. The accuracy of the local stress and strain fields are exam-ined in Sec. [START_REF] Dunant | A critical comparison of several numerical methods for computing effective properties of highly heterogeneous materials[END_REF] whereas the convergence rates of the various FFT methods are investigated in Sec. [START_REF] Willot | Elastic and electrical behavior of some random multiscale highly-contrasted composites[END_REF]. We conclude in Sec. [START_REF] Escoda | Estimation of local stresses and elastic properties of a mortar sample by FFT computation of fields on a 3D image[END_REF].

Microstructure and material elastic response

We are concerned with solving the equations of linear elasticity in a square or cubic domain Ω = [-1/2; 1/2] d in dimension d (d = 2 or 3):

σ i j (x) = C i j,kl (x)ε(x), ∂ i σ i j (x) ≡ 0, ε i j (x) = (1/2) ∂ i u j (x) + ∂ j u i (x) ,
(1) where ε(x) is the strain field, σ(x) the stress field, u(x) the displacement vector field, C(x) the local elasticity tensor and x is a point in Ω. Tensorial components refer to a system of Cartesian coordinates (e 1 ; e 2 ) in 2D and (e 1 ; e 2 ; e 3 ) in 3D. The material has an isotropic local elastic response that reads:

C i j,kl (x) = λ(x)δ i j δ kl + µ(x)(δ ik δ jl + δ il δ jk ), (2)
where δ is the Kronecker symbol and λ(x) and µ(x) are constant-per-phase Lamé's first and second coefficients. The local bulk modulus κ = λ+(2/d)µ and the elastic moduli take on values:

λ(x) = λ α , κ(x) = κ α , µ(x) = µ α ,
in phase α. For simplicity, we restrict ourselves to binary media and, by convention, α = 1 is the matrix and α = 2 the inclusions. Hereafter, we fix Poisson's ratios in each phase to ν 1 = ν 2 = 0.25 so that, in 3D and 2D [START_REF] Thorpe | Elastic moduli of two-dimensional composite continua with elliptic inclusions[END_REF], we have µ α /κ α = 0.6. The contrast of properties χ reads:

χ = κ 2 κ 1 = µ 2 µ 1 = λ 2 λ 1 , (3)
where 0 ≤ χ ≤ ∞. In the matrix, we also fix

κ 1 = 1 (d = 2 or 3), µ 1 = 0.6 (d = 2 or 3), λ 1 = 0.4 (d = 2), λ 1 = 0.6 (d = 3
), so that the local properties of the material are parametrized by one unique variable, the contrast of properties χ. In 3D, the Young modulus is E 1 = 3/2 in the matrix and E 2 = 3χ/2 in the inclusion. The medium is porous when χ = 0 and rigidly-reinforced when χ = ∞. Periodic boundary conditions are applied with the material subjected to an overall strain loading ε:

σ i j (x)n j (x) -# (x ∈ ∂Ω), ε kl (x) = ε kl , (4)
where n is the normal at the boundary ∂Ω of the domain Ω, oriented outward, -# denotes anti-periodicity and • denotes the spatial mean over Ω. The resulting effective elastic tensor C is computed from:

σ i j (x) = C i j,kl ε kl . (5)

Lippmann-Schwinger equation and FFT methods

Fourier methods are by principle based on the Lippmann-Schwinger equations. The latter follow from (1) and ( 4) as [START_REF] Milton | The Theory of Composites[END_REF]:

τ i j (x) = σ i j (x) -C 0 i j,kl ε kl (x), ε i j (x) = ε i j - x d d x G i j,kl (x -x)τ kl (x ), (6)
where we have introduced a homogeneous "reference" elasticity tensor C 0 and its associated polarization field τ and Green operator G. In the above we assume G = 0 so that ε = ε(x) holds. The Green operator has, in the Fourier domain, the closed form [START_REF] Kanaun | Self-consistent methods for composites[END_REF]:

G i j,kl (q) = q i q m C 0 m j,kn q n -1 q l sym , (7)
where q 0 are the Fourier wave vectors and the subscript sym indicates minor symmetrization with respect to the variables (i, j) and (k, l). Hereafter, we assume that C 0 is a symmetric, positive-definite, isotropic tensor defined by its bulk (κ 0 ) and shear (µ 0 ) moduli, or Lamé coefficient (λ 0 ). Accordingly, when q 0, the Green operator G is also symmetric definite and we have:

G i j,kl (q) = 1 µ 0       q i q l |q| 2 δ jk sym - λ 0 + µ 0 λ 0 + 2µ 0 q i q j q k q l |q| 4       . (8)
The "direct scheme" [START_REF] Moulinec | A fast numerical method for computing the linear and non linear mechanical properties of the composites[END_REF] consists in applying Eqs. ( 6) iteratively as:

ε k=0 ≡ ε, ε k+1 = ε -G * (σ -C 0 : ε k ) (k ≥ 0). (9)
In Moulinec and Suquet's method, the convolution product ( * ) above is computed as an algebraic product in the Fourier domain, making use of [START_REF]Morphhom software[END_REF]. Discrete Fourier transforms are used to switch between the space Ω and Fourier domain F . This amounts to representing the strain field as a trigonometric polynomial [START_REF] Vondřejc | An FFT-based Galerkin method for homogenization of periodic media[END_REF] of the form:

ε(x) = 1 L d q∈F ε(q) e iq•x , (10)
where ε(q) is the discrete Fourier transform of ε(x). Accordingly:

ε(q) = x∈Ω ε(x) e -iq•x . (11)
Similar forms are used for the stress and displacement fields. In practice, the domain Ω is discretized on a square or cubic grid of L d voxels and the operator G(q) in ( 7) is evaluated along equispaced Fourier modes:

q i = 2πm i L , m i = 1 -(L/2), ..., L/2, (L even), -(L -1)/2, ..., (L -1)/2, (L odd). ( 12
)
As noted in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF], when L is even, the relation:

G(q) * = G(-q), (13)
where G(q) * is the complex conjugate of G(q), is not verified when one of the components of q is equal to the highest frequency q i = π (i.e. m i = L/2). As a consequence, the backward Fourier transform of G(q)τ(q) used to compute the strain field has non-zero imaginary part even if τ(x) is purely real. To fix this problem, we follow [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] and set:

G(q) = C 0 -1 , if q i = π for some i. (14)
This choice enforces σ(q) = 0 at the concerned Fourier modes. In doing so, the strain field ε in (10) is not strictly-speaking irrotational, because of the lack of constraint at high Fourier modes for the strain field. We briefly mention another option that we explored in this work. It consists in forcing the symmetry by replacing G(q) with:

G(q) * + G(-q) 2 , ( 15
)
when one of the components of q equals π, which enforces ε(q) = 0 at the highest modes. Numerical experiments indicate that the choice for G(q) at the highest frequencies has little influence on the convergence rate, except at small resolution. When L < 128 pixels, faster convergence was achieved with the choice G(q) = C 0 -1 . Furthermore, in the 2D example studied in this work, the choice G(q) = C 0 -1 led to smaller oscillations, consistently with observations in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. The use of (15) was therefore not pursued further. We emphasize that, when L is odd, this discrepancy disappears and no special treatment is needed.

Refined FFT algorithms have been introduced to overcome the slow convergence rate of the direct scheme, observed for highly-contrasted composites, most notably the "accelerated scheme" [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] and "augmented Lagrangian" [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] methods. In this work, we use the extension of the accelerated scheme to elasticity [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF][START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF]:

ε k+1 = ε k + 2 C + C 0 -1 : C 0 : ε -ε k -G * C : ε k -C 0 : ε k . ( 16
)
The convergence rates of the accelerated and direct schemes depend on the choice of the reference tensor C 0 . For the accelerated scheme the optimal choice is [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF]:

κ 0 = √ κ 1 κ 2 , µ 0 = µ 1 µ 2 . ( 17
)
For the direct scheme, upper bounds on the eigenvalues of the Green operator suggest the choice [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF]:

κ 0 = β κ 1 + κ 2 , µ 0 = β µ 1 + µ 2 , (18)
with β = 1/2.

Discretization and approximation space

In this section, we derive the expression of a modified Green operator G that replaces G defined in [START_REF]Morphhom software[END_REF] and [START_REF] Vondřejc | An FFT-based Galerkin method for homogenization of periodic media[END_REF]. We give it in a form that includes previouslyproposed modified operators [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwisedisordered porous media[END_REF] and also introduce a new one.

Two dimensions

In the following, we assume that the strain and stress fields are defined on a grid of points in 2D, one per pixel. Eqs. [START_REF]GeoDict software[END_REF] and [START_REF] Moulinec | A fast numerical method for computing the linear and non linear mechanical properties of the composites[END_REF] are used to apply discrete Fourier transforms, but we do not postulate a representation in the continuum anymore. The equilibrium and strain admissibility conditions (1) are approximated by means of finite differences on which we apply the discrete transforms [START_REF]GeoDict software[END_REF] and [START_REF] Moulinec | A fast numerical method for computing the linear and non linear mechanical properties of the composites[END_REF]. In [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwisedisordered porous media[END_REF], this results in the following form:

k * i (q)σ i j (q) = 0, ε i j (q) = (1/2) k i (q)u j (q) + k j (q)u i (q) , (19)
where k and k * represent "discrete" gradient and divergence operators, respectively. In the centered scheme, one takes k equal to:

k C i (q) = i sin (q i ) , (20)
whereas in scheme [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwisedisordered porous media[END_REF], one chooses, for k:

k W i (q) = e iq i -1. ( 21
)
These expressions correspond, respectively, to the centered scheme:

∂ j σ i j (x) ≈ σ i j (x + e j ) -σ i j (x -e j ) 2 , ∂ j u i (x) ≈ u i (x + e j ) -u i (x -e j ) 2 , (22)
and to the forward-and-backward difference scheme:

∂ j σ i j (x) ≈ σ i j (x) -σ i j (x -e j ), ∂ j u i (x) ≈ u i (x + e j ) -u i (x). (23)
Using ( 19), the resulting Green operator reads:

G i j,kl (q) = k i (q) k m (q)C 0 m j,kn k * n (q) -1 k * l (q) sym , (24)
which is homogeneous in k, so that, with r i = k i /|k|:

G i j,kl (q) = (λ 0 + 2µ 0 ) r i r * l δ jk sym -λ 0 Re r i r * j Re r k r * l -µ 0 r i r j (r k r l ) * µ 0 2(λ 0 + µ 0 ) -λ 0 r 2 1 + r 2 2 2 , (25)
where Re(•) denotes the real part of the enclosed quantity. The denominator on the right-hand side is strictly positive due to the triangle inequality. Hereafter, the operator G is denoted by G C when k = k C and by G W when k = k W . Note that the operator G in ( 8) is recovered from (25) by setting k i (q) = iq i . Now the Green operator G is complex and follows the minor and major symmetries:

G i j,kl (q) = G ji,kl (q) = G i j,lk (q), G i j,kl (q) = G kl,i j (q) * . (26)
Furthermore we have:

G C,W (q) * = G C,W (-q), (27)
for the schemes ( 22) and ( 23), including when q i = π and when L is even, because k C,W i is real at the frequency q i = π. Therefore, the fix ( 14) in Sec. ( 3) is not necessary. However, a problem of a different nature arises when using the centered scheme [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwisedisordered porous media[END_REF] when L is even. Eq. ( 25) does not define the Green operator G C at the three frequencies q = (π; 0), (0; π) and (π; π), for which k C (q) = 0. This is because the second equation in [START_REF] Milton | The Theory of Composites[END_REF] has in general non-unique solutions for the displacement field u(x). Indeed, when L is even, the displacement is defined up to a linear combination of 2-voxels periodic fields. They are given by the following 8 independent fields:

v 1 m (x) = δ mn , v 2 m (x) = (-1) i δ mn , v 3 m (x) = (-1) j δ mn , v 4 m (x) = (-1) i+ j δ mn , x = i L ; j L , i, j = 1 -L 2 , ..., L 2 , n = 1, 2. ( 28
)
The operator G C remains finite when q approaches one of the modes (π; 0), (0; π) or (π; π), but can not be continuously extended at these modes. To fix this problem, we set, for the centered scheme:

G C i j,kl (q) = 0, if L is even and (q i = 0 or q i = π) for all i, (29)
which enforces ε(q) = 0 at the highest frequencies. The strain field ε(x) is accordingly admissible, and the stress field σ(x) is divergence-free, in the sense of [START_REF] Milton | The Theory of Composites[END_REF]. We explored the alternate choice G C (q) = C 0 -1 in [START_REF] Willot | Elastic behavior of composites containing boolean random sets of inhomogeneities[END_REF]. Almost identical convergence rates and oscillations were observed for the two options, however the choice

G C (q) = C 0 -1
does not produce an irrotational strain field and is not considered further. We emphasize that no special treatment is required for the operator G W at high frequencies since k W 0 when q 0. By substituting G with G C or G W in ( 9) and ( 16), we derive "direct" and "accelerated" schemes that solve [START_REF] Milton | The Theory of Composites[END_REF] or [START_REF] Kanaun | Self-consistent methods for composites[END_REF]. In the limit of very fine resolution, we have G C,W (q) ≈ G(q) when q → 0, which guarantees that the strain and stress fields do not depend on the employed discretization. This property holds for any choice of k such as k ∼ iq when q → 0.

On the one hand, derivatives are estimated more locally in the forward-andbackward scheme [START_REF] Thorpe | Elastic moduli of two-dimensional composite continua with elliptic inclusions[END_REF] than in the centered scheme [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwisedisordered porous media[END_REF], which is important along interfaces. On the other hand, the forward-and-backward scheme does not treat symmetrically the two angle bisectors e 1 + e 2 and e 1e 2 [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwisedisordered porous media[END_REF]. In a domain containing a single centered disc, the scheme produces fields that break the axial symmetries of the problem. In fact, the discretization ( 23) is actually one of four possible choices, all of them breaking the symmetries. Attempts to force the symmetry by averaging over the four Green operators or over the fields themselves, as proposed in [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwisedisordered porous media[END_REF], are not explored in this work. The fomer method indeed leads to less accurate "diffuse" local fields. The latter necessitates to run four different computations, in 2D, instead of one, which is cumbersome.

In the rest of this section, we derive a discrete scheme in 2D different from ( 22) and [START_REF] Kanaun | Self-consistent methods for composites[END_REF]. In this scheme, the displacement field is evaluated at the 4 corners of the pixels and the strain and stress fields are evaluated at the centers of the pixels. We first express these fields in the 45 o -rotated basis:

f 1 = e 1 + e 2 √ 2 , f 2 = e 2 -e 1 √ 2 , (30)
by:

u i = R iI u I , ε i j = R iI ε I J R J j , σ i j = R iI σ I J R J j , R iJ = 1 -2δ i1 δ J2 √ 2 , (31)
where uppercase indices refer to components in the rotated grid. We discretize (1) in the rotated basis by the centered differences (see Fig. 1):

σ I J (x) = C I J,KL (x) ε KL (x) , (32a) σ I1 (x) -σ I1 x - √ 2f 1 + σ I2 x + f 2 -f 1 √ 2 -σ I2 x -f 1 +f 2 √ 2 = 0, (32b)
ε KL (x) = 1 2 √ 2 u K x + f L √ 2 -u K x -f L √ 2 + u L x + f K √ 2 -u L x -f K √ 2 (32c)
where x lie at the centers of the pixels and x ± f I / √ 2 lie at the corners. Expressing back [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF] in the original Cartesian grid (e 1 ; e 2 ) and applying the backward discrete Fourier transform [START_REF]GeoDict software[END_REF] we arrive again at [START_REF] Brown | Discrete fourier transforms and their application to stress-strain problems in composite mechanics: a convergence study[END_REF] with the following expression for k:

k R i (q) = i 2 tan q i 2 1 + e iq 1 1 + e iq 2 . ( 33
)
We denote by G R the corresponding Green operator, derived by substituting k = k R in [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF]. The operator G R is real and also verifies:

G R (q) = G R (-q). (34)
However, when L is even, k R = 0 when q = (π; π) and G R is not defined by [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] at this frequency. Again, this is because (32c) gives the displacement field up to linear combinations of the 4 independent fields v 1,4 1,2 (see 28). Accordingly we set:

G R (q) = 0, when L is even, d = 2, q 1 = q 2 = π, (35)
which enforces strain compatibility and stress equilibrium, in the sense of [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF].

Figure 1: A pixel with edges parallel to the Cartesian axis (e 1 ; e 2 ). Superimposed: 45 o -rotated basis (f 1 ; f 2 ). The strain and stress fields are evaluated at the pixel center x (square). The displacement and the divergence of the stress field lie along the pixel corners (disks).

Three dimensions

We follow the same methodology in 3D. The equilibrium and strain admissibility conditions [START_REF] Brown | Discrete fourier transforms and their application to stress-strain problems in composite mechanics: a convergence study[END_REF] are unchanged, as well as the expression for the vectors k C,W in ( 20) and ( 21) resulting from ( 22) and [START_REF] Kanaun | Self-consistent methods for composites[END_REF]. In 3D, we also extend (33) as:

k R i (q) = i 4 tan q i 2 1 + e iq 1 1 + e iq 2 1 + e iq 3 , (36)
for the rotated scheme. The strain and stress fields are now evaluated at the centers of the voxels and the displacement field at their corners. Derivatives of the displacement are estimated by differences at opposite corners. For the strain components ε 11 and ε 12 :

ε 11 (x) ≈ 1 4 m=±1 n=±1 u 1 x + e 1 + me 2 + ne 3 2 -u 1 x - e 1 -me 2 -ne 3 2 ,(37a) ε 12 (x) ≈ 1 8 m=±1 n=±1 u 2 x + e 1 + me 2 + ne 3 2 -u 2 x - e 1 -me 2 -ne 3 2 + u 1 x + e 2 + me 1 + ne 3 2 -u 1 x - e 2 -me 1 -ne 3 2 , ( 37b
)
where x lie at the center of a voxel. The expression for the strain component ε 22 (resp. ε 33 ) is obtained after exchanging the indicia 1 and 2 (resp. 1 and 3) in (37a).

The component ε 23 (resp. ε 13 ) is derived from (37b) by exchanging the indicia 3 and 1 (resp. 3 and 2). Stress divergence is discretized in a similar manner. Its first component reads:

∂ i σ i1 (x) ≈ m=±1 n=±1 σ 11 x + e 1 + me 2 + ne 3 2 -σ 11 x - e 1 -me 2 -ne 3 2 +σ 12 x + e 2 + me 1 + ne 3 2 -σ 12 x - e 2 -me 1 -ne 3 2 +σ 13 x + e 3 + me 1 + ne 2 2 -σ 13 x - e 3 -me 1 -ne 2 2 , (38)
where x lie at one of the edges of a voxel. The components ∂ i σ i2 and ∂ i σ i3 are obtained from (38) by circular permutations of the indicia 1, 2 and 3. We note that (37) and ( 38) are the natural generalization of (32) to d = 3, expressed in the Cartesian basis (e 1 ; e 2 ; e 3 ). In 3D, Eq. ( 24) yields, for the Green operator:

G i j,kl (q) = λ 0 + 2µ 0 r i r * l δ jk sym + λ 0 r i r * l s jk sym -Re r i r * j Re r k r * l -µ 0 r i r j r * k r * l µ 0 2(λ 0 + µ 0 ) -λ 0 r 2 1 + r 2 2 + r 2 3 2
(39) where again r i = k i /|k| and s is the symmetric second-order tensor:

s j j = 4 Im(r i r * k ) 2 , s ji = -4 Im(r k r * j ) Im(r k r * i ) [instead of s jk = -4...], i j k i, (40)
with Im(•) the imaginary part of the enclosed complex quantity. Like in 2D, the operator G follows minor and major symmetries [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF].

Again, the operators G C , G W and G R are derived using the expression for G in (39) with k = k C , k W and k R , respectively. The symmetries [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF] and (34) are verified in 3D as well. But again, a special treatment is needed for G C and G R when L is even, at the modes q for which k C,R (q) = 0. Like in 2D, the displacement is undefined at these frequencies and the Fourier coefficients of the strain field are zero and so we set G C,R = 0 at these frequencies. More precisely, we apply (29) when d = 3 and, for the rotated scheme:

G R (q) = 0,
if L is even, d = 3 and q i = q j = π with i j.

The operators G C , G W and G R are, in 2D and 3D, periodic functions where, contrary to G, high frequencies are cut. Accordingly, we expect faster convergence rates for schemes using operators derived from finite differences and more exact local fields, as was previously observed in the conductivity problem [START_REF] Willot | Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF].

We also expect higher accuracy for the local fields when employing G R rather than the other discrete operators G W and G C . First, the operator G R is based on centered differences which are more precise than forward and backward differences, used in G W . Second, derivatives are evaluated more locally when using G R rather than G C . Indeed, the latter are computed at points separated by 2 voxels for

G C instead of √ 2 (in 2D) or √ 3 voxels (in 3D
) for G R . The above considerations guided the choice for the discretization schemes (37) and (38), leading to k R and G R . Clearly, many other choices are possible, and Eq. (39) gives a general class of Green operators based on finite-differences. The latter depend on the choice for the complex vector k. However, a systematic investigation of such discrete schemes is beyond the present study.

In the the rest of this study, we estimate the accuracy of the local fields and of the effective properties predicted by the various schemes, as well as their convergence rates. We denote by DS and AS the direct and accelerated schemes defined by ( 9) and ( 16) respectively, when G is used. We denote by DS C , DS W , DS R and AS C , AS W and AS R , the same algorithms obtained by substituting G with G C,W,R respectively. We emphasize that, for a given Green operator, the direct and accelerated schemes produce the same strain and stress fields, up to round-off errors.

Local strain and stress fields accuracy

Two-dimensions

Hereafter we consider the 2D 'four-cell' microstructure, where the periodic domain Ω is divided into 4 identical squares of surface fraction 25%. Its nontrivial solution with singular fields at the corners makes it a good benchmark for numerical schemes. Furthermore, the microstructure is discretized exactly at any resolution, provided L is even. In the following, we make use of a simplified version of the four-cell microstructure made of a single quasi-rigid square inclusion embedded in a matrix (Fig. 2). We set the contrast to χ = 10 3 . The material is subjected to the macroscopic strain loading:

ε i j = 1 2 δ i1 δ j2 + δ j1 δ i2 . (42)
We determine the strain and stress fields predicted by FFT schemes when using the Green operators G or G C,W,R . The fields are computed using the accelerated scheme ( 16) at discretizations L = 512, 1024 and 2048. Iterations are stopped when the strain and stress fields maximum variation over two iterations in any pixel is less than 2 10 -13 . These variations are the effect of round-off errors in double precision floating point numbers. These computations allow us to compare the effect of the discretization, independently of the algorithm used for convergence. We focus on the stress component σ 12 (x) parallel to the applied loading in a small region [-0.04; 0.04] 2 around the corner of the inclusion (Fig. 3). At low resolution L = 512, numerical methods predict values as large as 10.1 in a few pixels, because of the singularity of the stress field at the corner. To highlight the field patterns, we threshold out the values above 3.5, which amount to 0.24% of the pixels. Using the same color scale for all images, the smallest stress value, equal to 1.5, is shown in dark blue whereas the highest, equal to 3.5, is in dark red. Green, yellow and orange lie in-between.

As expected, in the limit of very fine resolution, all methods tend to the same local stress field, as shown by the similar field maps obtained at resolution L = 2048. However, use of the Green operator G leads to spurious oscillations along the interfaces of the inclusion, up to resolutions as big as 2048 2 pixels, a side-effect noticed in [START_REF] Willot | Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF] in conductivity. The oscillations do not disappear after computing local averages of the fields (not shown).

Strong oscillations are produced by schemes using G C as well, not only in the quasi-rigid inclusion, but also in the matrix. We observe checkerboard patterns in the former, and vertical and horizontal alignments in the latter, at resolution 1024 2 . These oscillations are greatly reduced by the use of G W . Still, due to the non-symmetric nature of G W , the stress is not correctly estimated along a line of width 1 pixel oriented upward from the inclusion corner. Similar patterns are observed, in other directions, along the three other corners of the inclusion (not shown). These issues are solved when using G R which produces a stress field that respects the symmetries of the problem. Furthermore, use of G R greatly reduces oscillations compared to G and G C .

Three-dimensions

In this section, we consider a 3D material analogous to the four-cell microstructure in 2D. We divide the periodic domain into 8 identical cubes of volume fraction 12.5%. One is the inclusion, the other 7 are the matrix. To highlight the symmetries of the problem, we assume the inclusion is centered in the domain Ω and contained in the region [-1/4; 1/4] 3 . Again, we apply a macroscopic strain loading of the form (42). The inclusion is quasi-rigid compared to the matrix with contrast of properties χ = 10 3 . We compute the strain and stress fields predicted by each Green operator using the accelerated scheme. As in Sec. (5.1) we let the iterations converge up to round-off errors in double precision. 4). The section is a cut parallel to one of the faces of the inclusion, normal to e 3 , of equation x 3 = -0.2461. The section intersects the inclusion, but is very close to the interface with the matrix. Again, to highlight the field patterns, we threshold out values of the field greater than 8.5, this time less than 0.04% of the voxels, and represent all field using the same color scale.

At high resolution L = 1024, the fields resulting from the use of G and G C,W,R are close to one another. However, stress patterns near the corners of the inclusion are less pronounced with G than with the other methods. At smaller resolutions L = 256 and L = 512, the stress fields predicted by G are notably different from the others, suggesting slower size-convergence with this operator. Furthermore, the field maps computed at resolution L = 512 confirms the results obtained in 2D: strong oscillations are observed inside the inclusion when using G and G C . The two methods produce artificial patterns directed vertically and horizontally, close to the interface. Conversely, the fields produced by G W and G R have the smallest oscillations, but that of G W are not symmetric. When L = 256, indeed, the stress field near the top-left corner of the inclusion stands out from that in the other corners. This effect only slowly disappears when L is increased. The solution resulting from the use of G R does not suffer from this problem. As in 2D, it produces symmetric fields. Furthermore, the latter are close to one another at all resolutions and contain almost no oscillations. L = 256 The section is parallel to one of the faces of the inclusion and close to the interface with the matrix.

L = 512 L = 1024 L = 2048 G G C G W G R
L = 512 L = 1024 G G C G W G R

Periodic array of spheres

Contrary to the previous sections, we now consider a microstructure without singularities (edges or corners) and focus on the effect of the Green operator discretization on the effective elastic properties. In the rest of this section, the elementary domain Ω contains one spherical inclusion of volume fraction 20%, so that the material is a periodic array of spheres. The spheres are very soft with contrast of properties χ = 10 -4 . We compute the effective elastic modulus C 11,11 produced by either G or G C,W,R at increasing resolutions L = 32, 64, 128, 256 and 512. Again, we use the accelerated scheme and iterations are stopped when the stress field maximum variation over two iterations in any pixel is less than 2 10 -10 . Results are shown in Fig. [START_REF] Dunant | A critical comparison of several numerical methods for computing effective properties of highly heterogeneous materials[END_REF] and are compared with the analytical estimate in [START_REF] Cohen | Simple algebraic approximations for the effective elastic moduli of cubic arrays of spheres[END_REF]. When the resolution increases, the effective elastic modulus C 11,11 increases up to a limit value that we estimate to about 1.208 ± 0.001, for all schemes. As observed in other studies [START_REF] Dunant | A critical comparison of several numerical methods for computing effective properties of highly heterogeneous materials[END_REF], very large systems are needed to compute this estimate at a high precision.

This is especially true of the Green operator G which has the slowest convergence with respect to the system size. At fixed resolution, the error on the predictions given by G is about 2 times larger than the one provided by G R , which, among all methods, gives the best estimate. The operators G C and G W stand inbetween. This is another indication of the benefits of the operator G R .

Convergence rate

Convergence rate with respect to stress equilibrium

In this section, we estimate the rates of convergence of the direct and accelerated schemes DS, DS C,W,R , AS and AS C,W,R , that use the various Green operators. All schemes enforce stress equilibrium at convergence only, therefore we follow [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF] and consider a criterion based on the L 2 -norm:

η = ||div(σ)|| 2 || σ || = 1 || σ || 1 |Ω| Ω d d x |div(σ)| 2 = 1 || σ || q |k(q) • σ(q)| 2 , (43)
where η 1 is the precision and the normalizing factor || σ || is the Frobenius norm: In (43) we set k = k C,W,R for the schemes using G C,W,R and k = iq when using the Green operator G, so that k • σ(q) is the divergence of the stress field in the Fourier domain, estimated according to the various discretization schemes.

|| σ || 2 = i, j σ i j (x) 2 .
We now estimate the convergence rates on a random microstructure. In the following, the domain Ω is a (periodized) Boolean model of spheres of resolution L = 64 and volume fraction 17%, below the percolation threshold of the spheres -of about 29% [START_REF] Willot | Elastic behavior of composites containing boolean random sets of inhomogeneities[END_REF]. To obtain meaningful comparisons, we use the same randomly-generated microstructure for all schemes. This particular configuration contains 743 spheres of diameter 5 voxels, about 13 times smaller than the size of Ω.

Taking ν 0 = ν 1 = ν 2 = 0.25 for the reference Poisson ratio, we compute numerically the number of iterations N(E 0 ) required to reach the precision η ≤ 10 -8 , for varying reference Young moduli E 0 , in the range 0 < E 0 < 1. We consider the Boolean model of spheres with contrast χ = 10 -5 and the various accelerated schemes AS and AS C,W,R (Fig. 6). Within the range 0 < E 0 0.03, the number of iterations N(E 0 ) is about the same for all accelerated schemes. When E 0 > 0.03, however, N(E 0 ) is a strongly increasing function of E 0 for scheme AS, contrary to the other schemes AS C,W,R . For the latter, N(E 0 ) decreases with E 0 up to a local minimum, beyond which variations are much less sensitive to E 0 . One unique local minimum around E 0 ≈ 0.09 is found for scheme AS R , whereas the schemes AS W and AS C exhibit two local minima.

The effect of the Poisson ratio is also investigated numerically. We let ν 0 = 0.25 ± 0.01 and 0.25 ± 0.05 for various values of E 0 with χ = 10 -5 and observe a strong increase of the number of iterations N(E 0 ), for the schemes AS and AS C,W,R . The same behavior is observed for the direct scheme DS and DS C,W,R with χ = 10 -2 . Therefore, in the following, we fix the Poisson ratio to ν 0 = 0.25 for the reference tensor, for all schemes and all contrast of properties χ. This leaves one parameter, E 0 , to optimize on. We use the gradient descent method to determine a local minimum of N(E 0 ) for arbitrary contrast and scheme DS, DS C,W,R , AS and AS C,W,R . As above, N(E 0 ) is the number of iterations necessary to reach η ≤ 10 -8 . We choose E 0 = 0.51(E 1 + E 2 ) for schemes DS, DS C,W,R and E 0 = √ E 1 E 2 for schemes AS and AS C,W,R as initial guess for E 0 , suggested by ( 18) and [START_REF] Willot | Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF]. At each step, we determine if E 0 is to be increased or decreased, by comparing N(E 0 ) with N(E 0 +δE 0 ) where δE 0 = 0.01E 0 . It frequently happens that N(E 0 ) = N(E 0 + δE 0 ). In that case, we compare the values of the precision η after N(E 0 ) iterations and follow the direction that minimizes η. Iterations are stopped whenever N(E 0 ) is unchanged after two descent steps.

The gradient descent method determines a local minimum rather than the global minimum, which is sub-optimal. To check the validity of the results, further numerical investigations are carried out for χ = 10 -2 , 10 2 and schemes DS and DS C,W,R . The method predicts the global minimum in these cases. This also holds for schemes AS and AS W,R with χ = 10 -5 , 10 5 , but not for scheme AS C with χ = 10 -5 . However, in this case the number of iterations N(E 0 ) are very similar at the two local minima, as shown in Fig. [START_REF] Willot | Elastic and electrical behavior of some random multiscale highly-contrasted composites[END_REF]. In the following, the results given by the gradient descent method are used as-is.

Results for the optimal reference E 0 are indicated in Fig. [START_REF] Escoda | Estimation of local stresses and elastic properties of a mortar sample by FFT computation of fields on a 3D image[END_REF]. For the direct scheme, the optimal reference follows (18) with 0.5003 ≤ β ≤ 0.509, independently of the Green operator used. Values of β smaller than 1/2 lead to nonconverging schemes. For the accelerated schemes, the situation is less simple, and differs depending on the Green operator in use. For the scheme AS with Green operator G, the choice ( 17) is optimal except in the region χ ≤ 10 -3 where the value of E 0 tend to a small constant of about 0.01. Similar behavior is found for the schemes AS C,R for which:

E 0 ≈ E 0 1 + √ χ, (44)
with E 0 1 = 0.07 for AS C and E 0 1 = 0.12 for AS R . Similar behavior has been observed numerically in [START_REF] Willot | Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF], in the context of conductivity. For the scheme AS W that uses G W , the optimal choice for E 0 follows the same pattern as above with E 0 1 = 0.7 except in the region 10 -4 ≤ χ ≤ 10 -1 . This behavior is an effect of the presence of two local minima, similar to that shown in Fig. [START_REF] Willot | Elastic and electrical behavior of some random multiscale highly-contrasted composites[END_REF] for χ = 10 -5 .

Convergence rates, computed with optimized reference, are represented in Fig. [START_REF]Morphhom software[END_REF] as a function of the contrast, in log-log scale. Results for χ = 0 (strictly porous media) have been included in the same graph (left point). As is wellestablished [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF][START_REF] Moulinec | Comparison of FFT-based methods for computing the response of composites with highly contrasted mechanical properties[END_REF], the number of iterations in the direct scheme DS scales as χ when χ 1 and 1/χ when χ 1. For the accelerated scheme AS, the number of iterations is smaller and follows √ χ when χ 1 and 1/ √ χ when χ 1, with one exception. At very high contrast of properties χ < 10 -6 , including at χ = 0, convergence is reached after a finite number of iterations, about 1, 300. This particular behavior is presumably sensitive to the value choosen for the requested precision η = 10 -8 . When χ < 1, the schemes DS C,W,R and AS C,W,R , that use G C,W,R , converge after a number of iterations not exceeding 430. As shown in Fig. [START_REF]Morphhom software[END_REF], the number of iterations is nearly constant in the range 0 ≤ χ ≤ 10 -5 . As expected, the accelerated schemes AS C,W,R are faster than the direct schemes DS C,W,R , with scheme 20 AS R proving the fastest. With this scheme and when χ < 1, the number of iterations is at most 168. Again, these results are qualitatively similar with that given in [START_REF] Willot | Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF] in the context of conductivity.

For rigidly-reinforced media (χ > 1), the number of iterations of schemes DS C,W,R and AS C,W,R follow the same powerlaw behaviors, with respect to χ, as that of DS or AS. In all considered schemes, the number of iterations continuously increases with the contrast. Differences are observed between the various accelerated schemes AS and AS C,W,R , with AS R the fastest. The use of Green operators associated to the problem for the strain field, as undertaken here, results in convergence properties that are worse in the region χ > 1 than when χ < 1. In this respect, benefits are to be expected from the use of dual Green operators [START_REF] Willot | Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF], associated to the problem for the stress fields.

Convergence rate with respect to the effective elastic moduli

In this section, we focus on the accelerated schemes AS and AS C,W,R and examine the rate of convergence of the various schemes with respect to the effective elastic moduli. We consider the same Boolean microstructure as given in Sec. ( 6 but discretized on higher resolution grids of 256 3 and 512 3 voxels. The volume fraction of the spheres are respectively 16.82% and 16.85%. For simplicity, the contrast of properties take on two values χ = 10 -4 and χ = 10 4 , so that the spheres are quasi-porous or quasi-rigid. We perform iterations of the schemes AS and AS C,W,R using the optimized reference moduli found in the previous section, on the lower resolution grid. We apply the macroscopic strain loading:

ε i j = δ i1 δ j1 .
At each iteration and for each scheme, we compute the elastic modulus C 11,11 , derived from the mean σ 11 of the stress component σ 11 . The convergence rate toward the elastic modulus is represented in Fig. [START_REF]CraFT software[END_REF] for quasi-porous spheres with L = 512 and in Fig. [START_REF]GeoDict software[END_REF] for quasi-rigid spheres with L = 256. In Fig. [START_REF]CraFT software[END_REF], for the sake of clarity, the elastic moduli are represented by symbols once every 5 iterations, except for the first five iterations of the scheme AS R which are all represented. Dotted lines are guide to the eyes. In the porous case, much better convergence is obtained with scheme AS R than with schemes AS and AS C,W , as shown in Fig. strong oscillations that are much reduced with AS R . After about 7 iterations, the estimate given by AS R is valid to a relative precision of 10 -2 . To achieve the same precision, more than 50 iterations are needed for schemes AS and AS C,W .

The situation is notably different for quasi-rigid spheres (Fig. 10). For all schemes, a much higher number of iterations is required to determine the elastic modulus C 11,11 with a precision of 10 -2 . The slower convergence rate follows that observed in Sec. [START_REF] Willot | Elastic and electrical behavior of some random multiscale highly-contrasted composites[END_REF], where convergence is much poorer for χ > 1 than for χ < 1, and where AS R is less advantageous compared to the other schemes . Nevertheless, in this case also, as shown in Fig. [START_REF]GeoDict software[END_REF], smaller oscillations are observed in the estimate for C 11,11 when using AS R rather than schemes AS or AS C,W .

Conclusion

In this work, a novel discretization method has been proposed in 2D and 3D for use in Fourier-based schemes. The core of the proposed scheme is a simple modification of the Green operator in the Fourier domain. The results obtained confirm those achieved in the context of conductivity [START_REF] Willot | Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF]. Compared to schemes using trigonometric polynomials as approximation space, or to other finite-differences methods, superior convergence rates have been observed in terms of local stress equilibrium, but also in terms of effective properties. More importantly, the solution for the local fields, predicted by the new discretization scheme is found to be more accurate than that of other methods, especially at the vicinity of interfaces. This property is important when applying FFT methods to solve more complex problems like large strain deformation [START_REF] Lahellec | Analysis of inhomogeneous materials at large strains using fast fourier transforms[END_REF]. The new method also provides better estimates for the effective elastic moduli. Furthermore, its estimates does not depend on the reference medium, because the scheme is based on a finite-differences discretization of continuum mechanics. Although not explored in this work, the modified Green operator can be used with most other FFT iterative solvers, like the "augmented Lagrangian" [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF] or with FFT algorithms that are less sensitive to the reference [START_REF] Zeman | Accelerating a FFT-based solver for numerical homogenization of a periodic media by conjugate gradients[END_REF][START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF], leading to the same local fields.

Figure 2 :

 2 Figure 2: Elementary periodic domain Ω = [1/2; 1/2] 2 containing a square inclusion with elastic moduli µ 2 , κ 2 (top-left, shown in white) embedded in a matrix (shown in gray) with elastic moduli µ 1 , κ 1 .

Figure 3 :

 3 Figure 3: Stress component σ 12 (x) predicted by the various FFT schemes at the three resolutions L = 512, 1024 and 2048 (left to right) in the region [-0.04; 0.04] 2 . The center of the region is the bottom-left corner of the square inclusion in Fig. (2).

Figure 4 :

 4 Figure 4: 2D section of the stress component σ 12 (x) along the plane x 3 = -0.2461 predicted by the various FFT schemes at the three resolutions L = 256, L = 512 and L = 1024 (left to right). The section is parallel to one of the faces of the inclusion and close to the interface with the matrix.

Figure 5 :

 5 Figure 5: Apparent elastic modulus C 11,11 estimated by FFT methods using the Green operators G and G C,W,R (black and red), at increasing resolution L. Orange: estimate in [28]; violet: estimate of the asymptotic modulus C 11,11 using FFT data.

5 Figure 6 :

 56 Figure 6: Number of iterations N(E 0 ) required to achieve convergence, as a function of the reference Young modulus E 0 , for the accelerated schemes AS and AS C,W,R using various Green operators. Convergence is achieved when the precision η = 10 -8 is reached. The microsctructure is a Boolean model of quasi-porous spheres with χ = 10 -5 .

Figure 7 :

 7 Figure 7: Optimal reference Young modulus E 0 as a function of the contrast of properties χ, for the various FFT methods, in log-log scale. Direct schemes: black and red (nearly superimposed to one another); accelerated schemes: blue and green. Results for a porous material (χ = 0) are indicated at the left of the graph. The material is a Boolean model of spheres with volume fraction 17%.

Figure 8 :

 8 Figure 8: Number of iterations as a function of the contrast of properties χ, for the various FFT methods, in log-log scale. Direct schemes: black and red; accelerated schemes: blue and green. Results for a porous material (χ = 0) are indicated at the left of the graph. The material is a Boolean model of spheres with volume fraction 17%.

Figure 9 :

 9 Figure

4 Figure 10 :

 410 Figure

Acknowledgements

The author thanks Carnot M.I.N.E.S for support through grant 20531.

[1] S.-B. Lee, R. A. Lebensohn, and A. D. Rollett. Modeling the viscoplastic micromechanical response of two-phase materials using fast fourier transforms. International Journal of Plasticity, 27(5):707-727, 2011.