Sényo Komlavi

Eloh

Alain Jacques

Stéphane Berbenni
email: stephane.berbenni@univ-lorraine.fr

Development of a new consistent discrete Green operator for FFT-based methods to solve heterogeneous problems with eigenstrains

Keywords: FFT-based method, Discrete Green Operators, eigenstrains, Eshelby Inclusion, inhomogeneties, dislocations

published or not. The documents may come

Development of a new consistent discrete Green operator for FFT-based methods to solve heterogeneous problems with eigenstrains

Komlavi Sényo Eloh, Alain Jacques, Stephane

Introduction

During the last decades, the Fast Fourier Transform (FFT)-based method for homogenization and computation of mechanical fields in periodic media received a surge of interest since the pioneering work of Moulinec andSuquet (1994, 1998). In contrast with the standard finite element method, which needs a more or less complex meshing procedure depending on studied heterogeneous microstructures [START_REF] Quey | Large-scale 3D random polycrystals for the finite element method: generation, meshing and remeshing[END_REF], the FFT-based method has the advantage to evaluate the local and overall responses directly on regular grids from novel destructive or non-destructive characterization techniques [START_REF] Lebensohn | Orientation image-based micromechanical modeling of subgrain texture evolution in polycrystalline copper[END_REF][START_REF] Willot | Microstructure-induced hotspots in the thermal and elastic responses of granular media[END_REF][START_REF] Pokharel | Instantiation of crystal plasticity simulations for micromechanical modelling with direct input from microstructural data collected at light sources[END_REF]. A clear advantage also lies in the easy access to efficient FFT packages [START_REF] Frigo | An adaptive software architecture for the FFT[END_REF]. For all these advantages, this method was used to solve different problems in mechanics of materials. First of all, FFT-based methods and algorithms were developed to compute the mechanical fields and the effective properties of linear elastic composites (Moulinec andSuquet, 1994, 1998;[START_REF] Müller | Mathematical vs. experimental stress analysis of inhomogeneities in solids[END_REF]. Heterogeneous materials with eigenstrains or thermal strains were treated with different FFT-based methods [START_REF] Dreyer | An approximate analytical 2D-solution for the stresses and strains in eigenstrained cubic materials[END_REF][START_REF] Vinogradov | An accelerated FFT algorithm for thermoelastic and non-linear composites[END_REF][START_REF] Liu | Analytical solution for elastic fields caused by eigenstrains in a half-space and numerical implementation based on FFT[END_REF][START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF][START_REF] Bertin | A FFT-based formulation for efficient mechanical fields computation in isotropic and anisotropic periodic discrete dislocation dynamics[END_REF][START_REF] Donegan | Simulation of residual stress and elastic energy density in thermal barrier coatings using fast Fourier transforms[END_REF][START_REF] Graham | Fast fourier transform discrete dislocation dynamics[END_REF][START_REF] Jacques | From modelling of plasticity in single crystal superalloys to high resolution X-rays three-crystals diffractometer peaks simulation[END_REF][START_REF] Wang | Analytical solutions for elastic fields caused by eigenstrains in two joined and perfectly bonded half-spaces and related problems[END_REF]. The method was also applied to non-linear composite materials [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF], rigid-viscoplastic/ elasto-viscoplastic polycrystals [START_REF] Lebensohn | N-site modeling of a 3d viscoplastic polycrystal using Fast Fourier Transform[END_REF][START_REF] Lee | Modeling the viscoplastic micromechanical response of two-phase materials using Fast Fourier Transforms[END_REF][START_REF] Lebensohn | Fast fourier transform-based modeling for the determination of micromechanical fields in polycrystals[END_REF][START_REF] Lebensohn | An elasto-viscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF][START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF], continuum dislocation/ disclination mechanics [START_REF] Berbenni | A numerical spectral approach for solving elasto-static Field Dislocation and G-Disclination Mechanics[END_REF][START_REF] Brenner | Numerical implementation of static Field Dislocation Mechanics theory for periodic media[END_REF][START_REF] Berbenni | A fast Fourier transform-based approach for generalized disclination mechanics within a couple stress theory[END_REF][START_REF] Djaka | Field Dislocation Mechanics for heterogeneous elastic materials: A numerical spectral approach[END_REF], deformation twinning [START_REF] Arulkumar | Numerical study of the stress state of a deformation twin in magnesium[END_REF][START_REF] Mareau | Micromechanical modelling of twinning in polycrystalline materials: Application to magnesium[END_REF], fracture [START_REF] Herrmann | Linear and elastic-plastic fracture mechanics revisited by use of Fourier transforms: theory and application[END_REF], damage [START_REF] Spahn | A multiscale approach for modeling progressive damage of composite materials using Fast Fourier Transforms[END_REF], nano-polycrystals [START_REF] Upadhyay | A higher order elasto-viscoplastic model using fast Fourier transforms: effects of lattice curvatures on mechanical response of nanocrystalline metals[END_REF], non local polycrystal plasticity [START_REF] Lebensohn | Numerical implementation of non-local polycrystal plasticity using fast Fourier transforms[END_REF] and finite deformation [START_REF] Lahellec | Analysis of inhomogeneous materials at large strains using fast Fourier transform[END_REF][START_REF] Eisenlohr | A spectral method solution to crystal elasto-viscoplasticity at finite strains[END_REF][START_REF] Shanthraj | Numerically robust spectral methods for crystal plasticity simulations of heterogeneous materials[END_REF][START_REF] Degeus | Finite strain FFT-based non-linear solvers made simple[END_REF][START_REF] Nagra | Efficient fast Fourier transform-based numerical implementation to simulate large strain behavior of polycrystalline materials[END_REF][START_REF] Vidyasagar | Deformation patterning in finite-strain crystal plasticity by spectral homogenization with application to magnesium[END_REF]. However, FFT-based solutions generally present two issues that several authors have attempted to overcome in the last two decades. The first one is the improvement of the numerical convergence of FFT-based algorithm for arbitrary mechanical contrasts and the second one -on which this paper is focused -is the numerical treatment of spurious oscillations when local fields are studied in the presence of materials discontinuities. Since the original iterative fixed-point algorithm known as the "basic scheme" (Moulinec andSuquet, 1994, 1998) was seen to be efficient for small mechanical contrasts only, several progresses were made through the "accelerated scheme" [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF], the "augmented Lagrangian scheme" [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] and the "polarization-based scheme" [START_REF] Monchiet | A polarization based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. Furthermore, [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] proposed a variational formulation of the Lippmann-Schwinger equation based on the Hashin-Shtrikman energy principle.

From this energy-based scheme, the same authors [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] introduced a Galerkin approximation of the Lippmann-Schwinger equation and studied the convergence for arbitrary contrasts. [START_REF] Zeman | Accelerating a FFT-based solver for numerical homogenisation of periodic media by conjugate gradients[END_REF] proposed an accelerated FFT-based scheme using the conjugate gradient method. It was shown that this algorithm accelerates the convergence rate for problems with high mechanical contrasts.

An alternative Galerkin approximation based on trigonometric collocation was also proposed by [START_REF] Vondřejc | An FFT-based Galerkin method for homogenization of periodic media[END_REF] which was well suited to be used with the conjugate gradient method. This gradient descent method was adapted to non linear behaviors using the Newton-Raphson algorithm in [START_REF] Gélébart | Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials[END_REF] and [START_REF] Kabel | Efficient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deformation[END_REF]. Recently, an efficient fast gradient descent method was proposed by [START_REF] Schneider | An FFT-based fast gradient method for elastic and inelastic unit cell homogenization problems[END_REF] for elastic and (non linear) inelastic heterogeneous materials which requires less memory than conjugate gradient methods.

