# Comptes Rendus
## Mécanique

Rémi Cornaggia, Marie Touboul and Cédric Bellis

FFT-based computation of homogenized interface parameters

Volume 350 (2022), p. 297-307

Published online: 4 July 2022

https://doi.org/10.5802/crmeca.119

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

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

Les Comptes Rendus. Mécanique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org

e-ISSN : 1873-7234
Short paper / *Note*

# FFT-based computation of homogenized interface parameters

Rémi Cornaggia<sup>a, *, a</sup>, Marie Touboul<sup>b</sup> and Cédric Bellis<sup>c</sup>

<sup>a</sup> Sorbonne Université, CNRS, UMR 7190, Institut Jean Le Rond d'Alembert, F-75005 Paris, France

<sup>b</sup> School of Mathematics, University of Manchester, Oxford Road, Manchester, M13 9PL, UK

<sup>c</sup> Aix Marseille Univ, CNRS, Centrale Marseille, LMA UMR 7031, Marseille, France

*E-mails*: remi.cornaggia@sorbonne-universite.fr (R. Cornaggia), marie.touboul@manchester.ac.uk (M. Touboul), bellis@lma.cnrs-mrs.fr (C. Bellis)

**Abstract.** The homogenization of microstructured interfaces requires solving specific problems posed on semi-infinite bands. To tackle these problems with existing FFT-based algorithms, a reformulation of these band problems into fully periodic cell problems, posed on bounded domains, is established. This is performed thanks to a Dirichlet-to-Neumann operator and a decomposition of the solution involving a boundary corrector, in a Fourier framework. A fixed-point algorithm and an example choice of corrector are proposed. Comparisons with other computational methods support this proposition.

**Keywords.** Homogenization, Dirichlet-to-Neumann, Cell problems, Band problems, FFT-based solvers.

Manuscript received 14 March 2022, revised 31 May 2022, accepted 7 June 2022.

## 1. Introduction

Homogenization of microstructured interfaces that are periodic along one direction, as proposed for antiplane shear waves by [1], results in non-standard transmission conditions across effective interfaces, see also [2]. These conditions are governed by a set of effective coefficients, that are the counterpart of well-known effective properties of fully-periodic materials. In the latter case, effective properties are obtained by solving so-called *cell problems* on representative cells with periodicity conditions on all boundaries, see e.g. [3] and the references therein. For microstructured interfaces, the effective coefficients are determined by solving what we call here a *band problem*, i.e. a problem posed on a semi-infinite band with periodicity conditions on its boundaries. Accurately solving this band problem is therefore crucial for the effective model to approximate the microstructured configuration, and numerical efficiency is required when multiple computations are performed, e.g. for topological optimization [2, 3].

<footnote>* Corresponding author.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Microstructured interface and homogenization principle and (b) representative band $Y_\infty$ and its partitioning.</figcaption></figure>

The main difficulty, compared to fully-periodic cell problems, consists in dealing with unbounded domains. Such a problem is encountered in the literature in waveguides, and transparent boundary conditions are proposed to restrain the computational domain to a bounded subset of the guide [4, 5]. The “restrained” problem is found to be very similar to a cell problem, except for the non-classical boundary conditions associated with the restriction.

This note shows how these tools can be applied to band problems, using classical cell problems solvers, namely FFT-based iterative solvers first proposed by [6] and now widely used and improved [7]. Extending FFT-based solvers to non-periodic boundary conditions has been proposed earlier, e.g. by [8] to account for Dirichlet conditions.

In Section 2, the band problem is first restrained to a bounded cell thanks to a Dirichlet-to-Neumann (DtN) operator, following a classical methodology [4]. Then, the solution of this problem is decomposed into a bi-periodic part (recovering the setting of cell problems) and a boundary corrector accounting for the DtN conditions. Using Fourier series expansion, we then build an explicit expression of the boundary corrector. An iterative FFT-based algorithm is then proposed as an example of implementation in Section 3, and numerical comparisons with two other methods are provided. Some perspectives conclude the paper.

## 2. Band problem and reformulation as a cell problem

