# Fourier-based schemes for computing the mechanical response of composites with accurate local fields

François Willot

► To cite this version:

François Willot. Fourier-based schemes for computing the mechanical response of composites with accurate local fields. Comptes rendus de l'Académie des sciences. Série IIb, Mécanique, 2015, 343 (3), pp.232-245. 10.1016/j.crme.2014.12.005. hal-01096757v2

HAL Id: hal-01096757

https://minesparis-psl.hal.science/hal-01096757v2

Submitted on 13 Jul 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fourier-based schemes for computing the mechanical response of composites with accurate local fields

François Willot

Mines ParisTech, PSL Research University, Centre for Mathematical Morphology, 35, rue St Honoré, 77300 Fontainebleau, France

Nota Bene: this is a “post-print” version of the published paper in which one typo (Eq. 40, page 11) has been corrected. Correction highlighted in blue; results are unchanged.

## Abstract

We modify the Green operator involved in Fourier-based computational schemes in elasticity, in 2D and 3D. The new operator is derived by expressing continuum mechanics in terms of centered differences on a rotated grid. Use of the modified Green operator leads, in all systems investigated, to more accurate strain and stress fields than using the discretizations proposed by Moulinec and Suquet (1994) or Willot and Pellegrini (2008). Moreover, we compared the convergence rates of the “direct” and “accelerated” FFT schemes with the different discretizations. The discretization method proposed in this work allows for much faster FFT schemes with respect to two criteria: stress equilibrium and effective elastic moduli.

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

## 1. 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 [2]. In FFT methods, the microstructure is defined by 2D or 3D images and the local stress and

Email address: francois.willot@ensmp.fr (François Willot)

Tel.: +33 1 64 69 48 07.
strain tensors are computed along each pixel or “voxel” in the image. Coupled with automatic or semi-automatic image segmentation techniques [3], this allows for the computation of the mechanical response of a material directly from experimental acquisitions, like focused ion beam or 3D microtomography techniques [4]. The latter often deliver images containing billions of voxels, for which FFT methods are efficient [5, 6]. This allows one to take into account representative volume elements of materials which are multiscale by nature such as concrete or mortar [7]. 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 [8, 9] but also as commercial ones [10].

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 [11] corresponds, under one technical assumption, to a particular choice of approximation space and optimization method [12] (see also [13]). 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 [14, 15].

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 [15, 13] 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 [13, 16]. In the context of conductivity, accurate local fields and improved convergence rates have also been obtained from modified Green operators based on finite-differences [17]. These results follow earlier works where continuum mechanics are expressed by centered [18, 19] or “forward and backward” finite differences [20].

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. (3). In Sec. (4), a general formulation of the Green operator is derived that incorporates methods in [20], and a new discretization scheme is proposed. The accuracy of the local stress and strain fields are examined.
ined in Sec. (5) whereas the convergence rates of the various FFT methods are investigated in Sec. (6). We conclude in Sec. (7).

## 2. Microstructure and material elastic response

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

$$ \sigma_{ij}(\mathbf{x}) = C_{ij,kl}(\mathbf{x})\varepsilon(\mathbf{x}), \quad \partial_i\sigma_{ij}(\mathbf{x}) \equiv 0, \quad \varepsilon_{ij}(\mathbf{x}) = (1/2)[\partial_i u_j(\mathbf{x}) + \partial_j u_i(\mathbf{x})], \quad (1) $$

where $\varepsilon(\mathbf{x})$ is the strain field, $\sigma(\mathbf{x})$ the stress field, $u(\mathbf{x})$ the displacement vector field, $C(\mathbf{x})$ the local elasticity tensor and $\mathbf{x}$ is a point in $\Omega$. Tensorial components refer to a system of Cartesian coordinates $(\mathbf{e}_1; \mathbf{e}_2)$ in 2D and $(\mathbf{e}_1; \mathbf{e}_2; \mathbf{e}_3)$ in 3D. The material has an isotropic local elastic response that reads:

$$ C_{ij,kl}(\mathbf{x}) = \lambda(\mathbf{x})\delta_{ij}\delta_{kl} + \mu(\mathbf{x})(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}), \quad (2) $$

where $\delta$ is the Kronecker symbol and $\lambda(\mathbf{x})$ and $\mu(\mathbf{x})$ are constant-per-phase Lamé's first and second coefficients. The local bulk modulus $\kappa = \lambda + (2/d)\mu$ and the elastic moduli take on values:

$$ \lambda(\mathbf{x}) = \lambda^{\alpha}, \quad \kappa(\mathbf{x}) = \kappa^{\alpha}, \quad \mu(\mathbf{x}) = \mu^{\alpha}, $$

in phase $\alpha$. For simplicity, we restrict ourselves to binary media and, by convention, $\alpha = 1$ is the matrix and $\alpha = 2$ the inclusions. Hereafter, we fix Poisson's ratios in each phase to $\nu^1 = \nu^2 = 0.25$ so that, in 3D and 2D [21], we have $\mu^\alpha/\kappa^\alpha = 0.6$. The contrast of properties $\chi$ reads:

$$ \chi = \frac{\kappa^2}{\kappa^1} = \frac{\mu^2}{\mu^1} = \frac{\lambda^2}{\lambda^1}, \quad (3) $$

where $0 \le \chi \le \infty$. In the matrix, we also fix $\kappa^1 = 1$ ($d = 2$ or $3$), $\mu^1 = 0.6$ ($d = 2$ or $3$), $\lambda^1 = 0.4$ ($d = 2$), $\lambda^1 = 0.6$ ($d = 3$), so that the local properties of the material are parametrized by one unique variable, the contrast of properties $\chi$. In 3D, the Young modulus is $E^1 = 3/2$ in the matrix and $E^2 = 3\chi/2$ in the inclusion. The medium is porous when $\chi = 0$ and rigidly-reinforced when $\chi = \infty$.

Periodic boundary conditions are applied with the material subjected to an overall strain loading $\bar{\varepsilon}$:

$$ \sigma_{ij}(\mathbf{x})n_j(\mathbf{x}) - \# (\mathbf{x} \in \partial\Omega), \quad \langle\varepsilon_{kl}(\mathbf{x})\rangle = \overline{\varepsilon_{kl}}, \quad (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 $\tilde{C}$ is computed from:

$$ \langle \sigma_{ij}(\mathbf{x}) \rangle = \tilde{C}_{ijkl} \overline{\varepsilon}_{kl}. \qquad (5) $$

### 3. 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 [22]:

$$ \tau_{ij}(\mathbf{x}) = \sigma_{ij}(\mathbf{x}) - C_{ijkl}^{0} \varepsilon_{kl}(\mathbf{x}), \quad \varepsilon_{ij}(\mathbf{x}) = \overline{\varepsilon}_{ij} - \int_{\mathbf{x}'} d^{d}(\mathbf{x}' - \mathbf{x}) G_{ij,kl}(\mathbf{x}' - \mathbf{x}) \tau_{kl}(\mathbf{x}'), \quad (6) $$

where we have introduced a homogeneous “reference” elasticity tensor $\tilde{C}^0$ and its associated polarization field $\tau$ and Green operator $\mathcal{G}$. In the above we assume $\langle \mathcal{G} \rangle = 0$ so that $\overline{\varepsilon} = \langle \varepsilon(\mathbf{x}) \rangle$ holds. The Green operator has, in the Fourier domain, the closed form [23]:

$$ G_{ij,kl}(\mathbf{q}) = \left\{ q_i \left[ q_m C_{m,j,k,n}^0 q_n \right]^{-1} q_l \right\}_{\text{sym}}, \qquad (7) $$

where $\mathbf{q} \neq 0$ are the Fourier wave vectors and the subscript $\text{sym}$ indicates minor symmetrization with respect to the variables $(i, j)$ and $(k, l)$. Hereafter, we assume that $\tilde{C}^0$ is a symmetric, positive-definite, isotropic tensor defined by its bulk $(\kappa^0)$ and shear $(\mu^0)$ moduli, or Lamé coefficient $(\lambda^0)$. Accordingly, when $\mathbf{q} \neq 0$, the Green operator $\mathcal{G}$ is also symmetric definite and we have:

$$ G_{ij,kl}(\mathbf{q}) = \frac{1}{\mu^0} \left[ \left( \frac{q_i q_l}{|\mathbf{q}|^2} \delta_{jk} \right)_{\text{sym}} - \frac{\lambda^0 + \mu^0}{\lambda^0 + 2\mu^0} \frac{q_i q_j q_k q_l}{|\mathbf{q}|^4} \right]. \qquad (8) $$

The “direct scheme” [11] consists in applying Eqs. (6) iteratively as:

$$ \varepsilon^{k=0} \equiv \overline{\varepsilon}, \quad \varepsilon^{k+1} = \overline{\varepsilon} - \mathcal{G} * (\sigma - \tilde{C}^0 : \varepsilon^k) \quad (k \ge 0). \qquad (9) $$

In Moulinec and Suquet’s method, the convolution product $(*)$ above is computed as an algebraic product in the Fourier domain, making use of (8). Discrete Fourier transforms are used to switch between the space $\Omega$ and Fourier domain $\mathcal{F}$. This amounts to representing the strain field as a trigonometric polynomial [12] of the form:

$$ \varepsilon(\mathbf{x}) = \frac{1}{L^d} \sum_{\mathbf{q} \in \mathcal{F}} \varepsilon(\mathbf{q}) e^{i\mathbf{q} \cdot \mathbf{x}}, \qquad (10) $$
where $\varepsilon(\boldsymbol{q})$ is the discrete Fourier transform of $\varepsilon(\boldsymbol{x})$. Accordingly:

$$ \varepsilon(\boldsymbol{q}) = \sum_{\boldsymbol{x} \in \Omega} \varepsilon(\boldsymbol{x}) e^{-i\boldsymbol{q} \cdot \boldsymbol{x}}. \qquad (11) $$

Similar forms are used for the stress and displacement fields. In practice, the domain $\Omega$ is discretized on a square or cubic grid of $L^d$ voxels and the operator $\mathbb{G}(\boldsymbol{q})$ in (7) is evaluated along equispaced Fourier modes:

$$ q_i = \frac{2\pi m_i}{L}, \qquad m_i = \begin{cases} 1 - (L/2), \dots, L/2, & (L \text{ even}), \\ -(L-1)/2, \dots, (L-1)/2, & (L \text{ odd}). \end{cases} \qquad (12) $$

As noted in [24], when $L$ is even, the relation:

$$ \mathbb{G}(\boldsymbol{q})^* = \mathbb{G}(-\boldsymbol{q}), \qquad (13) $$

where $\mathbb{G}(\boldsymbol{q})^*$ is the complex conjugate of $\mathbb{G}(\boldsymbol{q})$, is not verified when one of the components of $\boldsymbol{q}$ is equal to the highest frequency $q_i = \pi$ (i.e. $m_i = L/2$). As a consequence, the backward Fourier transform of $\mathbb{G}(\boldsymbol{q})\boldsymbol{\tau}(\boldsymbol{q})$ used to compute the strain field has non-zero imaginary part even if $\boldsymbol{\tau}(\boldsymbol{x})$ is purely real. To fix this problem, we follow [24] and set:

$$ \mathbb{G}(\boldsymbol{q}) = (\mathbb{C}^0)^{-1}, \qquad \text{if } q_i = \pi \text{ for some } i. \qquad (14) $$

This choice enforces $\boldsymbol{\sigma}(\boldsymbol{q}) = 0$ at the concerned Fourier modes. In doing so, the strain field $\varepsilon$ 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 $\mathbb{G}(\boldsymbol{q})$ with:

$$ \frac{\mathbb{G}(\boldsymbol{q})^* + \mathbb{G}(-\boldsymbol{q})}{2}, \qquad (15) $$