When materials discontinuities are present, numerical oscillations may occur, which can be due to Gibbs phenomenon or "aliasing" inherent to spectral methods. This question was first addressed by Moulinec andSuquet (1994, 1998) through a "truncated" Green operator, which simply consisted in disregarding the high-frequency term in the classic Green operator and replacing it by the inverse of the elastic stiffness tensor of the homogeneous reference medium. In order to avoid the artifacts related to highest frequencies, a second strategy used "attenuation factors" based on the discrete Fourier transform (DFT) of linear spline functions [START_REF] Kaßbohm | Improved approximations of Fourier coefficients for computing periodic structures with arbitrary stiffness distribution[END_REF] or low-pass filtering to dampen high frequency terms in the solutions [START_REF] Shanthraj | Numerically robust spectral methods for crystal plasticity simulations of heterogeneous materials[END_REF]. A third strategy first reported by [START_REF] Müller | Mathematical vs. experimental stress analysis of inhomogeneities in solids[END_REF], [START_REF] Dreyer | An approximate analytical 2D-solution for the stresses and strains in eigenstrained cubic materials[END_REF] and [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwise disordered porous media[END_REF] was to compute discrete Green operators based on Finite Difference discretization. This research direction was followed by many authors in the last few years [START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF][START_REF] Berbenni | A numerical spectral approach for solving elasto-static Field Dislocation and G-Disclination Mechanics[END_REF][START_REF] Willot | Fourier-based schemes for computing the mechanical response of composite with accurate local field[END_REF][START_REF] Berbenni | A fast Fourier transform-based approach for generalized disclination mechanics within a couple stress theory[END_REF][START_REF] Schneider | Computational homogenization of elasticity on a staggered grid[END_REF][START_REF] Lebensohn | Numerical implementation of non-local polycrystal plasticity using fast Fourier transforms[END_REF][START_REF] Schneider | Fft-based homogenization for microstructures discretized by linear hexahedral elements[END_REF][START_REF] Djaka | Field Dislocation Mechanics for heterogeneous elastic materials: A numerical spectral approach[END_REF][START_REF] Vidyasagar | Predicting the effective response of bulk polycrystalline ferroelectric ceramics via improved spectral phase field methods[END_REF][START_REF] Vidyasagar | Deformation patterning in finite-strain crystal plasticity by spectral homogenization with application to magnesium[END_REF]. While [START_REF] Müller | Mathematical vs. experimental stress analysis of inhomogeneities in solids[END_REF] and [START_REF] Dreyer | An approximate analytical 2D-solution for the stresses and strains in eigenstrained cubic materials[END_REF] considered a centered-finite difference-based scheme, [START_REF] Willot | Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwise disordered porous media[END_REF] reported forward and backward finite difference-based schemes. Among the different efficient recent schemes, a discrete Green operator based on a centered finite difference-based scheme on a rotated grid (i.e. a "rotated scheme") was recently proposed by [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composite with accurate local field[END_REF]. It was successfully applied to Field Dislocation Mechanics equations for heterogeneous materials in [START_REF] Djaka | Field Dislocation Mechanics for heterogeneous elastic materials: A numerical spectral approach[END_REF]. Higher-order centered difference-based schemes were also reported to be very efficient for different applications [START_REF] Vidyasagar | Predicting the effective response of bulk polycrystalline ferroelectric ceramics via improved spectral phase field methods[END_REF][START_REF] Vidyasagar | Deformation patterning in finite-strain crystal plasticity by spectral homogenization with application to magnesium[END_REF]. Another recent finite difference discretization scheme performed on a staggered grid was developed by [START_REF] Schneider | Computational homogenization of elasticity on a staggered grid[END_REF]. This refined scheme combined with DFT was seen to be efficient for three dimensional porous materials to give numerical solutions devoid of oscillations. Furthermore, a new discretization method based on linear hexahedral finite elements adapted for the basic scheme and the conjugate gradient method [START_REF] Schneider | Fft-based homogenization for microstructures discretized by linear hexahedral elements[END_REF] was reported and an equivalence has been demonstrated between the "rotated scheme" and hexahedral finite elements with reduced integration. In their Galerkin discretization of the Lippmann-Schwinger equation, [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] proposed a fourth strategy based on a "consistent discrete Green operator" which was adapted to a general variational framework based on the Hashin-Shtrikman energy principle. Later, it was however considered by [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] as too costly for three-dimensional computations and a filtering approach was finally used by these authors.

The objective of this paper is to propose new consistent periodized discrete Green operators named Discrete Green Operators (DGO), which are derived from the classical continuous Green tensor together with the Discrete Fourier Transform (DFT) in order to take explicitly into account the discreteness of the DFT. These new discrete operators are similar in spirit as the discrete Green operator developed by Brisard andDormieux (2010, 2012) but the expressions are different. These theoretical differences will be highlighted in the present paper. However, numerical comparisons with their operator will not be reported because they are out of the scope of the paper. It will be shown that the numerical oscillations near materials discontinuities originate from the commonly used Continuous Green Operator (CGO). This work is also focused on the development of a new DGO to compute the displacement field, which was recently already used to simulate X-rays diffraction peaks in other contributions [START_REF] Jacques | From modelling of plasticity in single crystal superalloys to high resolution X-rays three-crystals diffractometer peaks simulation[END_REF][START_REF] Eloh | The effect of crystal defects on 3D high-resolution diffraction peaks: A FFT-based method[END_REF]. Our starting point is the numerical resolution of the Lippmann-Schwinger equation using a DGO adapted to the numerical iterative scheme for heterogeneous materials with eigenstrains developed by [START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF].

The paper is organized as follows. In Section 2, the Lippmann-Schwinger equation for heterogeneous elastic materials with eigenstrains is recalled and its resolution in the Fourier space is reported. The iterative numerical scheme and the convergence criterion are described in section 2 as well. In Section 3, the DGO for displacements and strains are developed starting from a one-dimensional representation. The properties of these DGO are highlighted and are discussed in comparison with classic CGO.

In Section 4, several numerical results will be reported, such as the computation of the local stresses and displacements of composite materials with homogeneous or heterogeneous elasticity combined with a dilatational eigenstrain or an eigenstrain representing prismatic dislocation loops. The numerical stress and displacement solutions obtained with the DGO are discussed for cubic-shaped inclusions, spherical Eshelby and inhomogeneity problems. The results given by the DGO will be system-atically compared to analytical solutions and classic discretization method obtained with the CGO. In the whole manuscript, vectorial and tensorial fields will be written using bold characters. In some practical circumstances, tensors or vectors will be written with index notations using Cartesian coordinates and the Einstein's summation convention is used.

2 The Lippmann-Schwinger equation and its numerical spectral resolution

The Lippmann-Schwinger equation for eigenstrained heterogeneous materials

In this section, we present the field equations for a linear elastic problem with eigenstrains used to compute the response of a periodic composite medium. The periodic unit cell with volume Ω is subjected to a uniform overall strain tensor denoted E.

This overall strain corresponds to the spatial average of the strain field on the unit cell. Let us denote u, ε, σ, ε * the displacement, strain, Cauchy stress and eigenstrain fields, respectively.

Let us define a homogeneous linear elastic reference material with elastic stiffness C o such as at every position vector x ∈ Ω, the elastic stiffness tensor reads:

C(x) = C o + δC(x) (1)
where C is the fourth order tensor of spatially heterogeneous elastic moduli with classic major symmetries for components:

C ijkl = C jikl = C ijlk = C klij .
The tensor δC(x) denotes the spatial fluctuations of C(x) with respect to C o . For this problem, the balance of linear momentum without body force yields:

divσ(x) = 0 ∀x ∈ Ω (2)
where div is the divergence operator. Furthermore, from strain compatibility, the local total strain tensor ε(x) reads:

ε(x) = ∇ s u(x) (3)
where the symbol ∇ s denotes the symmetric part of the gradient operator ∇. In the present small strain setting, ε(x) additively decomposes as a sum of elastic strain ε e and eigenstrain ε * tensors as follows:

ε(x) = ε e (x) + ε * (x) = S(x) : σ(x) + ε * (x) (4)
where S = C -1 is the elastic compliance tensor. The eigenstrain field ε * is also called "stress-free strain" in the sense of [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion and related problems[END_REF]. It can represent thermal strains [START_REF] Vinogradov | An accelerated FFT algorithm for thermoelastic and non-linear composites[END_REF][START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF][START_REF] Donegan | Simulation of residual stress and elastic energy density in thermal barrier coatings using fast Fourier transforms[END_REF] or misfit strains [START_REF] Dingreville | From coherent to incoherent mismatched interfaces: A generalized continuum formulation of surface stresses[END_REF][START_REF] Jacques | From modelling of plasticity in single crystal superalloys to high resolution X-rays three-crystals diffractometer peaks simulation[END_REF]. The combination of Eq. 2, Eq. 3, Eq. 4 together with Eq. 1 and the introduction the Green's function method leads to an integral equation usually called "Lippmann-Schwinger equation" in the mechanics community, which reads:

ε(x) = E -(Γ * τ ) (x) ∀x ∈ Ω (5)
In Eq. 5, the symbol * denotes the spatial convolution product between two functions (in this case Γ and τ ). In this equation, Γ is the modified Green operator associated with C o and τ is the stress polarization field defined as:

τ (x) = δC(x) : ε(x) -C(x) : ε * (x) (6)
The numerical resolution of Eq. 5 is not an easy task for two main reasons. First, it contains a spatial convolution of two functions in real space, and the numerical full field computation of a convolution in the real space is complicated. Second, it is an implicit equation, which needs specific resolutions methods. The formal strain solution writes as a power series expansion [START_REF] Vinogradov | An accelerated FFT algorithm for thermoelastic and non-linear composites[END_REF]:

ε(x) = +∞ n=0 [-(Γ * δC)(x)] n : (E + (Γ * C : ε * ) (x)) (7)
To overcome these difficulties, Eq. 5 can be written in the Fourier space so that the spatial convolution is transformed to a simple product. For periodic media, the numerical Fast Fourier Transform (FFT)-based resolution initially proposed by Moulinec andSuquet (1994, 1998) is an efficient way to iteratively solve the implicit equation in order to obtain the strain solution.

Fourier space resolution

From the convolution theorem in the Fourier space [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF], the Lippmann-Schwinger equation (Eq. 5) can be transformed in the Fourier space as follows:

ε(ξ) = -Γ(ξ) : τ (ξ) ∀ξ = 0, ε(0) = E (8)
where Γ and τ are the Fourier transforms of Γ and τ , respectively. ξ represents the Fourier vector. In Eq. 8, Γ writes in index notations:

Γ ijql (ξ) = 1 4 ( G iq ξ j ξ l + G il ξ j ξ q + G jq ξ i ξ l + G jl ξ i ξ q ) (9)
where G is the Fourier transform of the elastic Green tensor defined as [START_REF] Mura | Micromechanics of defects in solids[END_REF]:

               G ij (ξ) = N ij (ξ) D(ξ) ∀ξ = 0 G ij (0) = 0 (10)
where N ij (ξ) denotes the rectangular components of the cofactor matrix related to the acoustic tensor

F ij = C o ijql ξ q ξ l and D(ξ) is the determinant of F ij . Due to the symmetry properties of C o ijql , N ij (ξ) satisfies: N ij (ξ) = N ji (ξ), therefore G ij = G ji
and Γ ijql = Γ jiql = Γ ijlq = Γ qlij . For isotropic elasticity using the two Lamé constants, the expression of Γ ijql was given for example in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF].

Using Eq. 3 and Eq. 8, the Fourier transform of the displacement field u reads:

u(ξ) = -B(ξ) : τ (ξ) (11)
where the expression of the third order tensor B writes in indicial form:

B ijq (ξ) = i 2 ( G ij ξ q + G iq ξ j ) ( 12
)
where i denotes the imaginary unit. The third order tensor B has the following symmetries: B ijq = B iqj due to the symmetry of the stress polarization tensor (τ ij = τ ji ).

Recently, [START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF] proposed an iterative scheme for eigenstrain problems of the augmented Lagrangian (AL) type [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] to solve Eq. 5 using

Eq. 8. The FFT-based numerical algorithm used for heterogeneous elastic materials with eigenstrains is now briefly presented.

Numerical Algorithm

The Discrete Fourier Transform (DFT) method is used. The unit cell with volume Ω is now discretized into a regular grid of N 1 × N 2 pixels in 2D problems or N 1 × N 2 × N 3 voxels in 3D problems. For 3D applications, the voxel coordinates and the corresponding discrete frequencies in Fourier space are respectively [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]:

x d (i 1 , i 2 , i 3 ) = (i 1 -1) T 1 N 1 , (i 2 -1) T 2 N 2 , (i 3 -1) T 3 N 3 , with i 1 = 1, ....N 1 i 2 = 1, ....N 2 i 3 = 1, ....N 3
and T i denotes the period of the unit cell in the i th direction. The discrete frequencies are given by:

ξ i = (- N i 2 + 1) 1 T i , (- N i 2 + 2) 1 T i , ..., - 1 T i , 0, 1 T i , ..., ( N i 2 -1) 1 T i , ( N i 2 ) 1 T i
if N i is even and

ξ i = - N i -1 2 1 T i , ..., - 1 T i , 0, 1 T i , ..., N i -1 2 1 T i if N i is odd.
The iterative fixed-point algorithm for heterogeneous elastic materials with eigenstrains and for macroscopic stress-free states (i.e. no applied overall stress) is now described following the "thermoelastic FFT"-based method developed by [START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF]:

Initialization : (a 0 ) E (0) = ε * (x d ) (b 0 ) ε o (x d ) = E (0) ∀x d ∈ Ω (c 0 ) Compute Γ(ξ j ) ∀ξ j Iteration(i + 1) ε i (x d ) and E i known (a) τ i+1 (x d ) = δC(x d ) : ε i (x d ) -C(x d ) : ε * (x d ) (b) τ i+1 (ξ j ) = F F T (τ i+1 (x d )) (c) ε i+1 (x d ) = E i -F F T -1 ( Γ(ξ j ) : τ i+1 (ξ j ))(x d )) (d) σ i+1 (x d ) = C(x d ) : (ε i+1 (x d ) -ε * (x d )) (e) E i+1 = ε i+1 (x d ) -C o -1 (x d ) : σ i+1 (x d ) (f ) Convergence test : ||ε i+1 (x d )-ε i (x d )|| ||E i+1 || ≤ e (13)
In the algorithm (Eq. 13), F F T and F F T -1 denote the Fast Fourier Transform and the Inverse Fast Fourier Transform respectively. The symbol . indicates a volume average over the unit cell and the symbol || . || denotes the Euclidian norm of a second-order tensor. At step (e), the macroscopic strain field is iteratively adjusted [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF][START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF]. In the inhomogeneous elastic case, the convergence criterion to stop the iterative procedure is defined on the strain field (see step (f)) when a given error e is reached. To have enough accuracy, this value should range from 10 -6 to 10 -10 depending on the problem. In the case where the tensor of the elastic stiffness tensor is not a function of the position (

C ijql (x d ) = C o ijql ,
i.e. homogeneous elasticity case), the algorithm converges in one iteration plus one more iteration to correct the applied strain field because the stress polarization is now explicit and equal to τ (x d ) = -C o : ε * (x d ) (using Eq. 6) with δC(x d ) = 0 at every point of the unit cell.

Difficulties linked to the Discrete Fourier Transform (DFT) method

Some difficulties appear in the use of the Discrete Fourier Transform (DFT) due to sampling interval [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF][START_REF] Press | Numerical recipes in C++[END_REF]. In the FFT-based algorithm (Eq. 13), the applicability of the DFT is based on the assumption that the relevant functions are periodic. If it is not the case and these functions are not bandwidth limited to less than the Nyquist critical frequency, the so-called "aliasing" effect [START_REF] Press | Numerical recipes in C++[END_REF] (pages 505-508) tends to make the solution inaccurate near the end points of the computational domain. Here, the continuous Green operators (CGO) are not periodic and this non periodicity becomes an issue for space limited computational domain. In the classic algorithm initially proposed by Moulinec andSuquet (1994, 1998), the modified Green tensor Γ involved in the Lippmann-Schwinger equation was not periodized. To fix the problem, [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] proposed a method which consists in using another expression of the Green operator Γ at the highest frequencies ξ j = ±( N j 2 -1) 1 T j (where T j denotes the period of the unit cell):

Γ = C o-1
At theses frequencies, the stress is forced to be zero when the convergence is reached.

For these authors, this choice allowed to avoid most of oscillations on the mechanical fields. More recently, [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] developed a periodized discrete Green operator adapted for their algorithm based on the Hashin-Shtrikman varia-tional approach. To our knowledge, a consistent periodized discrete Green operator (DGO) for the Lippmann-Schwinger equation has not been developed yet. The objective of the next section is to give some mathematical details on the development of such consistent DGO that will be applied to the thermoelastic FFT-based algorithm introduced in Eq. 13.

3 Development of consistent periodized discrete Green operators (DGO)

One-dimensional case

First, let us consider a one-dimensional problem for the sake of clarity. We denote e(x) a periodic function of a real variable x that can be identified as the strain field or the displacement field. This function is defined as a spatial convolution of two functions: t(x) a periodic function, which is similar to the stress polarization field and g(x) which is similar to a Green operator. By denoting d the spatial period of e(x) and t(x), T (ξ) the Fourier transform of t(x), we can write (here the notation j = √ -1 is adopted for the imaginary unit):

t(x) = ∞ k=-∞ T (ξ)exp(j2πξ h x) (14)
with ξ h = h d and:

T (ξ) = 1 d +d/2 -d/2 t(x)exp(-j2πξ h x)dx (15)
The function t(x) is periodic, then:

T (ξ) = 1 d d 0 t(x)exp(-j2πξ h x)dx
For N points along the x-axis, let us now suppose the following spatial discretization for their x position: 1). Let us assume that t(x) is approximated with the discrete stepwise uniform function The Fourier transform of t d (x i ) is denoted T d (ξ h ). Using the discrete piece-wise approximation of t(x) by t d (x i ) and the fact that:

x i = i d N + d 2N where i ∈ [0, N -1] (see Fig.
d 0 ≈ N -1 i=0 x i + d 2N x i -d 2N the Fourier transform of t d (x i ) reads: T d (ξ h ) = 1 d N -1 i=0 t d (x i ) x i +d/2N x i -d/2N exp(-j2πξ h x )dx (16)
Using the change of variables X = x -x i yields:

T d (ξ h ) = F T [t d (x i )](ξ) = 1 d N -1 i=0 t d (x i )exp(-j2πξ h x i ) d/2N -d/2N exp(-j2πξ h X)dX (17)
From the formula:

+d/2N -d/2N exp(-j2πξ h X)dX = d N sinc πξ h d N
we obtain:

T d (ξ h ) = 1 N sinc πξ h d N N -1 i=0 t d (x i )exp(-j2πξ h x i ) (18)
In Eq. 18, the Discrete Fourier Transform (DFT) of t d (x i ) can be identified:

T D (ξ h ) = N -1 i=0 t d (x i )exp(-j2πξ h x i ) (19)
Therefore, Eq. 18 becomes:

T d (ξ h ) = 1 N T D (ξ h )sinc π k N ( 20
)
where k varies from -∞ to +∞. Using a change of variables: k = mN + l, where l ∈ [0, N -1] and m varies from -∞ to +∞, the last equation writes:

T d (ξ h ) = (-1) m π(mN + l) T D (ξ h )sin π l N (21)
Eq. 21 defines the Fourier transform of a continuous periodic function as a function of the Fourier transform of the discrete stepwise function t d (x i ).

With the previous definitions, the Fourier transform of the function g(x) (Green operator) is denoted G(ξ). It is defined through the relationship between T (ξ) the Fourier transform of t(x) and E(ξ) the Fourier transform of e(x):

E(ξ) = G(ξ) T (ξ) (22)
with:

e(x) = ∞ k=-∞ E(ξ h )exp(j2πξ h x) (23)
where ξ h = k d denotes the discrete frequencies. Using the variable k = mN + l as introduced before, e(x) writes:

e(x) = N -1 l=0 ∞ m=-∞ E(ξ mN +l )exp(j2πξ mN x)exp(j2πξ l x) (24)
Here, e(x) is approximated by e d (x i ):

e d (x i ) = N -1 l=0 ∞ m=-∞ E(ξ mN +l )exp(j2πξ mN x i )exp(j2πξ l x i ) (25)
Using

x i = i d N + d 2N (
where i=0,1,.....N-1) leads to:

e d (x i ) = N -1 l=0 exp(j2πξ l x i ) ∞ m=-∞ (-1) m E(ξ mN +l ) (26)
Here, we can remark that if e d (x i ) has a finite value at all points x j in the periodic domain, then the infinite sum

∞ m=-∞ (-1) m E(ξ mN +l ) is also finite at all Fourier points ξ l with l ∈ [0, N -1]
. Now, we can use the definition of the inverse DFT to calculate e d (x i ):

e d (x i ) = 1 N N -1 l=0 E D (ξ l )exp(j2πξ l x i ) (27)
By comparing Eq. 26 and Eq. 27, the DFT of e(x) can be identified as:

E D (ξ l ) = N ∞ m=-∞ (-1) m E(ξ mN +l ) (28)
where E D (ξ l ) is finite at all Fourier points ξ l . For each discrete k = mN + l, we directly obtain from Eq. 22:

E(ξ mN +l ) = G(ξ mN +l ) T (ξ mN +l ) (29)
The expression of T (ξ mN +l ) is obtained from Eq. 21:

T d (ξ mN +l ) = (-1) m π(mN + l) T D (ξ mN +l )sin π l N (30)
Then, we have from the DFT:

T D (ξ mN +l ) = N -1 l=0 t d (x i )exp(-j2πξ mN +l x i ) = (-1) m T D (ξ l ) (31)
Finally, Eq. 30 reads:

T d (ξ mN +l ) = 1 π(mN + l) T D (ξ l )sin π l N (32)
By using Eq. 28, Eq. 29 and Eq. 32:

E D (ξ l ) = N ∞ m=-∞ G(ξ mN +l )(-1) m sin πl N π(mN + l) T D (ξ l ) (33)
From this last equation, we deduce the final expression of the periodized discrete Fourier transform of g(x) denoted G D (ξ l ):

G D (ξ l ) = N π sin πl N ∞ m=-∞ (-1) m mN + l G(ξ mN +l ) (34)
G D (ξ l ) represents a discrete Green operator (DGO) for this one-dimensional problem which converges to a finite value at all Fourier points

ξ l with l ∈ [0, N -1] due to the definition E D (ξ l ) = G D (ξ l ) T D (ξ l
) and the convergence of E D (ξ l ) to a finite value (see above). In the next sub-section, we will define two consistent periodized DGO for the three-dimensional problem to solve the Lippmann-Schwinger equation (i.e.

the strain solution) and the displacement field solution.

3.2 Consistent periodized discrete Green operators Γ D (ξ ijq ) and B D (ξ ijq )

The theoretical result given by Eq. 34 for the one-dimensional discrete configuration can be easily extended to obtain consistent three-dimensional periodized discrete Green operators (DGO) related to tensors Γ(ξ) and B(ξ) defined in section 2.2 through Eq. 9 and Eq. 12 respectively. With previous notations, the single index l in Eq. 34 can be replaced by three indices i, j, k to obtain two DGO expressions with triple infinite sums from the definitions:

ε D (ξ ijq ) = -Γ D (ξ ijq ) : τ D (ξ ijq ) ∀ξ ijq = 0 (35)
and:

u D (ξ ijq ) = -B D (ξ ijq ) : τ D (ξ ijq ) ∀ξ ijq = 0 (36)
such that:

Γ D (ξ ijq ) = A ijq ∞ m,n,p=-∞ (-1) m+n+p mN + i 1 nN + j 1 pN + q Γ(ξ mN +i,nN +j,pN +q ) (37)
and:

B D (ξ ijq ) = A ijq ∞ m,n,p=-∞ (-1) m+n+p mN + i 1 nN + j 1 pN + q B(ξ mN +i,nN +j,pN +q ) (38)
where:

A ijq = N π 3 sin πi N sin πj N sin πq N (39)
and ξ ijq are three-dimensional discrete frequencies with i ∈ [0,

N -1], j ∈ [0, N - 1], q ∈ [0, N -1].
Using the same arguments as for the one-dimensional case, it is straightforward to show that both infinite series Γ D (ξ ijq ) and B D (ξ ijq ) converge to finite values at all three-dimensional Fourier points ξ ijq . Numerical experiments (see section 4) will

show that these series converge slowly. However, it is observed that B D has a better convergence rate than Γ D . This can be explained by using Eq. 10 and the fact that N ij is a polynomial homogeneous function of degree 4 and D is a polynomial homogeneous function of degree 6, so according to Eq. 10 together with Eq. 9 or

Eq. 12, the Γ ijql components are homogeneous functions of degree 0 and the B ijq components are homogeneous functions of degree -1. The numerical computations of Γ D and B D are performed only once before performing the iterative algorithm and taking advantage of their symmetry properties:

Γ D ijql = Γ D jiql = Γ D ijlq = Γ D qlij ; Γ D ijql (-ξ) = Γ D ijql (ξ) and B D ijq = B D iqj ; B D ijq (-ξ) = B D ijq (ξ).
In section 4, a parametric study and examples will be proposed on three-dimensional computational grids to evaluate the convergence of the developed discrete consistent Green operators Γ D (ξ ijq ) and B D (ξ ijq ) that have the forms of infinite series by choosing a satisfactory truncation order to obtain accurate numerical results in comparisons with analytical solutions. Truncated discrete Green operators at the k-order will be denoted DGO k and the infinite sums on m, n, p will be replaced by finite sums on m, n, p in the interval [-k, k -1] where k is an integer. As shown in Appendix A, the values of the stepwise strain/displacement fields computed with DGO k at the center of the N 3 voxels are equal to those computed at the same points with the CGO (i.e. Γ) on a subdivided (2kN ) 3 grid. As the stress and strain fields have discontinuities at the boundaries of each of the N 3 voxels, a full computation with Γ on the (2kN ) 3 computational grid exhibits the usual Gibbs oscillations at the voxel edges. However, as the mechanical fields are actually computed at the centers of the N 3 voxels, i.e. far from the discontinuities, no oscillation occurs when using the DGO k .

Comparison with the consistent periodized discrete Green operator introduced

by Brisard andDormieux (2010, 2012) The discrete Green operator proposed in Eq. 37 is very similar in spirit to the consistent discrete Green operator developed by Brisard andDormieux (2010, 2012).

However, their operator has some major difference with Eq. 37. To compare both discrete Green operators, let us consider a two-dimensional computational grid as in Brisard andDormieux (2010, 2012). Therefore, Eq. 37 reduces to:

Γ D (ξ ij ) = N π 2 sin πi N sin πj N ∞ m,n=-∞ (-1) m+n mN + i 1 nN + j Γ(ξ mN +i,nN +j ) (40)
Using the definition of the sinc function, we have:

sinc π(mN + i) N = (-1) m N π sin πi N 1 mN + i = (-1) m i (mN + i) -1 sinc πi N (41)
Then, Eq. 40 can be rewritten as follows:

Γ D (ξ ij ) = (ij)sinc πi N sinc πj N ∞ m,n=-∞ (-1) m+n mN + i 1 nN + j Γ(ξ mN +i,nN +j ) (42)
which is different from the expression given by Brisard and Dormieux (2010) (p.

667):

Γ D (ξ ij ) = (i 2 j 2 )sinc 2 πi N sinc 2 πj N ∞ m,n=-∞ 1 (mN + i) 2 1 (nN + j) 2 Γ(ξ mN +i,nN +j ) (43)
By comparing Eq. 42 with Eq. 43, it is remarkable to see that both expressions are different in the sense that the discrete Green operator of Brisard andDormieux (2010, 2012) 7) and ( 12). Let us note that both discrete operators satisfy: Γ D (ξ 00 ) = 0. In this paper, the discrete Green operator derived by Brisard andDormieux (2010, 2012) will not be discussed anymore in the numerical examples since it is designed for a Galerkin based algorithm, i.e. starting from the weak form of the Lippmann-Schwinger equation instead of its strong form as in section 2.4.