In the two-dimensional plane, endowed with the classical Cartesian basis ($e_1, e_2$), we consider the interaction of scalar waves (e.g. antiplane shear waves) with a periodically microstructured interface (e.g. a row of inclusions, see Figure 1(a)). This interface is parallel to the $e_2$-direction, centered at $x_1 = 0$, occupies the domain $[-b/2, b/2] \times \mathbb{R}$ and is 1-periodic in $x_2$, and supports the variations of the normalized shear modulus $\mu$, relatively to $\mu = 1$ in the background domain $\{|x_1| > b/2\}$. The microstructure is therefore entirely characterized by the representative band $Y_\infty = \mathbb{R} \times S$ with $S = [-1/2, 1/2]$, see Figure 1(b).

In [1] it is shown that for long wavelengths, the microstructure can be replaced by non-trivial jump conditions on velocity and normal stress across an effective interface, see Figure 1(a). These conditions are written in terms of vector-valued coefficients $B = (B_1, B_2)$ and $C = (C_1, C_2)$ that are defined as:

$$ B = \lim_{x_1 \to +\infty} W - \lim_{x_1 \to -\infty} W \quad \text{and} \quad C = \int_{Y_\infty} \mu(\mathbf{x}) \partial_2 W(\mathbf{x}) \, d\mathbf{x}, \qquad (1) $$

in terms of the vector-valued field $W = (W_1, W_2)$, solution of the following band problem:

$$ \left\{ \begin{array}{l} \nabla \cdot (\mu[I + \nabla W]) = 0 \quad \text{in } Y_\infty, \\ W \text{ is 1-periodic in the } x_2 \text{ variable,} \\ \lim_{x_1 \to \pm\infty} \nabla W = 0, \end{array} \right. \qquad (2) $$

where $I$ is the second-order identity tensor.
**Remark 1.** It can be seen that the two components ($W_1, W_2$) solve uncoupled problems. Moreover, it is shown in [1] that $B_2 = -C_1$.

## 2.1. Modal decomposition and Dirichlet-to-Neumann operator

In this section, the band problem is reduced to a cell problem thanks to a modal decomposition in the homogeneous left and right halves of the band, and the introduction of the corresponding Dirichlet-to-Neumann (DtN) operators [4,5]. The band is first decomposed into two parts $Y^+$ and an interior cell $Y_b$ as:

$$Y_{\infty} = Y^{-} \cup Y_{b} \cup Y^{+},$$

with $Y^{-} = ]-\infty, -b/2] \times S$, $Y_b = [-b/2, b/2] \times S$ and $Y^{+} = [b/2, +\infty] \times S$, see Figure 1(b). The variations of $\mu$ are supported in the interior cell $Y_b$, so that the restriction of the solution to the right half-band $Y^+$, $W^+ := W|_{Y^+}$, is harmonic and satisfies

$$\begin{cases} \Delta W^+ = 0 & \text{in } Y^+, \\ W^+ \text{ is 1-periodic in the } x_2 \text{ variable,} \\ \lim_{x_1 \to +\infty} \nabla W^+ = 0. \end{cases}$$

Looking for $W^+$ as a separated variables function, and imposing the periodicity and condition at infinity, one obtains a classical decomposition onto a modal basis, see e.g. [5] or [1, Appendix B]:

$$W^{+}(x_1, x_2) = \sum_{n \in \mathbb{Z}} w_{n}^{+} e^{-|k_{n}|(x_{1}-b/2)} \phi_{n}(x_{2}), \quad \text{with } \phi_{n}(s) = e^{ik_{n}s} \text{ and } k_{n} = 2n\pi. \quad (3)$$

The Fourier modes $\phi_n$ are orthonormal with respect to the $L^2$ scalar product $(\cdot, \cdot)_S$ on a vertical section:

$$(\phi_p, \phi_q)_S = \delta_{pq} \quad \text{with} \quad (f, g)_S := \int_S f(x_2) \overline{g(x_2)} \, dx_2.$$

Looking at the section at $x_1 = b/2$, the modal coefficients $w_n^+$ can be expressed as:

$$w_n^+ = \left( W^+ \left( \frac{b}{2}, \cdot \right), \phi_n \right)_S = \int_S W^+ \left( \frac{b}{2}, x_2 \right) \overline{\phi_n(x_2)} \, dx_2. \quad (4)$$

Differentiating the modal decomposition (3) with respect to $x_1$ and using the expression (4) of the modal coefficients, one finally obtains the DtN operator $\Lambda^+$ linking the trace of $W$ and the normal derivative $\partial_1 W$ on the section $S$ at $x_1 = b/2$:

$$\partial_1 W \left( \frac{b}{2}, x_2 \right) = \Lambda^+ \left[ W \left( \frac{b}{2}, \cdot \right) \right] (x_2), \quad \text{where } \Lambda^+[f](x_2) := - \sum_{n \in \mathbb{Z}} (f, \phi_n)_S |k_n| \phi_n(x_2), \quad \forall x_2 \in S.$$

Similarly, the restriction $W^- = W|_{Y^-}$ to the left half-band $Y^-$ can be decomposed as:

$$W^{-}(x_1, x_2) = \sum_{n \in \mathbb{Z}} w_{n}^{-} e^{ik_{n}|(x_{1}+b/2)} \phi_{n}(x_{2}). \quad (5)$$

The left DtN relation (using again the outer normal derivative i.e. $-\partial_1 W$) is then found similarly for all $x_2 \in S$ as:

$$-\partial_1 W \left( -\frac{b}{2}, x_2 \right) = \Lambda^{-} \left[ W \left( -\frac{b}{2}, \cdot \right) \right] (x_2), \quad \text{where } \Lambda^{-}[f](x_2) := - \sum_{n \in \mathbb{Z}} (f, \phi_n)_S |k_n| \phi_n(x_2).$$

The DtN operator is found to be the same at both interfaces $x_1 = \pm b/2$, therefore we simply note:

$$\partial_n W(\pm b/2, \cdot) = \Lambda[W(\pm b/2, \cdot)], \quad \text{where } \Lambda = \Lambda^+ = \Lambda^-.$$

The band problem (2) can finally be rewritten as a *cell problem* on the bounded cell $Y_b$, with mixed periodic and DtN boundary conditions:

$$\begin{cases} \nabla \cdot (\mu[I + \nabla W]) = 0 & \text{in } Y_b, \\ W \text{ is 1-periodic in the } x_2 \text{ variable,} \\ \partial_n W(\pm b/2, \cdot) = \Lambda[W(\pm b/2, \cdot)]. \end{cases} \quad (6)$$
**Remark 2.** Introducing the above decompositions (3) and (5) into the definitions (1) of the effective coefficients also provides alternative expressions for these coefficients:

$$ B = w_0^+ - w_0^- = \int_S [W(b/2, x_2) - W(-b/2, x_2)] dx_2 \quad \text{and} \quad C = \int_{Y_b} \mu(x) \partial_2 W(x) dx, \quad (7) $$

that can be computed from the solution of problem (6).

## 2.2. Introducing a bi-periodic function as the new unknown

As stated in the introduction, we aim at using existing solvers designed to address cell problems with periodic boundary conditions. To account for the (non-periodic) DtN conditions in the horizontal direction, we introduce the following decomposition:

$$ W = W_{\#} + W_{\mathrm{B}}, \quad (8) $$

where

* $W_{\#}$ is a bi-periodic function, that admits a Fourier series expansion:

$$ W_{\#}(x_1, x_2) = \sum_{m \in \mathbb{Z}} \sum_{n \in \mathbb{Z}} w_{mn} \phi_m^b(x_1) \phi_n(x_2), \quad \text{with } \phi_m^b(x_1) := \phi_m(x_1/b) = e^{ik_m x_1/b}. \quad (9) $$

* $W_{\mathrm{B}}$ is a corrector introduced so that the sum $W_{\#} + W_{\mathrm{B}}$ satisfies the DtN boundary conditions (6). This corrector being periodic in the $x_2$ variable, it is natural to look for it as a Fourier series expansion in this variable:

$$ W_{\mathrm{B}}(x_1, x_2) = \sum_{n \in \mathbb{Z}} \psi_n(x_1) \phi_n(x_2), \quad (10) $$

in terms of Fourier modes $\phi_n$ and yet-to-be-chosen functions $\psi_n$.

Writing the boundary conditions $(\partial_n - \Lambda)W_{\mathrm{B}} = -(\partial_n - \Lambda)W_{\#}$ at $x_1 = \pm b/2$, recalling that $\Lambda\phi_n = -|k_n|\phi_n$ and projecting each equation onto the modes $\phi_n$, one shows that the functions $\psi_n$ satisfy two uncoupled Robin boundary conditions:

$$ x_1 = -\frac{b}{2}: \quad (-\partial_1 \psi_n + |k_n| \psi_n)(-b/2) = s_n^{(2)}/b - |k_n| s_n^{(1)}, \quad (11) $$

$$ x_1 = \frac{b}{2}: \quad (\partial_1 \psi_n + |k_n| \psi_n)(b/2) = -s_n^{(2)}/b - |k_n| s_n^{(1)}, \quad (12) $$

where the terms $s_n^{(1)}$ and $s_n^{(2)}$ are weighted sums of the Fourier coefficients of $W_{\#}$:

$$ s_n^{(1)} = \sum_{m \in \mathbb{Z}} (-1)^m w_{mn} \quad \text{and} \quad s_n^{(2)} = i \sum_{m \in \mathbb{Z}} (-1)^m k_m w_{mn}. \quad (12) $$

Note that if $W_{\#}$ is a real-valued function, then $w_{mn} = \overline{w(-m)n}$, so that $s_n^{(1)}$ and $s_n^{(2)}$ are also real.

As two conditions must be satisfied by each function $\psi_n$, one may look for linear combinations of two scalar-valued functions $(\psi_n^-, \psi_n^+)$, i.e. $\psi_n = z_n^+ \psi_n^+ + z_n^- \psi_n^-$, with $z_n^\pm$ being constant vector-valued coefficients. These functions must be chosen to ensure that the conditions (11) turn into an invertible linear system satisfied by the coefficients $z_n^\pm$ for each $n$. An example is given in Section 3.2 below. As a consequence, the corrector $W_{\mathrm{B}}$ is explicitly determined by the Fourier coefficients of the periodic function $W_{\#}$, i.e. there is a linear mapping $\mathcal{B}$ such that:

$$ W_{\mathrm{B}} = \mathcal{B} W_{\#}. \quad (13) $$

Finally, injecting the decomposition (8) and the mapping (13) into the original problem (6), the new unknown function $W_{\#}$ satisfies the PDE:

$$ \nabla \cdot (\mu[I + \nabla(W_{\#} + W_{\mathrm{B}})]) = 0, \quad \text{i.e. } \nabla \cdot (\mu[I + \nabla(I + \mathcal{B})W_{\#}]) = 0 \quad \text{in } Y_b, \quad (14) $$

along with bi-periodic boundary conditions, i.e. the initial band problem (2) has been reformulated into a (fully-periodic) cell problem on the periodic part $W_{\#}$ of the solution.
### 3. Example of implementation

Solving the new cell problem requires a numerical scheme to address the non-classical PDE (14), in which the “boundary correction” operator $\mathcal{B}$ has to be chosen as discussed above (i.e. the functions $\psi_n^{\pm}$ must be chosen). The choice of this operator may strongly affect the properties of a given solver (such as stability, convergence and efficiency).

In this part, an example of solver-corrector couple is proposed. A complete analysis of its performance is deferred to future work, but numerical illustrations are provided to assess the interest of the proposed choice.

#### 3.1. Fixed-point algorithm

Following [6], a reference homogeneous medium with constant material property $\mu_0$ is introduced, as well as the material fluctuation $\delta\mu = \mu - \mu_0$. Then (14) is rewritten in terms of a polarization field $\boldsymbol{\tau}$:

$$ \nabla \cdot (\mu_0 \nabla \boldsymbol{W}_{\#} + \boldsymbol{\tau}) = \mathbf{0}, \quad \text{with:} \quad \boldsymbol{\tau} = \delta\mu \nabla \boldsymbol{W}_{\#} + \mu(I + \nabla \mathcal{B} \boldsymbol{W}_{\#}). \quad (15) $$

When the polarization $\boldsymbol{\tau}$ is known, the first of these equations admits the following solution, equivalently expressed in space and in the Fourier domain:

$$ \boldsymbol{W}_{\#}(\boldsymbol{x}) = (\boldsymbol{\gamma}_0 \boldsymbol{\tau})(\boldsymbol{x}) \Leftrightarrow \hat{\boldsymbol{W}}_{\#}(\boldsymbol{\xi}) = \hat{\boldsymbol{\gamma}}_0(\boldsymbol{\xi}) \cdot \hat{\boldsymbol{\tau}}(\boldsymbol{\xi}), \quad (16) $$

where $\hat{f}$ is the Fourier transform of a function $f$, $\boldsymbol{\xi}$ is the Fourier variable, and $\boldsymbol{\gamma}_0$ is the Green operator for the displacement field in a bi-periodic cell associated with the reference medium $\mu_0$. This is a convolution operator in space whose Fourier transform is written:

$$ \hat{\boldsymbol{\gamma}}_0(\boldsymbol{\xi}) = i\boldsymbol{\xi}/(\mu_0|\boldsymbol{\xi}|^2) \quad \text{for } \boldsymbol{\xi} \neq \mathbf{0} \text{ and } \hat{\boldsymbol{\gamma}}_0(\mathbf{0}) = \mathbf{0}. $$

The value $\hat{\boldsymbol{\gamma}}_0(\mathbf{0})$ is set here to impose a null mean to the field $\boldsymbol{W}_{\#}$: this is just a convenient choice one can make since $\boldsymbol{W}$ and therefore $\boldsymbol{W}_{\#}$ are defined up to a constant irrelevant to the computation of the effective properties ($\mathbf{B}$, $\mathbf{C}$), see (7).

Then, using the property $\boldsymbol{\gamma}_0(\mu_0 \nabla \boldsymbol{u}) = -\boldsymbol{u}$, we rewrite formally (15) and (16) as:

$$ (i) \quad \boldsymbol{W}_{\#} = \boldsymbol{W}_{\#} + \boldsymbol{\gamma}_0 \mathbf{S}[\boldsymbol{W}_{\#}], \quad (ii) \quad \mathbf{S}[\boldsymbol{W}_{\#}] = \mu(I + \nabla(I + \mathcal{B})\boldsymbol{W}_{\#}), $$

where (i) stands for the equilibrium of the antiplane stress $\mathbf{S}$ and (ii) for the constitutive relation that defines this stress. Finally, the original problem is solved by applying standard fixed-point iterations to these equations, and alternating back and forth between space and Fourier domains to avoid any convolution product. The final algorithm is as follows:

$$ \text{Initialisation: } \boldsymbol{W}_{\#}^{(0)} = \mathbf{0}, \quad \mathbf{S}^{(0)} = \mu I, $$

Iterate $i+1$: $\hat{\boldsymbol{W}}_{\#}^{(i)}$ and $\hat{\boldsymbol{S}}^{(i)}$ are known

$$ \begin{align} (a) \quad \hat{\boldsymbol{W}}_{\#}^{(i+1)} &= \hat{\boldsymbol{W}}_{\#}^{(i)} + \hat{\boldsymbol{\gamma}}_0 \cdot \hat{\boldsymbol{S}}^{(i)} \tag{17} \\ (b) \quad \nabla \boldsymbol{W}_{\#}^{(i+1)} &= \mathcal{F}^{-1}(i\boldsymbol{\xi} \otimes \hat{\boldsymbol{W}}_{\#}^{(i+1)}) \notag \end{align} $$

$$ (c) \quad \nabla \mathcal{B} \boldsymbol{W}_{\#}^{(i+1)} \text{ is computed analytically from } \hat{\boldsymbol{W}}_{\#}^{(i+1)} $$

$$ (d) \quad \mathbf{S}^{(i+1)} = \mu(I + \nabla \boldsymbol{W}_{\#}^{(i+1)} + \nabla \mathcal{B} \boldsymbol{W}_{\#}^{(i+1)}) $$

$$ (e) \quad \hat{\boldsymbol{S}}^{(i+1)} = \mathcal{F}(\mathbf{S}^{(i+1)}) $$

$$ (f) \quad \text{Stopping criterion: exit if } \|\boldsymbol{W}^{(i+1)} - \boldsymbol{W}^{(i)}\|_{L^2(Y_b)} / \|\boldsymbol{W}^{(i)}\|_{L^2(Y_b)} < \varepsilon_{\text{FP}} $$

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the Fourier transform and its inverse, and $\varepsilon_{\text{FP}}$ is a user-defined tolerance. Note that the chosen criterion (f) assumes that the algorithm converges, which is not established yet, as discussed now.
**Remark 3.** The result of the iterations can be written as the Neumann series

$$W_{\#}^{(j)} = \sum_{i<j} (\mathcal{N}_0)^i (\gamma_0 \mu I), \quad \text{with } \mathcal{N}_0 = [I + \gamma_0(\mu \nabla(I + \mathcal{B}))]$$

which is convergent only if $\rho(\mathcal{N}_0) < 1$, where $\rho$ denotes the spectral radius. For “classical” cell problems, i.e. when $\mathcal{B} = 0$, this condition is enforced by choosing adequately the reference medium. In particular, choosing $\mu_0 = (\min \mu + \max \mu)/2$ ensures convergence, and even sometimes optimal convergence rate, see [9] and the references therein. Extending this analysis to the present case with $\mathcal{B} \neq 0$ is crucial to the complete validation of the proposed methodology, but is outside the scope of this note and will be considered in future work.

## 3.2. An example of corrector

We now propose a specific choice for the corrector $W_B$ given by (10). By choosing $\psi_0$ as a linear function, written explicitly in terms of the coefficient $B$ according to the definition (7), and the other functions $\psi_n^{\pm}$ for $n \neq 0$ as exponentials as in expressions (3) and (5), one defines:

$$W_B(x_1, x_2) = B \frac{x_1}{b} + \sum_{n \neq 0} \left( z_n^- e^{|k_n|(x_1+b/2)} + z_n^+ e^{-|k_n|(x_1-b/2)} \right) \phi_n(x_2). \quad (18)$$

Injecting this form into the Robin boundary conditions (11), one finds:

$$B = -s_0^{(2)}, \quad z_n^- = \frac{e^{-|k_n|b}}{2} \left( -\frac{1}{b|k_n|} s_n^{(2)} - s_n^{(1)} \right), \quad \text{and} \quad z_n^+ = \frac{e^{-|k_n|b}}{2} \left( \frac{1}{b|k_n|} s_n^{(2)} - s_n^{(1)} \right), \quad (19)$$

with the sums $s_n^{(1)}$ and $s_n^{(2)}$ given by (12).

**Remark 4.** Linear functions and exponential functions for all functions $\psi_n^{\pm}$ were also considered to build the corrector. The results are omitted here for brevity, but the choice above was clearly superior in terms of convergence speed of the fixed-point scheme, in the configuration of numerical tests presented below.

## 3.3. Numerical results

As proofs of concept for the proposed algorithm and choice of corrector, we compare the solutions obtained with the present approach to the ones given by two other methods, namely (i) the mode-matching method described in [1, Appendix B] and (ii) the finite element method on truncated bands. Two distinct microstructures are studied. Qualitative comparisons are provided by displaying chosen solutions, while quantitative comparisons are made on the evaluation of the homogenized coefficients for various material contrasts.

### 3.3.1. Discretisation of the proposed method

The various fields that intervene in algorithm (17) are discretized on a regular grid mapping the inner cell $Y_b$ and the Fourier transforms are performed using standard FFT libraries. In particular, the Fourier transform $\hat{W}_{\#}$ is discretized onto the set of coefficients $w_{mn}$ introduced in (9) and that are used to compute $W_B = \mathcal{B}W_{\#}$ in step (c) of (17). With the proposed corrector (18), the coefficient $B$ is explicitly given by (19), while the integral that defines $C$ in (7) is computed numerically as a discrete version of (7).

### 3.3.2. Rectangular inclusion: comparison with the mode-matching method

When the microstructured interface is made of rectangular inclusions aligned with the frame $(e_1, e_2)$, the mode-matching (MM) method can be used to solve the original band problem (2) and compute the homogenized coefficients as described in [1, Appendix B].
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Periodic part, (b) corrector and (c) total fields ($W_1, W_2$) on a 129² pixel grid, and (d) reference solution obtained with the mode matching method (MM) with 65 modes, for a stiff rectangular inclusion of size 0.8 × 0.5 and material property $\mu_i = 10$.</figcaption></figure>

<table><caption>**Table 1.** Iterations until convergence of the fixed-point scheme, for the rectangular inclusion and a tolerance $\epsilon_{FP} = 10^{-4}$ for the stopping criterion. The reference medium is chosen as $\mu_0 = (\min\mu + \max\mu)/2$, except for $\mu_i = 0.001$ (indicated by $(\star)$), for which fixing $\mu_0 = \max\mu = 1$ allowed convergence</caption><tr><td>Inclusion stiffness $\mu_i$</td><td>0.001</td><td>0.01</td><td>0.1</td><td>10</td><td>100</td><td>1000</td></tr><tr><td>Iterations</td><td>24 (★)</td><td>12</td><td>9</td><td>20</td><td>14</td><td>13</td></tr></table>

To tackle this configuration with the proposed framework, we consider a rectangular inclusion of lengths $(\ell_1, \ell_2) = (0.8, 0.5)$ embedded in a square interior cell $Y_b = [-0.5, 0.5]^2$. Both the inclusion and surrounding medium are homogeneous, with properties $\mu_i = 10$ and $\mu_m = 1$, respectively.

Figure 2(a–c) display the fields obtained with the proposed iterative algorithm on a square 129² pixels grid. We observe that (i) for $W_1$, the corrector is dominated by its linear part as seen from the amplitude of the fields and on their cross-sections displayed in Figure 3, and (ii) for $W_2$, there is no linear part since $B_2 = 0$ for such a symmetric configuration, see [1], and the corrector amplitude is much lower than this of the “periodic” part.

Moreover, the solution obtained with the MM method, using 65 modes and restricted to $Y_b$, is also plotted in Figure 2(d). An excellent qualitative agreement is observed. Figure 4 then provides a more quantitative comparison: for different material contrasts ranging from $10^{-3}$ to $10^3$, the homogenized coefficients are computed with both methods and compared. The relative difference stays below 5% in all cases. More surprisingly, the proposed iterative scheme converges in very few iterations, see Table 1, while the fixed-point scheme initially proposed by [6] and adapted to the present context is known to converge slowly for large material contrasts, see [9]. Without further convergence analysis, see Remark 3, we can only speculate that this might be a geometrical effect.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Cross-sections of the field $W_1$ displayed in Figure 2(a-c) for (a) $x_2 = 0$ (i.e. through the rectangular inclusion) and (b) $x_2 = -1/2$ (through the matrix).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Relative errors made on the homogenized coefficients $B_1$ and $C_2$ corresponding to a rectangular inclusion and for various material contrasts, with the mode matching method taken as a reference.</figcaption></figure>

### 3.3.3. Elliptic inclusion: comparison with the finite element method

For a second example, and to show the present method interest for non-rectangular geometries (where the mode matching method is not applicable), we consider an elliptic inclusion of semi-axes $(a_1, a_2) = (0.5, 0.2)$, whose major axis is tilted with an angle $\theta = 40^\circ$ w.r.t. $e_1$, and with material property $\mu_i = 10$. Figure 5(a-c) display the fields obtained with the proposed iterative algorithm on a square $129^2$ pixels grid. Since the cell is not symmetric with respect to the vertical axis $x_1 = 0$, there is one additional non-zero effective coefficient $B_2 = -C_1 \neq 0$, and the correctors (dominated by their linear part associated with **B**) are of the same magnitude than the periodic part for both components ($W_1, W_2$).

A finite element solution is also displayed in Figure 5(d), obtained thanks to the platform FreeFEM++ [10] on a truncated band $Y_L = [-L/2, L/2] \times S$ (with $L = 10$ for the displayed examples). Homogeneous Neumann boundary conditions $\partial_1 W = 0$ are imposed at the boundaries $x_1 = \pm L/2$.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Periodic part, (b) corrector and (c) total fields ($W_1, W_2$) on a 129² pixel grid, and (d) solution obtained with the finite element method (FE), for a stiff elliptic inclusion of semi-axes ($a_1, a_2$) = (0.5, 0.2), tilted with an angle $\theta$ = 40° and material property $\mu_i = 10$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Relative errors made on the homogenized coefficients $B_1, B_2 = -C_1$ and $C_2$ corresponding to an elliptic inclusion and for various material contrasts, with the finite element method taken as a reference.</figcaption></figure>

to approach the decaying condition $\lim_{x_1 \to \pm\infty} \nabla W = 0$. We refer to [2, Section 4.3.1.4] for more detail on this approach.

Again, a very good qualitative agreement between both methods is observed in the moderate contrast case $\mu_i = 10$. For various material contrasts, Figure 6 shows a reasonable discrepancy between both methods, although the error on the coefficient $B_2 = -C_1$ increases for stiffer inclusions. In Table 2 are presented the iterations needed to reach convergence. This time the expected characteristic slow-down of the fixed-point scheme is clearly observed for high contrasts, either for stiffer or softer inclusions. In the latter case, the choice of the reference property $\mu_0 = (\max\mu + \min\mu)/2$ fails to ensure convergence, but the choice $\mu_0 = \max\mu$ works. These preliminary observations confirm the need for a detailed convergence analysis as discussed in Remark 3.
<table><caption>**Table 2.** Iterations until convergence of the fixed-point scheme, for the elliptic inclusion and a tolerance $\epsilon_{FP} = 10^{-4}$ for the stopping criterion. The reference medium is chosen as $\mu_0 = (\min\mu + \max\mu)/2$, except for $\mu_i = 0.001$ and $\mu_i = 0.01$ (indicated by $(\star)$), for which fixing $\mu_0 = \max\mu = 1$ allowed convergence</caption><tr><td>Inclusion stiffness $\mu_i$</td><td>0.001</td><td>0.01</td><td>0.1</td><td>10</td><td>100</td><td>1000</td></tr><tr><td>Iterations</td><td>183 (★)</td><td>64 (★)</td><td>11</td><td>20</td><td>66</td><td>266</td></tr></table>

**4. Conclusion**

The main purpose of this note was to present a reformulation of semi-infinite band problems into fully periodic cell problems. This is performed thanks to a Dirichlet-to-Neumann operator and a decomposition involving a boundary corrector, in a Fourier-based framework, having in mind a FFT-based implementation. As proofs of concept, a fixed-point algorithm and an example choice a corrector are given. The resulting procedure produces relevant results and is successfully employed for topological optimization of microstructured interfaces in [2, Chap. 4]. The established framework (modal decomposition, DtN operator, Fourier representation) should be easily extended to other material or geometrical configuration where similar band problems are of interest.

The key point of the convergence of the fixed-point algorithm is not settled yet. An analysis tackling this point is necessary, and should also pave the way on extensions to other algorithms, using the rich literature on the FFT-based algorithms, see e.g. [9] or the recent review [7].

Finally, other problems posed on domains with non-fully-periodic boundary conditions, arising in other homogenization procedures, could be considered. One may think e.g. to problems posed on a half-band to account for a rigid wall or void (Dirichlet or Neumann boundary condition) [11], band problems arising when looking for transmission conditions between microstructured materials [12, Chap. 4], or more classical Dirichlet or Neumann boundary conditions on bounded cells, where our approach based on a corrector function should be compared to alternative propositions [8].

**Conflicts of interest**

The authors declare no conflicts of interest.

**Dedication**

The manuscript was written through contributions from all authors. All authors have given approval to the final version of the manuscript.

**References**

[1] J.-J. Marigo, A. Maurel, K. Pham, A. Sbitti, "Effective dynamic properties of a row of elastic inclusions: The case of scalar shear waves", *J. Elast.* **128** (2017), no. 2, p. 265-289.

[2] M. Touboul, "Acoustic and elastic wave propagation in microstructured media with interfaces: homogenization, simulation and optimization", PhD Thesis, Aix-Marseille Université, 2021, https://tel.archives-ouvertes.fr/tel-03411353/.

[3] R. Cornaggia, C. Bellis, "Tuning effective dynamical properties of periodic media by FFT-accelerated topological optimization", *Int. J. Numer. Methods Eng.* **121** (2020), no. 14, p. 3178-3205.

[4] I. Harari, I. Patlashenko, D. Givoli, "Dirichlet-to-Neumann maps for unbounded wave guides", *J. Comput. Phys.* **143** (1998), no. 1, p. 200-223.
[5] A.-S. Bonnet-Ben Dhia, G. Legendre, "An alternative to Dirichlet-to-Neumann maps for waveguides", C. R. Mat. **349** (2011), no. 17-18, p. 1005-1009.

[6] H. Moulinec, P. Suquet, "A numerical method for computing the overall response of nonlinear composites with complex microstructure", Comput. Meth. Appl. Mech. Eng. **157** (1998), no. 1, p. 69-94.

[7] M. Schneider, "A review of nonlinear FFT-based computational homogenization methods", Acta Mech. **232** (2021), no. 6, p. 2051-2100.

[8] L. Gélébart, "A modified FFT-based solver for the mechanical simulation of heterogeneous materials with Dirichlet boundary conditions", C. R. Méc. **348** (2020), no. 8-9, p. 693-704.

[9] H. Moulinec, F. Silva, "Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials", Int. J. Numer. Methods Eng. **97** (2014), no. 13, p. 960-985.

[10] F. Hecht, "New development in FreeFem++", J. Numer. Math. **20** (2012), no. 3-4, p. 251-265.

[11] A. Maurel, J.-J. Marigo, K. Pham, "Effective boundary condition for the reflection of shear waves at the periodic rough boundary of an elastic body", *Vietnam J. Mech.* **40** (2018), no. 4, p. 303-323.

[12] S. Fliss, "Wave propagation in periodic media: mathematical analysis and numerical simulation", 2019, Habilitation à diriger des recherches, Université Paris Sud (Paris 11), https://hal.archives-ouvertes.fr/tel-02394976.