when one of the components of $\boldsymbol{q}$ equals $\pi$, which enforces $\varepsilon(\boldsymbol{q}) = 0$ at the highest modes. Numerical experiments indicate that the choice for $\mathbb{G}(\boldsymbol{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 $\mathbb{G}(\boldsymbol{q}) = (\mathbb{C}^0)^{-1}$. Furthermore, in the 2D example studied in this work, the choice $\mathbb{G}(\boldsymbol{q}) = (\mathbb{C}^0)^{-1}$ led to smaller oscillations, consistently with observations in [24]. 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” [25] and “augmented Lagrangian” [26] methods. In this work, we use the extension of the accelerated scheme to elasticity [26, 27]:

$$ \varepsilon^{k+1} = \varepsilon^k + 2(C + C^0)^{-1} : C^0 : [\bar{\varepsilon} - \varepsilon^k - G^* (C : \varepsilon^k - C^0 : \varepsilon^k)]. \quad (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 [25, 27]:

$$ \kappa^0 = \sqrt{\kappa^1 \kappa^2}, \quad \mu^0 = \sqrt{\mu^1 \mu^2}. \quad (17) $$

For the direct scheme, upper bounds on the eigenvalues of the Green operator suggest the choice [27]:

$$ \kappa^0 = \beta(\kappa^1 + \kappa^2), \quad \mu^0 = \beta(\mu^1 + \mu^2), \quad (18) $$

with $\beta = 1/2$.

## 4. Discretization and approximation space

In this section, we derive the expression of a modified Green operator $G'$ that replaces $G$ defined in (8) and (12). We give it in a form that includes previously-proposed modified operators [20] and also introduce a new one.

### 4.1. 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. (10) and (11) 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 (10) and (11). In [20], this results in the following form:

$$ k_i^*(q)\sigma_{ij}(q) = 0, \quad \varepsilon_{ij}(q) = (1/2)[k_i(q)u_j(q) + k_j(q)u_i(q)], \quad (19) $$

where $k$ and $k^*$ represent “discrete” gradient and divergence operators, respectively. In the centered scheme, one takes $k$ equal to:

$$ k_i^C(q) = i \sin(q_i), \quad (20) $$
whereas in scheme [20], one chooses, for $k$:

$$k_i^W(\boldsymbol{q}) = e^{iq_i} - 1. \qquad (21)$$

These expressions correspond, respectively, to the centered scheme:

$$\partial_j \sigma_{ij}(\mathbf{x}) \approx \frac{\sigma_{ij}(\mathbf{x} + \mathbf{e}_j) - \sigma_{ij}(\mathbf{x} - \mathbf{e}_j)}{2}, \qquad \partial_j u_i(\mathbf{x}) \approx \frac{u_i(\mathbf{x} + \mathbf{e}_j) - u_i(\mathbf{x} - \mathbf{e}_j)}{2}, \qquad (22)$$

and to the forward-and-backward difference scheme:

$$\partial_j \sigma_{ij}(\mathbf{x}) \approx \sigma_{ij}(\mathbf{x}) - \sigma_{ij}(\mathbf{x} - \mathbf{e}_j), \qquad \partial_j u_i(\mathbf{x}) \approx u_i(\mathbf{x} + \mathbf{e}_j) - u_i(\mathbf{x}). \qquad (23)$$

Using (19), the resulting Green operator reads:

$$G'_{ij,kl}(\boldsymbol{q}) = \left\{ k_i(\boldsymbol{q}) [k_m(\boldsymbol{q}) C^0_{m,j,Kn} k_n^*(\boldsymbol{q})]^{-1} k_l^*(\boldsymbol{q}) \right\}_{\text{sym}}, \qquad (24)$$

which is homogeneous in $k$, so that, with $r_i = k_i/|k|$:

$$G'_{ij,kl}(\boldsymbol{q}) = \frac{(\lambda^0 + 2\mu^0) (r_i r_l^* \delta_{jk})_{\text{sym}} - \lambda^0 \text{Re}(r_i r_j^*) \text{Re}(r_k r_l^*) - \mu^0 r_i r_j (r_k r_l)^*}{\mu^0 [2(\lambda^0 + \mu^0) - \lambda^0 |r_1^2 + r_2^2|^2]}, \qquad (25)$$

where $\text{Re}(\cdot)$ 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 $\mathbb{G}'$ is denoted by $\mathbb{G}^C$ when $k = k^C$ and by $\mathbb{G}^W$ when $k = k^W$. Note that the operator $\mathbb{G}$ in (8) is recovered from (25) by setting $k_i(\boldsymbol{q}) = iq_i$. Now the Green operator $\mathbb{G}'$ is complex and follows the minor and major symmetries:

$$G'_{ij,kl}(\boldsymbol{q}) = G'_{ji,kl}(\boldsymbol{q}) = G'_{ij,lk}(\boldsymbol{q}), \qquad G'_{ij,kl}(\boldsymbol{q}) = [G'_{kl,ij}(\boldsymbol{q})]^*. \qquad (26)$$

Furthermore we have:

$$\mathbb{G}^{C,W}(\boldsymbol{q})^* = \mathbb{G}^{C,W}(-\boldsymbol{q}), \qquad (27)$$

for the schemes (22) and (23), including when $q_i = \pi$ and when $L$ is even, because $k_i^{C,W}$ is real at the frequency $q_i = \pi$. Therefore, the fix (14) in Sec. (3) is not necessary. However, a problem of a different nature arises when using the centered scheme (20) when $L$ is even. Eq. (25) does not define the Green operator $\mathbb{G}^C$ at the three frequencies $\boldsymbol{q} = (\pi; 0), (0; \pi)$ and $(\pi; \pi)$, for which $k^C(\boldsymbol{q}) = 0$. This is because the second equation in (22) has in general non-unique solutions for the displacement field $\boldsymbol{u}(\mathbf{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_m^1(\mathbf{x}) = \delta_{mn}, \quad v_m^2(\mathbf{x}) = (-1)^j \delta_{mn}, \quad v_m^3(\mathbf{x}) = (-1)^j \delta_{mn}, \quad v_m^4(\mathbf{x}) = (-1)^{i+j} \delta_{mn}, \\ \mathbf{x} = \left(\frac{i}{L}; \frac{j}{L}\right), \quad i, j = 1 - \frac{L}{2}, \dots, \frac{L}{2}, \quad n = 1, 2. \quad (28)$$

The operator $\mathbb{G}^C$ remains finite when $q$ approaches one of the modes $(\pi; 0), (0; \pi)$ or $(\pi; \pi)$, but can not be continuously extended at these modes. To fix this problem, we set, for the centered scheme:

$$G_{ijkl}^{\mathbb{G}}(\mathbf{q}) = 0, \quad \text{if } L \text{ is even and } (q_i = 0 \text{ or } q_i = \pi) \text{ for all } i, \quad (29)$$

which enforces $\varepsilon(\mathbf{q}) = 0$ at the highest frequencies. The strain field $\varepsilon(\mathbf{x})$ is accordingly admissible, and the stress field $\sigma(\mathbf{x})$ is divergence-free, in the sense of (22). We explored the alternate choice $\mathbb{G}^C(\mathbf{q}) = (\mathbb{C}^0)^{-1}$ in (29). Almost identical convergence rates and oscillations were observed for the two options, however the choice $\mathbb{G}^C(\mathbf{q}) = (\mathbb{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 $\mathbb{G}^W$ at high frequencies since $\mathbf{k}^W \neq 0$ when $\mathbf{q} \neq 0$.

By substituting $\mathbb{G}$ with $\mathbb{G}^C$ or $\mathbb{G}^W$ in (9) and (16), we derive “direct” and “accelerated” schemes that solve (22) or (23). In the limit of very fine resolution, we have $\mathbb{G}^{C,W}(\mathbf{q}) \approx \mathbb{G}(\mathbf{q})$ when $\mathbf{q} \to 0$, which guarantees that the strain and stress fields do not depend on the employed discretization. This property holds for any choice of $\mathbf{k}$ such as $\mathbf{k} \sim i\mathbf{q}$ when $\mathbf{q} \to 0$.

On the one hand, derivatives are estimated more locally in the forward-and-backward scheme (21) than in the centered scheme (20), which is important along interfaces. On the other hand, the forward-and-backward scheme does not treat symmetrically the two angle bisectors $\mathbf{e}_1 + \mathbf{e}_2$ and $\mathbf{e}_1 - \mathbf{e}_2$ [20]. 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 [20], are not explored in this work. The former 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 (23). 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°-rotated basis:

$$ \mathbf{f}_1 = \frac{\mathbf{e}_1 + \mathbf{e}_2}{\sqrt{2}}, \quad \mathbf{f}_2 = \frac{\mathbf{e}_2 - \mathbf{e}_1}{\sqrt{2}}, \qquad (30) $$

by:

$$ u_i = R_{iI}u_I, \quad \varepsilon_{ij} = R_{iI}\varepsilon_{IJ}R'_{jI}, \quad \sigma_{ij} = R_{iI}\sigma_{IJ}R'_{jI}, \quad R_{iJ} = \frac{1 - 2\delta_{i1}\delta_{j2}}{\sqrt{2}}, \qquad (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):

$$ \sigma_{IJ}(\mathbf{x}) = C_{IJKL}(\mathbf{x}) \varepsilon_{KL}(\mathbf{x}), \qquad (32a) $$

$$ \sigma_{I1}(\mathbf{x}) - \sigma_{I1}(\mathbf{x} - \sqrt{2}\mathbf{f}_1) + \sigma_{I2}(\mathbf{x} + \frac{\mathbf{f}_2 - \mathbf{f}_1}{\sqrt{2}}) - \sigma_{I2}(\mathbf{x} - \frac{\mathbf{f}_1 + \mathbf{f}_2}{\sqrt{2}}) = 0, \qquad (32b) $$

$$ \varepsilon_{KL}(\mathbf{x}) = \frac{1}{2\sqrt{2}} \left[ u_K(\mathbf{x} + \frac{\mathbf{f}_L}{\sqrt{2}}) - u_K(\mathbf{x} - \frac{\mathbf{f}_L}{\sqrt{2}}) + u_L(\mathbf{x} + \frac{\mathbf{f}_K}{\sqrt{2}}) - u_L(\mathbf{x} - \frac{\mathbf{f}_K}{\sqrt{2}}) \right] (32c) $$

where $\mathbf{x}$ lie at the centers of the pixels and $\mathbf{x} \pm \mathbf{f}_I / \sqrt{2}$ lie at the corners. Expressing back (32) in the original Cartesian grid ($\mathbf{e}_1, \mathbf{e}_2$) and applying the backward discrete Fourier transform (10) we arrive again at (19) with the following expression for $\mathbf{k}$:

$$ k_i^R(\mathbf{q}) = \frac{i}{2} \tan\left(\frac{q_i}{2}\right) (1 + e^{i q_1}) (1 + e^{i q_2}). \qquad (33) $$

We denote by $\mathbb{G}^R$ the corresponding Green operator, derived by substituting $\mathbf{k} = \mathbf{k}^R$ in (25). The operator $\mathbb{G}^R$ is real and also verifies:

$$ \mathbb{G}^R(\mathbf{q}) = \mathbb{G}^R(-\mathbf{q}). \qquad (34) $$

However, when $L$ is even, $\mathbf{k}^R = 0$ when $\mathbf{q} = (\pi; \pi)$ and $\mathbb{G}^R$ is not defined by (25) at this frequency. Again, this is because (32c) gives the displacement field up to linear combinations of the 4 independent fields $v_{1,2}^{1,4}$ (see 28). Accordingly we set:

$$ \mathbb{G}^R(\mathbf{q}) = 0, \quad \text{when } L \text{ is even, } d = 2, q_1 = q_2 = \pi, \qquad (35) $$

which enforces strain compatibility and stress equilibrium, in the sense of (32).
<figure><img src="image_2.png" /><figcaption>Figure 1: A pixel with edges parallel to the Cartesian axis ($e_1$; $e_2$). Superimposed: 45°-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).</figcaption></figure>

## 4.2. Three dimensions

We follow the same methodology in 3D. The equilibrium and strain admissibility conditions (19) are unchanged, as well as the expression for the vectors $k^{C,W}$ in (20) and (21) resulting from (22) and (23). In 3D, we also extend (33) as:

$$k_i^R(\boldsymbol{q}) = \frac{i}{4} \tan\left(\frac{q_i}{2}\right) (1 + e^{iq_1})(1 + e^{iq_2})(1 + e^{iq_3}), \quad (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 $\varepsilon_{11}$ and $\varepsilon_{12}$:

$$\varepsilon_{11}(\boldsymbol{x}) \approx \frac{1}{4} \sum_{\substack{m=\pm 1 \\ n=\pm 1}} \left[ u_1 \left( \boldsymbol{x} + \frac{\boldsymbol{e}_1 + m\boldsymbol{e}_2 + n\boldsymbol{e}_3}{2} \right) - u_1 \left( \boldsymbol{x} - \frac{\boldsymbol{e}_1 - m\boldsymbol{e}_2 - n\boldsymbol{e}_3}{2} \right) \right], (37a)$$

$$\varepsilon_{12}(\boldsymbol{x}) \approx \frac{1}{8} \sum_{\substack{m=\pm 1 \\ n=\pm 1}} \left[ u_2 \left( \boldsymbol{x} + \frac{\boldsymbol{e}_1 + m\boldsymbol{e}_2 + n\boldsymbol{e}_3}{2} \right) - u_2 \left( \boldsymbol{x} - \frac{\boldsymbol{e}_1 - m\boldsymbol{e}_2 - n\boldsymbol{e}_3}{2} \right) \right. \\ \left. \qquad + u_1 \left( \boldsymbol{x} + \frac{\boldsymbol{e}_2 + m\boldsymbol{e}_1 + n\boldsymbol{e}_3}{2} \right) - u_1 \left( \boldsymbol{x} - \frac{\boldsymbol{e}_2 - m\boldsymbol{e}_1 - n\boldsymbol{e}_3}{2} \right) \right], \quad (37b)$$

where $x$ lie at the center of a voxel. The expression for the strain component $\varepsilon_{22}$ (resp. $\varepsilon_{33}$) is obtained after exchanging the indicia 1 and 2 (resp. 1 and 3) in (37a). The component $\varepsilon_{23}$ (resp. $\varepsilon_{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:

$$ \partial_i \sigma_{i1}(\mathbf{x}) \approx \sum_{\substack{m=\pm 1 \\ n=\pm 1}} \left[ \begin{aligned} & \sigma_{11} \left( \mathbf{x} + \frac{\mathbf{e}_1 + m\mathbf{e}_2 + n\mathbf{e}_3}{2} \right) - \sigma_{11} \left( \mathbf{x} - \frac{\mathbf{e}_1 - m\mathbf{e}_2 - n\mathbf{e}_3}{2} \right) \\ & + \sigma_{12} \left( \mathbf{x} + \frac{\mathbf{e}_2 + m\mathbf{e}_1 + n\mathbf{e}_3}{2} \right) - \sigma_{12} \left( \mathbf{x} - \frac{\mathbf{e}_2 - m\mathbf{e}_1 - n\mathbf{e}_3}{2} \right) \\ & + \sigma_{13} \left( \mathbf{x} + \frac{\mathbf{e}_3 + m\mathbf{e}_1 + n\mathbf{e}_2}{2} \right) - \sigma_{13} \left( \mathbf{x} - \frac{\mathbf{e}_3 - m\mathbf{e}_1 - n\mathbf{e}_2}{2} \right) \end{aligned} \right], \quad (38) $$

where $\mathbf{x}$ lie at one of the edges of a voxel. The components $\partial_i \sigma_{i2}$ and $\partial_i \sigma_{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 ($\mathbf{e}_1$; $\mathbf{e}_2$; $\mathbf{e}_3$).

In 3D, Eq. (24) yields, for the Green operator:

$$ G'_{ijkl}(\mathbf{q}) = \frac{\left(\lambda^0 + 2\mu^0\right) \left(r_i r_l^* \delta_{jk}\right)_{\text{sym}} + \lambda^0 \left[\left(r_i r_l^* s_{jk}\right)_{\text{sym}} - \text{Re}\left(r_i r_j^*\right) \text{Re}\left(r_k r_l^*\right)\right] - \mu^0 r_i r_j r_k^* r_l^*}{\mu^0 \left[2(\lambda^0 + \mu^0) - \lambda^0 \left|r_1^2 + r_2^2 + r_3^2\right|^2\right]} \quad (39) $$

where again $r_i = k_i/|k|$ and $s$ is the symmetric second-order tensor:

$$ s_{jj} = 4 \operatorname{Im}(r_i r_k^*)^2, \quad s_{ji} = -4 \operatorname{Im}(r_k r_j^*) \operatorname{Im}(r_k r_i^*) \quad [\text{instead of } s_{jk} = -4\dots], \quad i \neq j \neq k \neq i, \quad (40) $$

with $\operatorname{Im}(\cdot)$ the imaginary part of the enclosed complex quantity. Like in 2D, the operator $\mathcal{G}'$ follows minor and major symmetries (26).

Again, the operators $\mathcal{G}^C$, $\mathcal{G}^W$ and $\mathcal{G}^R$ are derived using the expression for $\mathcal{G}'$ in (39) with $\mathbf{k} = \mathbf{k}^C, \mathbf{k}^W$ and $\mathbf{k}^R$, respectively. The symmetries (27) and (34) are verified in 3D as well. But again, a special treatment is needed for $\mathcal{G}^C$ and $\mathcal{G}^R$ when $L$ is even, at the modes $\mathbf{q}$ for which $\mathbf{k}^{C,R}(\mathbf{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 $\mathcal{G}^{C,R} = 0$ at these frequencies. More precisely, we apply (29) when $d = 3$ and, for the rotated scheme:

$$ \mathcal{G}^R(\mathbf{q}) = 0, \quad \text{if } L \text{ is even, } d = 3 \text{ and } q_i = q_j = \pi \text{ with } i \neq j. \quad (41) $$

The operators $\mathcal{G}^C$, $\mathcal{G}^W$ and $\mathcal{G}^R$ are, in 2D and 3D, periodic functions where, contrary to $\mathcal{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 [17].
We also expect higher accuracy for the local fields when employing $\mathbb{G}^R$ rather than the other discrete operators $\mathbb{G}^W$ and $\mathbb{G}^C$. First, the operator $\mathbb{G}^R$ is based on centered differences which are more precise than forward and backward differences, used in $\mathbb{G}^W$. Second, derivatives are evaluated more locally when using $\mathbb{G}^R$ rather than $\mathbb{G}^C$. Indeed, the latter are computed at points separated by 2 voxels for $\mathbb{G}^C$ instead of $\sqrt{2}$ (in 2D) or $\sqrt{3}$ voxels (in 3D) for $\mathbb{G}^R$.

The above considerations guided the choice for the discretization schemes (37) and (38), leading to $k^R$ and $\mathbb{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 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 $\mathbb{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 $\mathbb{G}$ with $\mathbb{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.

# 5. Local strain and stress fields accuracy

## 5.1. Two-dimensions

Hereafter we consider the 2D ‘four-cell’ microstructure, where the periodic domain $\Omega$ is divided into 4 identical squares of surface fraction 25%. Its non-trivial 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 $\chi = 10^3$. The material is subjected to the macroscopic strain loading:

$$ \overline{\varepsilon}_{ij} = \frac{1}{2} (\delta_{i1}\delta_{j2} + \delta_{j1}\delta_{i2}). \qquad (42) $$

We determine the strain and stress fields predicted by FFT schemes when using the Green operators $\mathbb{G}$ or $\mathbb{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 \times 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 $\sigma_{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 [17] 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$.

## 5.2. *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 $\Omega$ 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 $\chi = 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.
<figure><img src="image_3.png" /><figcaption>Figure 2: Elementary periodic domain $\Omega = [1/2; 1/2]^2$ containing a square inclusion with elastic moduli $\mu^2, \kappa^2$ (top-left, shown in white) embedded in a matrix (shown in gray) with elastic moduli $\mu^1, \kappa^1$.</figcaption></figure>

A 2D section of the stress component $\sigma_{12}(x)$ is represented in Fig. (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 $\mathcal{G}$ and $\mathcal{G}^{C,W,R}$ are close to one another. However, stress patterns near the corners of the inclusion are less pronounced with $\mathcal{G}$ than with the other methods. At smaller resolutions $L = 256$ and $L = 512$, the stress fields predicted by $\mathcal{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 $\mathcal{G}$ and $\mathcal{G}^C$. The two methods produce artificial patterns directed vertically and horizontally, close to the interface. Conversely, the fields produced by $\mathcal{G}^W$ and $\mathcal{G}^R$ have the smallest oscillations, but that of $\mathcal{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 $\mathcal{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.
<figure><img src="image_4.png" /><figcaption>Figure 3: Stress component $\sigma_{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).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: 2D section of the stress component $\sigma_{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.</figcaption></figure>
### 5.3. 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 $\Omega$ 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 $\chi = 10^{-4}$. We compute the effective elastic modulus $\tilde{C}_{11,11}$ produced by either $\mathbb{G}$ or $\mathbb{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 \cdot 10^{-10}$. Results are shown in Fig. (5) and are compared with the analytical estimate in [28]. When the resolution increases, the effective elastic modulus $\tilde{C}_{11,11}$ increases up to a limit value that we estimate to about $1.208 \pm 0.001$, for all schemes. As observed in other studies [5], very large systems are needed to compute this estimate at a high precision.

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

## 6. Convergence rate

### 6.1. Convergence rate with respect to stress equilibrium

In this section, we estimate the rates of convergence of the direct and accelerated schemes DS, DSC, W, R, AS and ASC, W, R, that use the various Green operators. All schemes enforce stress equilibrium at convergence only, therefore we follow [27] and consider a criterion based on the $L^2$-norm:

$$ \eta = \frac{\|\mathrm{div}(\boldsymbol{\sigma})\|_2}{\|\langle\boldsymbol{\sigma}\rangle\|} = \frac{1}{\|\langle\boldsymbol{\sigma}\rangle\|} \sqrt{\frac{1}{|\Omega|} \int_{\Omega} d^d\boldsymbol{x} |\mathrm{div}(\boldsymbol{\sigma})|^2} = \frac{1}{\|\langle\boldsymbol{\sigma}\rangle\|} \sqrt{\sum_q |\boldsymbol{k}(q) \cdot \boldsymbol{\sigma}(q)|^2}, \quad (43) $$

where $\eta \ll 1$ is the precision and the normalizing factor $\|\langle\boldsymbol{\sigma}\rangle\|$ is the Frobenius norm:

$$ \|\langle\boldsymbol{\sigma}\rangle\|^2 = \sum_{i,j} \langle\sigma_{ij}(\boldsymbol{x})\rangle^2. $$
<figure><img src="image_6.png" /><figcaption>Figure 5: Apparent elastic modulus $\tilde{C}_{11,11}$ estimated by FFT methods using the Green operators $\mathbb{G}$ and $\mathbb{G}^{C,W,R}$ (black and red), at increasing resolution $L$. Orange: estimate in [28]; violet: estimate of the asymptotic effective modulus $\tilde{C}_{11,11}$ using FFT data.</figcaption></figure>

In (43) we set $\boldsymbol{k} = \boldsymbol{k}^{C,W,R}$ for the schemes using $\mathbb{G}^{C,W,R}$ and $\boldsymbol{k} = i\boldsymbol{q}$ when using the Green operator $\mathbb{G}$, so that $\boldsymbol{k} \cdot \boldsymbol{\sigma}(\boldsymbol{q})$ is the divergence of the stress field in the Fourier domain, estimated according to the various discretization schemes.

We now estimate the convergence rates on a random microstructure. In the following, the domain $\Omega$ 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% [29]. 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 $\Omega$.

Taking $v^0 = v^1 = v^2 = 0.25$ for the reference Poisson ratio, we compute numerically the number of iterations $N(E^0)$ required to reach the precision $\eta \le 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 $\chi = 10^{-5}$ and the various accelerated schemes AS and $\text{AS}_{C,W,R}$ (Fig. 6). Within the range $0 < E^0 \le 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 $\text{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 \approx 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 $v^0 = 0.25 \pm 0.01$ and $0.25 \pm 0.05$ for various values of $E^0$ with $\chi = 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 $\chi = 10^{-2}$. Therefore, in the following, we fix the Poisson ratio to $v^0 = 0.25$ for the reference tensor, for all schemes and all contrast of properties $\chi$. 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 $\eta \le 10^{-8}$. We choose $E^0 = 0.51(E^1 + E^2)$ for schemes DS, DS$_{C,W,R}$ and $E^0 = \sqrt{E^1 E^2}$ for schemes AS and AS$_{C,W,R}$ as initial guess for $E^0$, suggested by (18) and (17). At each step, we determine if $E^0$ is to be increased or decreased, by comparing $N(E^0)$ with $N(E^0+\delta E^0)$ where $\delta E^0 = 0.01 E^0$. It frequently happens that $N(E^0) = N(E^0 + \delta E^0)$. In that case, we compare the values of the precision $\eta$ after $N(E^0)$ iterations and follow the direction that minimizes $\eta$. 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 $\chi = 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 $\chi = 10^{-5}, 10^5$, but not for scheme AS$_C$ with $\chi = 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. (6). 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. (7). For the direct scheme, the optimal reference follows (18) with $0.5003 \le \beta \le 0.509$, independently of the Green operator used. Values of $\beta$ smaller than $1/2$ lead to non-converging 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 $\mathbb{G}$, the choice (17) is optimal except in the region $\chi \le 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 \approx E_1^0 + \sqrt{\chi}, \qquad (44)$$

with $E_1^0 = 0.07$ for AS$_C$ and $E_1^0 = 0.12$ for AS$_R$. Similar behavior has been
<figure><img src="image_7.png" /><figcaption>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 $\eta = 10^{-8}$ is reached. The microstructure is a Boolean model of quasi-porous spheres with $\chi = 10^{-5}$.</figcaption></figure>

observed numerically in [17], in the context of conductivity. For the scheme AS$_W$ that uses $\mathbb{G}^W$, the optimal choice for $E^0$ follows the same pattern as above with $E_1^0 = 0.7$ except in the region $10^{-4} \le \chi \le 10^{-1}$. This behavior is an effect of the presence of two local minima, similar to that shown in Fig. (6) for $\chi = 10^{-5}$.

Convergence rates, computed with optimized reference, are represented in Fig. (8) as a function of the contrast, in log-log scale. Results for $\chi = 0$ (strictly porous media) have been included in the same graph (left point). As is well-established [26, 30], the number of iterations in the direct scheme DS scales as $\chi$ when $\chi \gg 1$ and $1/\chi$ when $\chi \ll 1$. For the accelerated scheme AS, the number of iterations is smaller and follows $\sqrt{\chi}$ when $\chi \gg 1$ and $1/\sqrt{\chi}$ when $\chi \ll 1$, with one exception. At very high contrast of properties $\chi < 10^{-6}$, including at $\chi = 0$, convergence is reached after a finite number of iterations, about 1,300. This particular behavior is presumably sensitive to the value chosen for the requested precision $\eta = 10^{-8}$.

When $\chi < 1$, the schemes DS$_{C,W,R}$ and AS$_{C,W,R}$, that use $\mathbb{G}^{C,W,R}$, converge after a number of iterations not exceeding 430. As shown in Fig. (8), the number of iterations is nearly constant in the range $0 \le \chi \le 10^{-5}$. As expected, the accelerated schemes AS$_{C,W,R}$ are faster than the direct schemes DS$_{C,W,R}$, with scheme
<figure><img src="image_8.png" /><figcaption>Figure 7: Optimal reference Young modulus $E^0$ as a function of the contrast of properties $\chi$, 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 ($\chi = 0$) are indicated at the left of the graph. The material is a Boolean model of spheres with volume fraction 17%.</figcaption></figure>

$AS_R$ proving the fastest. With this scheme and when $\chi < 1$, the number of iterations is at most 168. Again, these results are qualitatively similar with that given in [17] in the context of conductivity.

For rigidly-reinforced media ($\chi > 1$), the number of iterations of schemes $DS_{C,W,R}$ and $AS_{C,W,R}$ follow the same powerlaw behaviors, with respect to $\chi$, 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 $\chi > 1$ than when $\chi < 1$. In this respect, benefits are to be expected from the use of dual Green operators [17], associated to the problem for the stress fields.

## 6.2. 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).
<figure><img src="image_9.png" /><figcaption>Figure 8: Number of iterations as a function of the contrast of properties $\chi$, for the various FFT methods, in log-log scale. Direct schemes: black and red; accelerated schemes: blue and green. Results for a porous material ($\chi = 0$) are indicated at the left of the graph. The material is a Boolean model of spheres with volume fraction 17%.</figcaption></figure>

but discretized on higher resolution grids of 256³ and 512³ voxels. The volume fraction of the spheres are respectively 16.82% and 16.85%. For simplicity, the contrast of properties take on two values $\chi = 10^{-4}$ and $\chi = 10^4$, so that the spheres are quasi-porous or quasi-rigid.

We perform iterations of the schemes AS and AS<sub>C,W,R</sub> using the optimized reference moduli found in the previous section, on the lower resolution grid. We apply the macroscopic strain loading:

$$ \bar{\varepsilon}_{ij} = \delta_{i1}\delta_{j1}. $$

At each iteration and for each scheme, we compute the elastic modulus $\tilde{C}_{11,11}$, derived from the mean $\langle\sigma_{11}\rangle$ of the stress component $\sigma_{11}$. The convergence rate toward the elastic modulus is represented in Fig. (9) for quasi-porous spheres with $L = 512$ and in Fig. (10) for quasi-rigid spheres with $L = 256$. In Fig. (9), 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<sub>R</sub> which are all represented. Dotted lines are guide to the eyes. In the porous case, much better convergence is obtained with scheme AS<sub>R</sub> than with schemes AS and AS<sub>C,W</sub>, as shown in Fig. (9). The estimate for $\tilde{C}_{11,11}$ predicted by AS and AS<sub>C,W</sub> present
<figure><img src="image_10.png" /><figcaption>Figure 9: Estimate of the elastic modulus $\tilde{C}_{11,11}$ as a function of the number of iterations performed, for a 3D Boolean model of quasi-porous spheres. Black symbols: accelerated schemes AS, $AS_C$, $AS_W$; red: scheme $AS_R$ (orange: Hashin and Shtrikman upper bound).</figcaption></figure>

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 $\tilde{C}_{11,11}$ with a precision of $10^{-2}$. The slower convergence rate follows that observed in Sec. (6), where convergence is much poorer for $\chi > 1$ than for $\chi < 1$, and where $AS_R$ is less advantageous compared to the other schemes. Nevertheless, in this case also, as shown in Fig. (10), smaller oscillations are observed in the estimate for $\tilde{C}_{11,11}$ when using $AS_R$ rather than schemes AS or $AS_{C,W}$.

## 7. 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 [17]. Compared to schemes using trigonometric polynomials as approximation space, or to other finite-differences
<figure><img src="image_11.png" /><figcaption>Figure 10: Estimate of the elastic modulus $\tilde{C}_{11,11}$ as a function of the number iterations performed, for a 3D Boolean model of quasi-rigid spheres. Black lines: accelerated schemes AS, AS<sub>C</sub>, AS<sub>W</sub>; red: scheme AS<sub>R</sub>.</figcaption></figure>

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 [31]. The new method also provides better estimates for the effective elastic moduli. Furthermore, its estimates do 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” [32] or with FFT algorithms that are less sensitive to the reference [14, 15], leading to the same local fields.

## 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.
[2] J. Li, S. Meng, X. Tian, F. Song, and C. Jiang. A non-local fracture model for composite laminates and numerical simulations by using the fft method. *Composites Part B: Engineering*, 43(3):961–971, 2011.

[3] M. Faessel and D. Jeulin. Segmentation of 3D microtomographic images of granular materials with the stochastic watershed. *Journal of Microscopy*, 239(1):17–31, 2010.

[4] F. Willot, L. Gillibert, and D. Jeulin. Microstructure-induced hotspots in the thermal and elastic responses of granular media. *International Journal of Solids and Structures*, 50(10):1699–1709, 2013.

[5] C. F. Dunant, B. Bary, A. B. Giorla, C. Pniguel, J. Sanahuja, C. Toulemonde, A. B. Tran, F. Willot, and J. Yvonnet. A critical comparison of several numerical methods for computing effective properties of highly heterogeneous materials. *Advances in Engineering Software*, 58:1–12, 2013.

[6] F. Willot and D. Jeulin. Elastic and electrical behavior of some random multiscale highly-contrasted composites. *International Journal for Multiscale Computational Engineering: special issue on Multiscale modeling and uncertainty quantification of heterogeneous materials*, 9(3):305–326, 2011.

[7] J. Escoda, F. Willot, D. Jeulin, J. Sanahuja, and C. Toulemonde. Estimation of local stresses and elastic properties of a mortar sample by FFT computation of fields on a 3D image. *Cement and Concrete Research*, 41(5):542–556, 2011.

[8] Morphhom software, http://cmm.ensmp.fr/morphhom, accessed October 4, 2014.

[9] CraFT software, http://craft.lma.cnrs-mrs.fr, accessed October 4, 2014.

[10] GeoDict software, http://www.geodict.de, accessed October 4, 2014.

[11] H. Moulinec and P. Suquet. A fast numerical method for computing the linear and non linear mechanical properties of the composites. *Comptes rendus de l'Academie des sciences, Série II*, 318(11):1417–1423, 1994.

[12] J. Vondřejc, J. Zeman, and I. Marek. An FFT-based Galerkin method for homogenization of periodic media. *Computers & Mathematics with Applications*, 68(3):156–173, 2014.
[13] S. Brisard and L. Dormieux. Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites. *Computational Methods for Applied Mechanical Engineering*, 217(220):197-212, 2012.

[14] J. Zeman, J. Vodrejc, J. Novak, and I. Marek. Accelerating a FFT-based solver for numerical homogenization of a periodic media by conjugate gradients. *Journal of Computational Physics*, 229(21):8065-8071, 2010.

[15] S. Brisard and L. Dormieux. FFT-based methods for the mechanics of composites: A general variational framework. *Computational Materials Science*, 49(3):663-671, 2010.

[16] S. Brisard and F. Legoll. Periodic homogenization using the lippmann-schwinger formalism, 2014. http://arxiv.org/abs/1411.0330.

[17] F. Willot, B. Abdallah, and Y.-P. Pellegrini. Fourier-based schemes with modified green operator for computing the electrical response of heterogeneous media with accurate local fields. *International Journal for Numerical Methods in Engineering*, 98(7):518-533, 2014.

[18] W. H. Müller. Mathematical vs. experimental stress analysis of inhomogeneities in solids. *Journal de Physique*, 6(C1):139-148, 1996.

[19] C.M. Brown, W. Dreyer, and W.H. Müller. Discrete fourier transforms and their application to stress-strain problems in composite mechanics: a convergence study. *Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences*, 458:1967-1987, 2002.

[20] F. Willot and Y.-P. Pellegrini. Fast Fourier transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwise-disordered porous media. In *D. Jeulin, S. Forest (eds)*, "Continuum Models and Discrete Systems CMDS 11", pages 443-449, Paris, 2008. École des Mines.

[21] M. F. Thorpe and P. N. Sen. Elastic moduli of two-dimensional composite continua with elliptic inclusions. *Journal of the Acoustical Society of America*, 77(5):1674-1680, 1985.

[22] G. W. Milton. *The Theory of Composites*. Cambridge Univ. Press, Cambridge, 2002.
[23] S. Kanaun and V. Levin. *Self-consistent methods for composites*. Springer, Dordrecht, The Netherlands, 2008.

[24] H. Moulinec and P. Suquet. A numerical method for computing the overall response of nonlinear composites with complex microstructure. *Computer Methods in Applied Mechanics and Engineering*, 157(1):69–94, 1998.

[25] D. J. Eyre and G. W Milton. A fast numerical scheme for computing the response of composites using grid refinement. *The European Physical Journal Applied Physics*, 6(1):41–47, 1999.

[26] J.-C. Michel, H. Moulinec, and P. Suquet. A computational scheme for linear and non-linear composites with arbitrary phase contrast. *International Journal for Numerical Methods in Engineering*, 52(1-2):139–160, 2001.

[27] H. Moulinec and F. Silva. Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials. *International Journal for Numerical Methods in Engineering*, 97(13):960–985, 2014.

[28] I. Cohen. Simple algebraic approximations for the effective elastic moduli of cubic arrays of spheres. *Journal of the Mechanics and physics of Solids*, 52(9):2167–2183, 2004.

[29] F. Willot and D. Jeulin. Elastic behavior of composites containing boolean random sets of inhomogeneities. *International Journal of Engineering Science*, 47(2):313–324, 2009.

[30] H. Moulinec and P. Suquet. Comparison of FFT-based methods for computing the response of composites with highly contrasted mechanical properties. *Physica B: Condensed Matter*, 338(1-4):58–60, 2003.

[31] N. Lahellec, J.-C. Michel, H. Moulinec, and P. Suquet. Analysis of inhomogeneous materials at large strains using fast fourier transforms. In *Proc. IUTAM Symposium on Computational Mechanics of Solid Materials at Large Strains*, pages 247–258. Kluwer Academic Publishers, 2001.

[32] J.-C. Michel, H. Moulinec, and P. Suquet. A computational method based on augmented Lagrangians and fast fourier transforms for composites with high contrast. *Computer Modelling in Engineering & Sciences*, 1(2):79–88, 2000.