Results and discussion

In the forthcoming applications, the unit cell with volume Ω is discretized in 128 ×

128 × 128 voxels. The heterogeneous material is only subjected to an eigenstrain field (which is a stress-free strain). As only eigenstrains are imposed, the overall applied stress is set to zero. The overall strain is unknown and updated at each step following the algorithm described in section 2. The detailed eigenstrain fields will be specified for each application. In this section, the FFT-based method using the classical continuous Green operator will be denoted CGO (Continuous Green Operator) [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] and the FFT-based method with the proposed discrete Green operator (computed with a given truncation order k) will be denoted DGO k (Discrete Green Operator). First, in section 4.1, a cubic-shaped inclusion with eigenstrain field is studied in homogeneous isotropic elasticity and the obtained stress fields are calculated using the DGO with different orders k as well as the CGO approximation. The numerical results will be compared to analytical solutions considered as periodic solutions started from the solution given by Li and Anderson (2001) (see Appendix B). In section 4.2, the case of a thin square-shaped inclusion with one voxel in the thickness is considered assuming homogeneous isotropic elasticity and an eigenstrain representing a prismatic dislocation loop. In this case, we will carefully examine the displacement profiles. The displacement vector is numerically computed with the CGO and DGO k approximations and is compared to analytical solutions (see Appendix B). In section 4.3, a spherical Eshelby-type problem is reconsidered to assess the present numerical method based on DGO k . In section 4.4, a spherical inhomogeneity solution with eigenstrain is considered and discussed using low contrasted elastic moduli between phases as in [START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF]. The presence of the eigenstrain field is at the origin of the development of inelastic In order to quantify the convergence of periodized DGO k , the relative difference is calculated between the analytical solution and the numerical solutions for σ 33 obtained at different locations from DGO k or CGO. This relative difference is defined as:

=| σ num 33 -σ anal 33 σ anal 33 |
The values of are computed at the specific points O, A, B, C of the inclusion. Fig. 5 shows the relative difference plotted as a function of k using logarithmic scales. It is observed that decreases as k increases. Besides, the value of calculated with the DGO k approximation reaches an admissible value (i.e. < 5×10 -3 ) for a truncation of the DGO corresponding to k = 81. have enough accuracy for this example. However, the computational time is more important with DGO k : it is multiplied by (2k) 3 (due to 3D infinite sum containing in Eq. 37) but in practice, this multiplicative factor is less important. For example, the computational time for computed CGO on our cluster (Intel Xeon X5550

2.66GHz CPU 46 nodes, 48Go per node) is 7 seconds and the computational time for DGO 81 is 116913 seconds. The numerical burden is thus very important but DGO 81 is computed only one time and can be stored for later use. In addition, it is not always necessary to set k = 81. In Fig. 5, it is shown that it is possible to have also acceptable results for lower k, except at the cubic-shaped center.

In addition to accuracy, the most important advantage of DGO k is its beneficial effect on the oscillations that appear on some components of the stress field (e.g.

Eshelby-like prismatic dislocation loop in a homogeneous isotropic elastic periodic medium

Let us consider the same homogeneous isotropic elastic constants and uniform eigenstrain tensor as in the previous case but here a square-shaped inclusion is discretized in 32 × 32 × 1 voxels corresponding to an Eshelby-like squared prismatic dislocation loop. Here, the studied mechanical field is the displacement field, which means that the third order tensor B D is exploited (Eq. 38). The aim of this example is to study the convergence of the associated DGO k when k increases, its accuracy on the displacement field and its capability to completely suppress oscillations on the displacement field at the materials discontinuity between the inclusion and the matrix along the line d 3 in the z-direction (see Fig. 2). The eigenstrain field (ε * ij = 0 for all i, j except ε * 33 = 1) is responsible for an equivalent Burgers vector for this prismatic loop equal to b = 2a 3 ε * 33 where 2a 3 is the thickness of the inclusion in the z-direction (i.e. the voxel size) as for the analytical solution (see Appendix B). DGO k or CGO. These relative differences are defined as: the DGO k is more accurate than the CGO from k = 5 where 1 < 10 -3 . Regarding 3 (Fig. 8(b)), the solution given by the DGO k is more accurate than the CGO for all k values with 3 < 10 -3 .

1 =| u num 1 -u anal 1 u anal 1 | 3 =| u num 3 -u anal
Furthermore, on almost all components, the displacement field giving by the FFTsolution using the CGO based on B shows numerical oscillations due to the physical displacement discontinuity embodied by the Burgers vector describing strong plastic incompatibility at the interface between the inclusion and the matrix. These oscillations are completely removed by replacing B with the proposed periodized discrete operator B D and the DGO 21 (k=21) solution is a good compromise between accuracy as seen in Fig. 7 and computational cost (k remains low enough).

To finish the discussion on this example, let us focus on the impact of the Fourier for two different configurations: a unit cell constituted of 128 × 128 × 128 voxels with a centered square inclusion of 32 × 32 × 1 voxels and a unit cell constituted of 32 × 32 × 32 voxels with a centered square inclusion of 8 × 8 × 1 voxels. Fig. 9 shows the normalized component u 3 /b of the displacement field along the line d 3 for these two configurations. We first notice that the CGO exhibits numerical oscillations close to the discontinuity of u 3 which should theoretically be equal to b. The magnitude of these oscillations decreases as the number of voxels increases. However, these oscillations are not completely removed with the CGO and are always present even for a good resolution of 128 × 128 × 128 voxels. The new periodized discrete operator B D here embodied by the DGO 21 allows to remove all oscillations for both resolutions. This result shows a non negligible advantage of the DGO over the CGO to give accurate results devoid of oscillations even for low Fourier grid resolutions to describe a given physical problem.

Spherical inclusion problem with eigenstrains in a homogeneous isotropic elastic periodic medium

In this section, we evaluate the accuracy of the proposed DGO solution by comparing the numerical results to the reference Eshelby solution for a spherical inclusion problem with eigenstrain in a homogeneous isotropic elastic periodic medium. According to the Eshelby's theoretical result, the stress is uniform inside an ellipsoidal inclusion [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion and related problems[END_REF] whereas it is non uniform outside an ellipsoidal inclusion [START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF][START_REF] Mura | Micromechanics of defects in solids[END_REF]. Let us note that for the Eshelby's inclusion solution, the ellipsoidal inclusion (here spherical) is embedded in an infinite matrix. In order to verify that an inclusion radius of 14 voxels for a unit cell of 128 × 128 × 128 voxels is a good approximation, an interesting study was reported by [START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF] who used classic CGO. This study showed that this choice is reasonable and leads to a good match between the numerical FFT-based solution for periodic medium and Eshelby analytical solutions for infinite medium. Here, the same inclusion size is chosen to test the DGO 81 as it was shown in section 4.1 to be a relevant choice for the cubic-shaped inclusion.

Here, we use the same isotropic elastic constants and eigenstrain tensor as [START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF] in order to make comparisons with their work. The homogeneous Young's modulus is 65.4GPa and the homogeneous Poisson ratio is 0.42. The eigenstrain field contained in the inclusion is defined as a pure dilatational eigenstrain:

ε * ij = 0.005 if i = j and ε * ij = 0 if i = j.
The stress field along a centerline parallel to the x-axis is represented for three

x i = i d N .
cases in Fig. 10: CGO, DGO 81 and Eshelby's solution. The constant value of the stress component σ 11 for the interior points predicted by Eshelby's analytical result is -375.9M P a. Both FFT-methods based on CGO and DGO 81 give approximately the same value: -374M P a. Thus, the corresponding error is about 0.5% which shows a good accuracy without oscillation. This figure also shows a special effect at the interface between inclusion and matrix where a large deviation from the Eshelby's solution is observed. This error is not due to oscillations as observed before for the cubic-shaped inclusion near discontinuities but is due to so-called "dangling voxels" which are isolated voxels located at the poles of the spherical inclusion. This discretization problem was explained by [START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF] who suppressed the problem by increasing the total number of voxels of the unit cell up to 1024 × 1024 × 1024 voxels. Here, the problem is solved by changing the spatial definition of the spherical inclusion with regard to the spatial discretization: a voxel should be assigned to the inclusion if its distance from the center is less than the given radius of the inclusion instead of less than or equal as frequently used corresponding to the spatial discretization:

x i = i d N .
Therefore, the problem is naturally solved with the real space discretization proposed in Fig. 1:

x i = i d N + d 2N .
As example, Fig. 11 shows the stress components σ 11 and σ 33 obtained with the present spatial discretization described in Fig. 1. It is shown that the "dangling voxels" effect is eliminated and the CGO and the DGO 81 solutions for σ 11 and σ 33 give similar and accurate results for these stress components.

Spherical inhomogeneity problem with eigenstrains

Here, the same spherical configuration as in section 4.3 is considered but with elastic inhomogeneity in addition to eigenstrain. An inhomogeneity radius of 14 voxels is still used for a unit cell of 128 × 128 × 128 voxels. The isotropic elastic stiffness tensor of the inhomogeneity is different from the one of the matrix phase. Like in [START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF], the bulk modulus of the matrix phase is fixed and the bulk given by the Eshelby's analytical solution and computed with CGO and DGO 81 . Here, the following spatial discretization is used:

x i = i d N + d 2N
modulus of the inhomogeneity phase is defined by the ratio R = K m /K i = 2, which means that the inhomogeneity i is more compliant than the matrix m. Here, the case R = 0.5, i.e. an inhomogeneity stiffer than the matrix, is also investigated to test the numerical convergence for both inverse contrasts, see Fig. 13. The isotropic Young's modulus of the matrix phase is E m = 65.4GPa and the Poisson ratio for matrix and inhomogeneity phases is ν m = ν i = 0.42. This corresponds to a small mechanical contrast that can be treated with the fixed-point algorithm described in Eq. 13.

The elastic constants of the reference medium are chosen as:

K 0 = β(K m + K i ),
ν 0 = ν m = ν i = 0.42 with β = 0.5. A special attention is given to the number of iterations needed for convergence and to the stress solution accuracy using the CGO or the DGO. At each iteration, the error for testing convergence is reported using the strain field criterion as defined in the step (f) of Eq. 13. The iterative algorithm stops when the error is less than or equal to e = 1 × 10 -6 to have enough accuracy on stress fields. The obtained mechanical fields are compared to available analytical solutions. Analytical solutions for spherical inhomogeneity problems in infinite medium were derived by [START_REF] Mura | Micromechanics of defects in solids[END_REF] and [START_REF] Ju | Effective elastoplastic behavior of metal matrix composites containing randomly located aligned spheroidal inhomogeneities. Part i: micromechanics-based formulation[END_REF]. along a centered line parallel to the x-axis for the spherical inhomogeneity problem with eigenstrain, given by analytical solutions [START_REF] Mura | Micromechanics of defects in solids[END_REF][START_REF] Ju | Effective elastoplastic behavior of metal matrix composites containing randomly located aligned spheroidal inhomogeneities. Part i: micromechanics-based formulation[END_REF]) and numerically computed with CGO and DGO 81 . Here, R = K m /K i = 2.

In Fig. 12, both FFT-based results given by the CGO and the DGO 81 are compared with the analytical solutions for normalized stress field components σ 11 /E m and σ 33 /E m along a centered line parallel to the x-axis. The analytical uniform value of the stress field in the inhomogeneity is -317M P a for both components. Both CGO and DGO 81 give approximately same numerical uniform values inside the inhomogeneity for these components, which is -313.5M P a. This corresponds to a relative error of 1% which is acceptable.

The number of iterations before reaching convergence is 15 iterations regardless of the used Green operator. Furthermore, Fig. 13 describes the error as a function of iteration number up to 50 iterations. For this mechanical contrast, it is observed that the errors calculated at each iteration with CGO and DGO 81 are very close to each other. However, beyond 20 iterations at same iteration number, the error obtained with DGO 81 is a little bit smaller than the one with CGO. This observation means that in the heterogeneous case, the accuracy with DGO 81 is better than with CGO. cation order k should be increased for higher contrasts.

In addition, 2D sections of the stress component σ 12 normalized by 2µ m computed with CGO and DGO 81 are reported on Fig. 15 to study stress patterns close to the interface between the matrix and the inhomogeneity. The section is a cut parallel to (y, z) plane intersecting the center of the inhomogeneity along the x-axis. Following also reported by [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composite with accurate local field[END_REF]. They are clearly due to numerical oscillations which spoil the stress solution as already observed in Fig. 14. In contrast, the application of DGO 81 allows to completely remove these artificial patterns, see Fig.

T D (ξ h ) = N k-1 m=-k (-1) m T d 2kN (ξ mN +h ) (A.7)
We also restrict the T D 2kN sum (Eq. A.6) to its first term, i.e. Γ truncated to -kN, ., kN -1. The inverse DFT of the truncated Γ, computed on 2kN points, will exhibit Gibbs oscillations at discontinuities, i.e. at points y l = x i + d 2N . However, these oscillations should be damped for y l = x i , where t d is actually computed, thus the inverse DFT of T D k is free of the Gibbs oscillations.

The whole proof of section 3.1 can now be performed again for T D 2kN instead of T D . As a conclusion, the discrete operator DGO k with truncation order k is thus equivalent to use the CGO (i.e. the classic Γ operator) on a 2k times finer resolution. The main advantage to use the DGO k is thus both decreases of memory storage and CPU time for a same accuracy as compared to the CGO.

B Derivations of analytical stress and displacement components for periodic distributions of cubic-shaped inclusions with eigenstrains in isotropic elasticity

Stress components

First of all, the analytical stress solution for the stress field due to an eigenstrain ε * ij in a cubic-shaped isotropic elastic inclusion embedded in an infinite isotropic elastic medium was reported by [START_REF] Chiu | Stress-field due to initial strains in a cuboid surrounded by an infinite elastic space[END_REF]; [START_REF] Li | A compact solution for the stress field from a cuboidal region with a uniform transformation strain[END_REF]. The cuboidal inclusion dimensions are 2a 1 , 2a 2 , 2a 3 in the x, y, z directions, respectively. We consider the following homogeneous elastic properties in the medium: shear modulus µ, Young's modulus E and Poisson ratio ν. In the case where the only non zero eigenstrain component is ε * 33 , the analytical solution reads (see eq. ( 5) in [START_REF] Li | A compact solution for the stress field from a cuboidal region with a uniform transformation strain[END_REF]):

C n2 C n3 C n1 Rn -C n1 C n2 C n3 [C 2 n2 +C 2 n3 ]Rn ] σ 33 (x) Eε * 33 = 1 8π(1-ν 2 ) 8 n=1 (-1) n [-2 arctan C n1 C n3 C n2 Rn + C n1 C n2 C n3 [C 2 n1 +C 2 n3 ]Rn -2 arctan C n2 C n3 C n1 Rn + C n1 C n2 C n3 [C 2 n2 +C 2 n3 ]Rn ] (B.1)
In this expression, the distance R n between the position vector x(x 1 , x 2 , x 3 ) and the nth corner of the cuboidal inclusion (n = 1, ..., 8) is defined as:

R n (x) = C 2 n1 + C 2 n2 + C 2 n3
where the components C nj are given by:

[C nj ] =                                             
x 1 -a 1 x 2 -a 2 x 3 -a 3

x 1 + a 1 x 2 -a 2 x 3 -a 3

x 1 + a 1 x 2 + a 2 x 3 -a 3

x 1 -a 1 x 2 + a 2 x 3 -a 3

x 1 -a 1 x 2 + a 2 x 3 + a 3

x 1 -a 1 x 2 -a 2 x 3 + a 3

x 1 + a 1 x 2 -a 2 x 3 + a 3

x 1 + a 1 x 2 + a 2 x 3 + a 3                                              j = 1, 2, 3
In order to derive the analytical stress solution to be compared with the FFT-based results, an infinite periodic distribution of cuboidal inclusions with the same periods as in the FFT-based calculations should be considered, which theoretically corresponds to an infinite sum of individual stress contributions using Eq. (B.1). In order to reduce the computations, we first calculate in the unit cell the total stress given by the superposition of 27 × 27 × 27 periodically distributed cuboidal inclusions. Then, this truncated infinite sum solution is corrected by substracting the average residual stress given by f C ijkl ε * kl where f is the volume fraction of periodic eigenstrained inclusions in the infinite medium. Such truncation and correction were considered as satisfactory enough to be used an analytical solution that has been compared to FFT results for the local stress field (see Figs. 3 and4).

Displacement components

Let us consider again the cuboidal configuration with same dimensions and same elastic moduli, eigenstrain as before. Assuming that the eigenstrain field ε * ij is due to an array of rectangular dislocation loops with general equivalent Burgers vector denoted b(b 1 , b 2 , b 3 ), this eigenstrain is responsible for a displacement vector field u(x) that can be calculated (after lengthy derivations) from the Burgers displacement formula, see [START_REF] Hirth | Theory of Dislocations[END_REF]. Therefore, the first component U 1 (x) related to the x direction writes:

u 1 (x) = 1 8π(1 -ν) 8 n=1 (-1) n {ε * 11 [(2ν -3)(C n3 log(R n + C n3 ) + C n2 log(R n + C n2 ))+ 2(1 -ν)C n1 arctan C n2 C n3 C n1 R n ] + ε * 22 [(1 -2ν)(C n2 log(R n + C n3 )) + 2ν(C n1 arctan C n2 C n3 C n1 R n -C n3 log(R n + C n2 )] + ε * 33 [(1 -2ν)(C n3 log(R n + C n2 )) + 2ν(C n1 arctan C n2 C n3 C n1 R n - C n2 log(R n + C n3 )]} (B.
2)

The eigenstrains are related to the Burgers vector components as follows:

ε * 11 = b 1 /(2a 1 ) ε * 22 = b 2 /(2a 2 ) ε * 33 = b 3 /(2a 3 )
The others components of the displacements fields are obtained by cyclic permutation. Let us note that for the numerical applications presented in this paper, the only non zero eigenstrain component is ε * 33 , i.e. only b 3 = 0.

In order to derive the analytical displacement solution to be compared with the FFT-based results, an infinite periodic distribution of cuboidal inclusions with the same periods as in the FFT-based calculations should be considered, which theoretically corresponds to an infinite sum of individual displacement contributions using Eq. (B.2). In order to reduce the computations, we first calculate in the unit cell the total displacement given by the superposition of a

  t d (x i ) derived in a consistent way with the introduced spatial discretization as shown in Fig.1.

Figure 1 .

 1 Figure 1. Definition of t d (x i ) with respect to t(x). This configuration is consistent with the spatial discretization explained in the text.

  contains the same terms as for the present discrete Green operator but these ones are raised to the power squared. A fine comparison with the paper of Brisard and Dormieux (2010)(see the exact DFT expression of the stress polarization on top right p. 666 in[START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF]) shows that our derivations in section 3.1 is the same as their derivation up to Eq. 21 but is different from Eqs. 21 to 34. This difference results from the use the Lippmann-Schwinger equation which is a convolution as in Eq. 22 while in[START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] an energy formulation based on the Hashin-Shtrikman variational principle was derived with the help of the Plancherel theorem, see their equations (

4. 1

 1 Cubic-shaped eigenstrained inclusion in a homogeneous isotropic elastic periodic medium For this case, the isotropic elastic constants of both matrix and inclusion are the following: The Young's modulus is E m = E i = 333.4GP a and the Poisson ratio is ν m = ν i = 0.26. This approximately corresponds to the room temperature elastic constants for a single crystalline Ni-based superalloy. The eigenstrain ε * ij is assumed to be present in the inclusion and uniform in it defined as: ε * ij = 0 for all i, j exceptε * 33 = 1.The inclusion shape is cubic and is centered in the unit cell. The dimensions of the cubic-shaped inclusion are 32 × 32 × 32 voxels. This inclusion is schematically represented in Fig.2with the Cartesian basis. Some specific material points are considered to examine the spatial stress variations at different locations away or near materials discontinuities. These points are denoted O, A, B, C in Fig.2and their positions are given in terms of voxel numbers from the center O of the unit cell (or the cube). Therefore, the voxel coordinates of these points are O(0, 0, 0), A(16, 0, 0), B(16, 16, 16), C(0, 16, 16). Three different specific lines denoted d 1 , d 2 , d 3 are also reported on Fig.2. These lines are specifically designed to study stress profiles along their path as a function of the distance from materials discontinuities (inclusion/matrix interfaces, junctions between two sides of the cube, edge lines, cube corners, etc.)

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Representation of the cubic-shaped inclusion. The line d 1 is defined as parallel to the x-axis through the center of the unit cell, the line d 2 is defined as parallel to the x-axis through an edge line, the line d 3 is defined as paralle to the z-axis through the center of the unit cell. Materials points O, A, B, C are some specific points used in analysis of numerical results.

Fig. 5

 5 Fig. 5 also shows the relative difference computed with the CGO. At the center of the inclusion (i.e. where there is no materials discontinuity) (Fig. 5(a)), is close to 10 -4 with the CGO which is similar to the result given by the DGO k with k = 101. However, the prediction given by the CGO is not as good as point O for the other points where materials discontinuities are present such as points B, A, C as observed on Figs. 5(b)(c)(d) respectively. For all these points, the relative difference with the analytical solution obtained with DGO 81 is much smaller than the one with CGO. It means that the proposed periodized discrete Green operator Γ D (Eq. 37) improves the accuracy of the classical FFT-based method algorithm based on the CGO in the vicinity of materials discontinuities and k = 81 reflects a good compromise to

Figure 5 .

 5 Figure 5. Relative difference computed at some specific points O (figure (a)), B (figure (b)), A (figure (c)), C (figure (d)). The specific points were defined in Fig. 2 of the cubic-shaped inclusion with CGO and DGO k (for different values of k). Solid lines represent the CGO solution and dotted lines represent the DGO k solutions.

σFigure 6 .

 6 Figure 6. The shear stress component σ 12 normalized with 2µ m (where µ m is the elastic shear modulus) for the cubic-shaped inclusion configuration with homogeneous isotropic elasticity computed with CGO, DGO 81 and FFT-solution with CGO "truncated" at the highest frequencies (Moulinec and Suquet, 1998) (named "modified CGO" in the figure).

Fig. 7 Figure 7 .

 77 Fig. 7 shows the normalized components u 3 /b and u 1 /b along the line d 3 obtained from the approximation of CGO, DGO k (for different values of k) and the analytical solution (see Appendix B). From k = 21, the values of these normalized displacements are almost identical, which shows a faster convergence of B D in comparison with the results obtained with Γ D in section 4.1. This trend was expected as discussed in section 3.

|Figure 8 .

 8 Figs. 8(a)(b) shows the values of 1 and 3 for the DGO k and the CGO methods. These values are computed at the center of the inclusion (point O) and plotted as a function of k using logarithmic scales. Regarding 1 (Fig. 8(a)), the solution given by

Figure 9 .

 9 Figure 9. Effect of the discretization level (Fourier grid resolution) on the accuracy of the displacement field component u 3 /b along d 3 for a flat cubic-shaped inclusion (Eshelbylike squared prismatic dislocation loop) given by CGO, and DGO 21 for two resolutions: 32×32×32 voxels (8×8×1 voxels for the inclusion) and 128×128×128 voxels (32×32×1 voxels for the inclusion).

Figure 10 .

 10 Figure 10. Normalized stress field component σ 11 /E m along a centered line parallel to the x-axis for the spherical inclusion with eigenstrain, given by the Eshelby's analytical solution and computed with CGO and DGO 81 . Here, the classic spatial discretization is used: x i = i d N .

Figure 11 .

 11 Figure 11. Normalized stress field components σ 11 /E m (figure (a)) and σ 33 /E m (figure (b)) along a centered line parallel to the x-axis for the spherical inclusion with eigenstrain,

Figure 12 .

 12 Figure 12. Normalized stress field components σ 11 /E m (figure (a)) and σ 33 /E m (figure (b))

Figure 13 .

 13 Figure13. The convergence error computed with the strain field criterion for both CGO and DGO 81 methods as a function of iteration number up to 50 iterations for two mechanical contrasts: R = K m /K i = 2 and R = K m /K i = 0.5.

Fig. 14 Figure 14 .

 1414 Fig. 14 represents the shear component σ 12 normalized by 2µ m computed with CGO and DGO 81 . The figure shows strong oscillations along the x-axis with the CGO.These oscillations are completely removed by the application of the proposed DGO 81 which also well matches the analytical solution especially at materials discontinuities.As a conclusion of this inhomogeneous case, the substitution of the classical continuous Green operator (CGO) by the proposed discrete Green operator (DGO 81 ) allows to remove spurious oscillations on mechanical fields and does not modify the convergence rate of the iterative scheme for the chosen mechanical contrast. In this work, the study of more important varying stiffness contrasts were not investigated as well as the optimized choice of the elastic stiffness reference medium. The two main reasons are that: (i) the present FFT-based algorithm is applied to Ni-based superalloys which exhibit a low mechanical contrast between phases, (ii) the trun-

Fig. 15

 15 Fig. 15(a) using Point Data or Fig. 15(c) using Cell Data with the visualization software ParaView[START_REF] Utkarsh | The ParaView Guide: A Parallel Visualization Application[END_REF], it is seen that the CGO produces artificial 2D wavy patterns on this stress component which are more pronounced inside the inhomogeneity and near material discontinuities. Such kind of artefacts using the CGO was

Figure 15 .

 15 Figure 15. 2D sections of the stress component σ 12 normalized by 2µ m computed with CGO and DGO 81 . Using the visualization software ParaView (Utkarsh, 2015), Point Data are used to represent the result of CGO on figure (a) and DGO 81 on figure (b); Cell Data are used to represent the result of CGO on figure (c) and DGO 81 on figure (d).

  15(b) using Point Data or Fig. 15(d) using Cell Data. The eight pairs of stress hot spots (see Fig. 15(b)) result from the discrete voxelization for this spherical inhomogeneity, which gives a different rendering using Cell Data instead of Point Data with the software ParaView (Utkarsh, 2015). A new consistent periodized discrete Green operator (abbreviated DGO k ) is proposed to replace the classic continuous modified Green operator (CGO) in Fast Fourier Transform-based methods to numerically solve the Lippmann-Schwinger equation for periodic heterogeneous media with eigenstrain fields. The advantage of this new consistent Green operator is that it is well adapted to iterative schemes based on the Lippmann-Schwinger equation as the one developed by Anglin et al.(2014). Through the different three-dimensional configurations treated in this paper considering different inclusion shapes, eigenstrains and elastic properties (cubicshaped inclusion, spherical Eshelby inclusions and spherical inhomogeneity), this periodized DGO is very effective to suppress numerical oscillations near materials discontinuities and is computed once at the beginning of the iterative algorithm by truncating the series in the expression of the DGO up to a given order k. A same type of DGO is also developed to obtain the displacement solution for the same types of materials. In this paper, it is successfully applied to the case of the displacement generated by a prismatic dislocation loop described as an Eshelby-like square inclusion with one voxel thickness. It is also observed that its convergence to accurate results devoid of oscillations is faster that the strain-based DGO due to a different Green operator.Further studies will be focused on the extension of the present DGO method for dislocation mechanics and to other more refined FFT-based schemes than fixed-point algorithms, i.e. fast gradient methods. Optimal choices of the DGO truncation order k and reference medium for heterogeneous materials containing rigid inclusions and voids will be investigated in a future work. Furthermore, the numerical method based on DGO will be applied to simulate Synchrotron high-resolution diffraction Therefore, there are two alternative descriptions of the step function t d , either with initial N points, or with 2kN points (i.e. finer discretization). T D and T D 2kN can both be expressed as infinite sums of components of T d . The values of t d computed by inverse DFT should be equal for y l = x i . Using the inverse DFT of T D to compute the step function t d can also be seen as an alternative way to compute the DFT of T D 2kN with the same accuracy but on a subset of the y l points.

  n [-2ν ln(2R n + 2C n3 ) + C n3 Rn ] n [-ln(2R n + 2C n2 ) -C n2 Rn + n [-ln(2R n + 2C n1 ) -C n1 Rn + n [-2ν arctan C n1 C n3 C n2 Rn -C n1 C n2 C n3

Acknowledgements

The authors thank G. Ribarik (Department of Materials Physics, Eötvös University, Budapest, Hungary) for useful advices. SB also thanks R. A.

Lebensohn (Los Alamos Nat. Lab., USA) and P. Suquet (LMA, CNRS Marseille, France) for fruitful discussions on discrete operators. This work was supported by the French State through the program "Investment in the future" (LabEx "DAMAS" referenced as ANR-11-LABX-0008-01).

 [START_REF] Eloh | The effect of crystal defects on 3D high-resolution diffraction peaks: A FFT-based method[END_REF]. To be well understood, such diffraction peaks need very accurate predictions of local displacement fields in the course of plastic deformation, which are only reachable with the new DGO k .

Appendix

A Equivalence between DGO k on a set of N points and CGO on a set of 2kN points

As in section 3.1, a discrete stepwise function t d (x i ) is considered with

, where i varies from 0 to N -1. With N points, the discrete Fourier transform (DFT) T D of t d is given in Eq. 19:

The DFT denoted T D was given in Eq. 28 as a function of the Fourier transform T d :

Choosing a finer discretization, the period d can be divided into 2kN points y l = ld 2kN where l varies from 0 to 2kN -1. Let us now define the DFT denoted T D 2kN as follows:

The zero shift of the y l points is chosen so that the x i points are included in the y l set. Using the inverse FT yields:

Then, from exp(j2πξ 2knN y l ) = 1, we get:

Now, T D 2kN is calculated as follows:

finite number of 27 × 27 × 27 periodically distributed cuboidal inclusions. Then, this truncated infinite sum solution is corrected by substracting the overall displacement given by f ε * ij x j where f is the volume fraction of periodic eigenstrained inclusions in the infinite medium. Such truncation and correction were considered as satisfactory enough to be used an analytical solution that has been compared to FFT results for the local displacement field (see Fig. 7).