# An accurate and robust FFT-based solver for transient diffusion in heterogeneous materials

Lionel Gelebart

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

## To cite this version:

Lionel Gelebart. An accurate and robust FFT-based solver for transient diffusion in heterogeneous materials. Comptes Rendus. Mécanique, 2025, 353, pp.113-125. 10.5802/crmeca.281. cea-04937180

HAL Id: cea-04937180

https://cea.hal.science/cea-04937180v1

Submitted on 9 Feb 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Mécanique

Lionel Gélébart

An accurate and robust FFT-based solver for transient diffusion in heterogeneous materials

Volume 353 (2025), p. 113-125

Online since: 9 January 2025

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

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Research article / *Article de recherche*

# An accurate and robust FFT-based solver for transient diffusion in heterogeneous materials

*Un solveur FFT précis et robuste pour la diffusion transitoire dans les matériaux hétérogènes*

Lionel Gélébart<sup>●,a</sup>

**Abstract.** The purpose of the present letter is to propose an efficient, accurate and robust FFT-based solver for transient diffusion in heterogeneous materials with “realistic” BC, taking advantage of two recent advances in terms of boundary conditions and finite difference schemes to overcome their actual limitations (periodic BC and spurious oscillations). It is an essential step towards couplings between mechanics and other physics (such as the diffusion of species) through FFT-based solvers. Discrete Trigonometric Transform are used to implement non-periodic boundary conditions, and a finite difference (FD) scheme recently proposed by Finel is advantageously compared to the common hexahedral FD scheme. “Accurate” refers to two properties: accurate in term of locality with a small size of Finite Difference pencil to capture fluctuations around material heterogeneities, and accurate in term of precision with the absence of spurious spatial oscillations (at least in the reported cases with well-separated inclusions). The “robustness” is here associated to the stability of the solver, especially associated to the implicit time integration method. The description of the method focuses on thermal diffusion but applies to any other similar diffusion process (with the same type of parabolic equation). As a by-product, the FD scheme proposed by Finel is introduced in a more general (for mixing finite different schemes) and simple way (no introduction of FCC subgrids), extending its usage to any type of grid parity (not only even grids).

**Résumé.** Cette note propose un solveur FFT efficace, précis et robuste pour la diffusion transitoire dans des matériaux hétérogènes avec des conditions aux limites (CL) dites « réalistes ». Cette approche exploite deux avancées récentes : le traitement des conditions aux limites non périodiques et l'utilisation d'un nouveau schéma de différences finies limitant les oscillations parasites. Il s'agit d'une étape essentielle vers le couplage entre la mécanique et d'autres phénomènes physiques (comme la diffusion d'espèces) à l'aide de solveurs FFT. Les Transformées Trigonométriques Discrètes sont utilisées pour implémenter des conditions aux limites non périodiques, tandis qu'un schéma de différences finies (DF) récemment proposé par Finel est comparé avantageusement au schéma DF classique hexaédrique.

Le terme « précis » se réfère à deux aspects : (1) la précision locale, grâce à un schéma de DF de petite taille permettant de capturer les fluctuations autour des hétérogénéités du matériau, et (2) la précision en termes de fidélité, avec l'absence d'oscillations spatiales parasites (au moins dans les cas étudiés avec des inclusions bien séparées). La « robustesse » fait référence à la stabilité du solveur, notamment grâce une méthode d'intégration temporelle implicite. Bien que la description se concentre sur la diffusion thermique, la méthode s'applique à tout processus de diffusion similaire (régie par une équation parabolique du même type).

En complément, le schéma DF proposé par Finel est présenté de manière plus générale (pour mélanger différents schémas de différences finies) et simplifiée (sans introduction de sous-grilles FCC), élargissant également son usage à tout type de parité de grille (et pas seulement aux grilles paires).

<footnote><sup>a</sup> Université Paris Saclay, CEA, SRMA, 91191, Gif/Yvette, France</footnote>

<footnote>E-mail: lionel.gelebart@cea.fr</footnote>
**Keywords.** FFT, Heterogeneous materials, Transient diffusion, Finite difference, Tetrahedra.

**Mots-clés.** FFT, Matériaux hétérogènes, Diffusion transitoire, Différences finies, Tétraèdres.

*Manuscript received 23 October 2024, revised and accepted 18 December 2024.*

# 1. Introduction

Among many other applications, the development of an FFT-based solver for transient diffusion in heterogeneous material is essential to study the diffusion of species coming from the surface of the material and their coupling with mechanics. For such coupled simulations, the concentration of atoms induces a deformation of the crystal taken into account through an eigenstrain (stress-free strain) directly proportional to the concentration of atoms [1, 2]. Therefore, an accurate description of the concentration field is required. Besides, the boundary conditions generally used in this context are Dirichlet boundary conditions (applied concentration at the surface of the material, e.g. [3, 4]) or Neumann BC (applied flux, e.g. [5]). Hence, for such simulations, the finite element method is an appropriate numerical tool (e.g. [3–5]). However, the method is quite heavy in terms of memory or computation time, especially when considering 3D complex microstructures. FFT-based solvers for heterogeneous materials, initially proposed in the context of static mechanics by Moulinec and Suquet [6], provide an efficient alternative to standard FE solvers. The method introduces a homogeneous material on which the solution of an auxiliary problem is directly obtained by using Discrete Fourier Transforms (named FFT, for Fast Fourier Transform, in the following) through the application of the discrete Green operator. This efficient resolution of the auxiliary problem is used as a step of a fixed-point algorithm [6]. Described by similar partial differential equations, the method applies similarly to stationary diffusion problems on heterogeneous materials. Besides their efficiency, these solvers suffer from two important drawbacks in the present context. First, using Discrete Fourier Transforms, the boundary conditions (BC) are implicitly periodic. To circumvent this limitation, various teams have recently investigated the question for diffusion [7–9] and mechanics [10–12], taking advantage of Discrete Trigonometric Transforms (DTT), using both Finite Difference or Galerkin based spatial approximations. Second, the original method [6] suffers from spurious oscillations attributed to the Gibbs phenomenon. Using discrete Green operators built on Finite Difference schemes drastically reduces these artefacts. Among them, the most popular scheme, initially proposed by Willot and referred as the rotated scheme [13], is equivalent to the use of hexaedral linear finite elements with reduced integration (with a single Gauss point), noted HEX8R [14]. However, if the HEX8R is insensitive to Gibbs phenomenon, it comes with another kind of oscillations known as hourglass phenomenon [15]. Various techniques propose to stabilize hourglass [15], even in the context of FFT-based solvers [16], but they introduce an additional numerical parameter to be adjusted. Very recently, Finel proposed another type of Finite Difference scheme and its FFT-based implementation for static mechanics, completely removing the spurious oscillations observed with HEX8R [17], without any additional numerical parameter.

The purpose of the present letter is to propose an efficient, accurate and robust FFT-based solver for transient diffusion in heterogeneous materials with “realistic” BC, taking advantage of two recent advances in terms of boundary conditions and finite difference schemes to overcome actual limitations (periodic BC and spurious oscillations). “Accurate” refers to two properties: accurate in term of locality with small sizes of Finite Difference pencils to capture fluctuations around material heterogeneities, and accurate in term of precision with the absence of spurious spatial oscillations. The “robustness” is here associated to the stability of the solver, especially associated to the implicit time integration method. The description of the method focuses
on thermal diffusion but applies to any other similar diffusion process (with the same type of parabolic equation). As a by-product, the FD scheme proposed by Finel [17] is introduced in a more general (for mixing finite different schemes) and simple way (no introduction of FCC subgrids), extending its usage to any type of grid parity (not only even grids).

## 2. Method

### 2.1. Implicit FFT-based resolution with non-periodic boundary conditions

The volume equations governing the transient thermal diffusion problem reads:

$$ \begin{cases} \rho c_p \frac{\partial T}{\partial t} + \text{div}(\mathbf{q}) = 0 \\ \mathbf{q} = -\underline{K} \cdot \nabla T \end{cases} \qquad (1) $$

with $\rho$, $c_p$ and $\underline{K}$, the density, specific heat capacity and thermal conductivity; material parameters that can be heterogeneous over the domain $\Omega$, a parallelepiped domain of size $L \times L \times L$ m$^3$. $T$ and $\mathbf{q}$ are respectively the temperature and the heat flux. A set of “realistic” initial and boundary conditions is proposed. The initial temperature $T_0$ is assumed homogeneous except on the face $x=0$, where $T(x=0) = T^{\text{left}}$. The boundary conditions are Dirichlet BC on faces $x=0$ and $x=L$ with respectively $T(x=0) = T^{\text{left}}$ and $T(x=L) = T_0$, for any time. The four other faces are isolated (i.e. Neumann BC with $\mathbf{q} \cdot \mathbf{n} = 0$, $\mathbf{n}$ being the outer normal vector to the face).

Equation (2), with $dT$ as unknown field, discretizes Equation (1) in space and time, with a $\theta$-method in time, purely implicit when $\theta=1$ (Euler implicit) and purely explicit when $\theta=0$ (Euler explicit). The case $\theta=0.5$ corresponds to the implicit Cranck-Nicholson method. The method is implicit as soon as $\theta > 0$. The equation reads:

$$ \rho c_p \frac{dT}{dt} - \text{div}_D(\underline{K} \cdot \nabla_D (T^t + \theta dT)) = 0, \qquad (2) $$

with $T^t$ the temperature field, known at the beginning of the time step ($t$), $T^t + dT$, the temperature at the end of the time step ($t+dT$), and $dT$ the unknown field. The subscript (D) corresponds to a finite differences (FD) approximation for operators div and $\nabla$. The domain $\Omega$ is regularly discretized in small elements (or voxels). For the two FD schemes considered below, the temperature is defined at voxels' corners, as well as $\text{div}_D(\mathbf{q})$, and the heat flux $\mathbf{q}$ is defined at voxels' centres, as $\nabla_D T$. The grid of centres has a size $n^3$ and the grid of corners has a size $(n+1)^3$.

Introducing homogeneous quantities ($\rho c_p)_0$ and $k_0$ allows for defining a fix-point algorithm on $dT$:

$$ \frac{(\rho c_p)_0}{dT} dT - \theta k_0 \Delta_D dT = \frac{(\rho c_p)_0 - \rho c_p}{dT} dT + \text{div}_D(\underline{K} \cdot \nabla_D (T^t + \theta dT) - \theta k_0 \nabla_D (dT)) \qquad (3) $$

Actually, for a given field $dT$, the right hand-side (rhs) term can be evaluated in real space using the FD scheme (D). Noting $p$ the rhs term, and assuming periodic BC for $dT$ (non-periodic BC are considered later), the left hand-side (lhs) term is inverted using Discrete Fourier Transforms to obtain $dT$ as follows:

$$ \hat{dT} = \frac{\hat{p}}{\frac{(\rho c_p)_0}{dT} + \theta k_0 \|\xi_D\|^2} = \hat{G}_0(\xi) \hat{p} \qquad (4) $$

where $\xi_D$ are “modified” wave-vectors that are functions of the classical wave-vectors $\xi$. The expression of $\xi_D$ depends on the chosen FD scheme (D) (various expressions can be found for different FD schemes in [18]).
In the case of non-periodic boundary conditions, expression (4) still holds but with Discrete Fourier Transform performed on 2 (or 4) times extended signals, evaluated by Discrete Trigonometric Transforms (DTTs). The correspondence between the DTTs and the types of BC are reported in Appendix A. For a detailed description, readers are invited to consult various recent papers considering the extension of FFT-based solvers to non-periodic BC [7–9] and especially [7] written for various FD schemes, including the HEX8R scheme used below, and for any type of BC.

Regarding the BC applied in the present case, $dT$ is assigned to null Dirichlet BC (on faces $x=0$ and $x=L$), and null Neumann BC on the other faces. As explained in [7], the null Dirichlet BC on faces $x=0$ and $x=L$ are enforced during derivation in real space by assuming an anti-symmetric field $dT$ for the points located outside the domain (i.e. $x < 0$ and $x > L$). The equivalent in Fourier space is achieved by using the Discrete Sine Transform of type I (DST1). For the four other faces, the null flux is enforced during derivation by assuming a symmetric field $dT$ on the points located outside the domain (i.e. $y < 0$ and $y > L$, $z < 0$ and $z > L$). The equivalent in Fourier space is achieved using the Discrete Cosine Transform of type I (DCT1). Note that for the two FD schemes considered below (HEX8R and TETRA2), the gradient terms $\nabla_D(\cdot)$ of Equation (3) are evaluated at voxels’ centres from the temperature field defined at voxels’ nodes. Hence, their evaluation with these FD schemes do not rely on quantities located outside the domain. On the contrary, the divergence term $\text{div}_D(\cdot)$ of Equation (3) is evaluated (with HEX8R or TETRA2 scheme) at voxels’ nodes from the quantity $(K \cdot \nabla_D(T^t + \theta dT) - \theta k_0 \nabla_D(dT))$ defined at voxel’s centres. Hence, the evaluation of the $\text{div}_D(\cdot)$ term relies on a quantity defined outside the domain using appropriate symmetries. For the quantity $K \cdot \nabla_D(T^t + \theta dT) - \theta k_0 \nabla_D(dT)$, the symmetries are those of $\nabla_D(dT)$ (for each component of the gradient, the derivation in a given direction inverses the type of symmetry in that direction (symmetry $\leftrightarrow$ anti-symmetry) and let unchanged the type of symmetry in the two other directions, see [7]).

## 2.2. Focus on two finite difference schemes

The presentation above introduces the $\theta$-method in time without any precision on the FD scheme used in space. Focus is put on small size Finite Difference pencils: the widely adopted HEX8R scheme [13, 14] and the recently proposed tetrahedral-based scheme [17], denoted TETRA2 in the following (the “2” standing for the two tetrahedras used to define the scheme).

For the HEX8R scheme, the gradient is evaluated at a voxel’s centre from the 8 values at voxel’s corners and the divergence at a voxel’s corner from the 8 values defined at the 8 neighbouring voxel’s centres. In both cases, the finite difference derivation can be deduced from the approximated volume averages given below (where $\Omega_D = \Omega_H$ is an hexahedra depicted on Figure 1, whose vertices are whether voxels’ corners, for the temperature gradient, or the voxels’ centres for the flux divergence):

$$ \left\{ \begin{aligned} \nabla T_D &= \frac{1}{V_D} \iiint_{\Omega_D} \nabla T \, dv = \frac{1}{V_D} \iint_{\partial \Omega_D} T n \, ds \\ \text{div}_D(\boldsymbol{q}) &= \frac{1}{V_D} \iiint_{\Omega_D} \text{div}(\boldsymbol{q}) \, dv = \frac{1}{V_D} \iint_{\partial \Omega_D} \boldsymbol{q} \cdot n \, ds \end{aligned} \right. \qquad (5) $$

The approximation of integrals, on each face of $\partial\Omega_D$, by a discrete sum on nodal quantities allows defining the complete FD scheme. Note that in that case, the same finite difference scheme is used for gradient and divergence, the only difference comes from the location of the derived quantity (whether a grid of corners, for the temperature, or a grid of centres, for the flux). The usage of the HEX8R scheme in FFT-based methods is well documented and the expression of the associated modified wave-vectors, necessary for the FFT-based implementation, can be found in [13, 18].
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Definition of the volumes used for the different Finite Difference schemes (see Equation (5)).</figcaption></figure>

The TETRA2 scheme corresponds to a very recent proposition by Finel [17], in the context of numerical periodic homogenization for mechanics. Introduced in a rather complex way, a simplified description is proposed below. The main idea is that the temperature gradient is evaluated on two grids, both located at voxels' centres (same location as with the HEX8R scheme). Referring to a finite element perspective, they can be regarded as two Integration Point (IP) grids (instead of a single one for HEX8R). For each IP grid, the temperature gradient is evaluated with one of the two tetrahedral FD schemes, T1 and T2 (see Figure 1). As the two tetrahedras are defined within an hexaedra, this finite difference scheme has the same support as HEX8R (i.e. the same points are used).

In practice, the two temperature gradient fields, $\nabla_{T1} T$ and $\nabla_{T2} T$, are evaluated using a tetrahedral Finite difference (using Equation (5)) with $\Omega_D = \Omega_{T1}$ or $\Omega_D = \Omega_{T2}$, the two tetrahedras defined on Figure 1). Taking into account the two associated fluxes $q_{T1}$ and $q_{T2}$ (with $q_{T1} = -K \cdot \nabla_{T1} T$), the discrete flux divergence is defined from a mixture of the two fluxes, $q_{T1}$ and $q_{T2}$, as follows:

$$ \mathrm{div}_{T}(\boldsymbol{q}) = \frac{1}{2}(\mathrm{div}_{T1}(\boldsymbol{q}_{T2}) + \mathrm{div}_{T2}(\boldsymbol{q}_{T1})), \quad (6) $$

so that, the implicit Equation (2) becomes:

$$ \rho c_p \frac{dT}{dt} - \frac{1}{2} (\mathrm{div}_{T1}(K \cdot \nabla_{T2}(T^t + dT)) + \mathrm{div}_{T2}(K \cdot \nabla_{T1}(T^t + dT))) = 0, \quad (7) $$

and finally, the fix-point algorithm corresponding to Equation (3) reads:

$$ \frac{(\rho c_p)_0}{dt} dT - \theta k_0 \Delta_{T1T2} dT = \frac{(\rho c_p)_0 - \rho c_p}{dt} dT + \frac{1}{2} (\mathrm{div}_{T2}(K \cdot \nabla_{T1}(T^t + dT)) - \theta k_0 \nabla_{T1}(dT)) + \mathrm{div}_{T1}(K \cdot \nabla_{T2}(T^t + dT)) - \theta k_0 \nabla_{T2}(dT)) \quad (8) $$

As previously, evaluating in real space the r.h.s term of Equation (8), denoted $p$, from a given $dT$ field, is the first step of the fix-point algorithm, and the second step consists of inverting the $lhs$ term in Fourier space to determine a new $dT$, used to evaluate a new rhs term and so on. The notation $\Delta_{T1T2}$ stands for $1/2(\mathrm{div}_{T1}(\nabla_{T2}(\cdot)) + \mathrm{div}_{T2}(\nabla_{T1}(\cdot)))$. In Fourier space, the divergence and gradient operators reads:

$$ \begin{aligned} \widehat{\nabla_{T1} T} &= i \hat{T} \boldsymbol{\xi}_{T1} = i \hat{T} d_{T1} \exp \left( \frac{i(\xi_x + \xi_y + \xi_z)}{2} \right) \\ \widehat{\mathrm{div}_{T1}(q)} &= i \hat{q} \cdot \boldsymbol{\xi}_{T1} = i \hat{q} \cdot d_{T1} \exp \left( -\frac{i(\xi_x + \xi_y + \xi_z)}{2} \right) \end{aligned} \quad (9) $$

The definition of $d_{T1}$ can be found in [17] and it must be emphasized that $d_{T1}$ and $d_{T2}$ are complex conjugate ($d_{T1} = \overline{d_{T2}}$). Hence, the application of the operator $\Delta_{T1T2}$ in Fourier space
corresponds to a multiplication by the term $-\overline{d_{T1}} \cdot \overline{d_{T1}} = -\|d_{T1}\|^2 (= -\|d_{T2}\|^2)$. Finally, the rhs term of Equation (8) can be easily inverted, replacing $\|\xi_D\|^2$ by $\|d_{T1}\|^2$ in Equation (4).

**Remarks.** This presentation is rather different from the presentation proposed by Finel [17]. It is believed simplified (without the introduction of FCC subgrids) and offers a quite general and versatile framework for mixing finite differences operators. In addition, this description (as well as its implementation) is independent of the grid size parity (odd or even number of voxels), whereas Finel limits its usage to even grids.

## 2.3. Additional details

To reduce the number of implicit iterations, the initialization of the fix point algorithm, at each time step, consists of a linear extrapolation from the solution fields $dT_{n-2}$ and $dT_{n-1}$ evaluated at the last two times steps.

In the present implementation, the fix-point algorithm can be accelerated by an Anderson acceleration technique [19] as proposed in [20] for accelerating FFT-based solvers in static mechanics. However, the benefit of the acceleration was not so clear in the case presented in Section 4, and numbers of iterations in Table 1 are given for the classical fix-point algorithm.

<table><caption>**Table 1.** Number of iterations as a function of the step size and the explicit or implicit method (up to the time of interest $t = 300dT_{CFL}$)</caption><thead><tr><th></th><th>Number of steps</th><th>HEX8R</th><th>TETRA2</th></tr></thead><tbody><tr><td>Explicit dt = dt<sup>*</sup><sub>CFL</sub>/2</td><td>600</td><td>600<sup>*</sup></td><td>600<sup>*</sup></td></tr><tr><td>Implicit dt = dt<sup>*</sup><sub>CFL</sub>/2</td><td>600</td><td>2876</td><td>3257</td></tr><tr><td>Implicit dt = dt<sup>*</sup><sub>CFL</sub></td><td>300</td><td>2406</td><td>2712</td></tr><tr><td>Implicit dt = 2dt<sup>*</sup><sub>CFL</sub></td><td>150</td><td>2246</td><td>2584</td></tr></tbody></table>

The choice made for the reference thermal conductivity comes from the proposition of Moulinec for isotropic materials [6]: $k_0 = 1/2(\min(k) + \max(k))$. If materials are anisotropic, an “isotropized” thermal conductivity can be taken into account. Note that the choice of $k_0$ does not modify the solution but the convergence rate of the algorithm.

For the sake of brevity in Section 4, only the conductivity is considered heterogeneous and ($\rho c_p$) is assumed homogeneous. Note that simulations performed with heterogeneous ($\rho c_p$) demonstrate similar results regarding spurious oscillations.

Finally, the relative error between two successive iterates is used as convergence criterion with a low threshold $\varepsilon_0$ to ensure the quality of the results:

$$ \frac{\|dT^k - dT^{k-1}\|_2}{\|dT^k\|_2} < \varepsilon_0 = 10^{-6} \qquad (10) $$

## 3. Validation on a homogeneous case

The algorithm and its implementation is partly validated on the homogeneous case for which the analytical solution is known. Considering the initial and boundary conditions described in Section 2.1, the 3D problem reduces to a 1D problem whose solution reads:

$$ T(x, t) = \left(1 - \frac{x}{L}\right) T^{\text{left}} + \sum_{k \in \mathbb{N}^*} -\frac{2T^{\text{left}}}{k\pi} \exp\left(-\left(\frac{k\pi}{L}\right)^2 \frac{K}{\rho c_p} t\right) \sin\left(\frac{k\pi}{L}x\right) \quad (11) $$

<footnote>* For the explicit method, each iteration is a direct evaluation (no iterative back and forth in Fourier space).</footnote>
The parameters used for the validation are $L = 1$ m, $K = 10 \text{ W} \cdot \text{m}^{-1} \cdot \text{K}^{-1}$, $\rho = 1000 \text{ kg} \cdot \text{m}^{-3}$, $c_p = 1000 \text{ J} \cdot \text{kg}^{-1} \cdot \text{K}^{-1}$, $T^{\text{left}} = 100$ K and $T^0 = 0$ K. For the numerical simulation, the 3D cubic domain is discretized in $32^3$ voxels (even if the solution is known invariant in directions $y$ and $z$).

### 3.1. Stability of the methods

The numerical experiments performed with the HEX8R and TETRA2 schemes combined with the explicit Euler method ($\theta = 0$) were stable as soon as $dt \le dt_{\text{CFL}}/2$, with $dt_{\text{CFL}} = \rho c_p dx^2/K$ and $dx$ the voxel dimension, and systematically stable with the implicit Euler method ($\theta = 1$). These observations are consistent with mathematical results obtained with an equivalent 1D algorithm, a $\theta$-method in time and centred second derivative in space (i.e. $d^2T/dx^2 \sim (T_{n+1} + T_{n-1} - 2T_n)/dx^2$). Actually, considering the loading and the homogeneous material, the solution is uniaxial and, in that context, both 3D schemes (HEX8R and TETRA2) turns out equivalent to this standard and well-known 1D algorithm that is conditionally stable for $\theta = 0$ (stable if $dt \le dt_{\text{CFL}}/2$), and unconditionally stable if $\theta \ge 1/2$.

### 3.2. Temperature profiles comparison

The times of interest for that comparison are $2dt_{\text{CFL}}$, $6dt_{\text{CFL}}$ and $18dt_{\text{CFL}}$, and the considered time steps are $dt_{\text{CFL}}/2$, the limit of stability for the explicit method, and $dt_{\text{CFL}}$, only accessible to the implicit method. Figure 2(left) plots the temperature profiles comparing analytical and numerical profiles. As expected from the remark in the previous subsection, HEX8R and TETRA2 schemes provide exactly the same results ("+" and "x" symbols are superimposed), in good agreement with the analytical profiles. The normalized error between analytical and numerical results, $(T_{\text{an.}} - T_{\text{num.}})/T^{\text{left}}$, is plotted on Figure 2(right). For time $t = 2dt_{\text{CFL}}$, results for the explicit method ($\theta = 0$), used with $dt = dt_{\text{CFL}}/2$ (its limit of stability), exhibits an error of ~6% that reduces to ~3% with the implicit method ($\theta = 1$). With a larger time of interest, $t = 6dt_{\text{CFL}}$ (dashed lines), the error decreases for both methods but the implicit scheme remains more precise (still with $dt = dt_{\text{CFL}}/2$). Using a twice-larger step size (i.e. $dt = dt_{\text{CFL}}$) the explicit method provides a similar precision w.r.t. the explicit method (with a twice-smaller step). As the purpose is more the validation of the solver and its implementation than a detailed comparison between the two methods (explicit vs implicit), widely explored in handbooks for this simple 1D case, the next step focuses on the heterogeneous case, for which HEX8R and TETRA2 schemes are no longer equivalent.

## 4. Application to the heterogeneous case

A simple case of application is proposed with a single spherical inclusion of radius 0.3 m at the position (0.15 m, 0.4 m, 0.6 m) within a cubic domain of size $1 \times 1 \times 1 \text{ m}^3$, discretized in $32 \times 32 \times 32$ voxels. The mass density and specific heat capacity are assumed homogeneous (same values as previously), while the thermal conductivity is set to $10^4 \text{ W} \cdot \text{m}^{-1} \cdot \text{K}^{-1}$. Compared to the conductivity of the matrix, $10 \text{ W} \cdot \text{m}^{-1} \cdot \text{K}^{-1}$, the contrast is quite high (1000). As spurious oscillations generally correlate with high contrasts, this case should allow for discriminating the two FD schemes, which were equivalent in the homogeneous case used in the previous section.

### 4.1. Stability of the method

The stability condition for the explicit method ($\theta = 0$) has been explored empirically on this example: whether combined with the HEX8R or the TETRA2 scheme, the method is stable as
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Transient diffusion in a homogeneous material. (left) Numerical temperature profiles compared to analytical solution for $t = 2dt_{CFL}$, $6dt_{CFL}$ and $18dt_{CFL}$. (right) Normalized error $(T_{an.} - T_{num.})/T^{left}$, for time $t = 2dt_{CFL}$ with continuous lines, and $t = 6dt_{CFL}$ with dashed lines. In the legend, $dt = 1/2$ stands for $dt = dt_{CFL}/2$, and $dt = 1$ for $dt = dt_{CFL}$.</figcaption></figure>

soon as the step size $dt \le dt_{CFL}^*/2$, with $dt_{CFL}^* = \rho c_p dx^2/\max(K)$, and unstable otherwise. This condition also proved valid for simulations with a low conductivity inclusion. Note that this definition of $dt_{CFL}^*$ allows for using the same stability condition than the homogeneous case with $dt_{CFL}$. For heterogeneous $\rho c_p$ the condition proved valid for $dt_{CFL}^* = \min(\rho c_p) dx^2/K$. For both $K$ and $\rho c_p$ heterogeneous, a similar definition for $dt_{CFL}^*$, if exists, is still under investigation. Finally, all the tests performed with the implicit method ($\theta = 1$) remained stable.

## 4.2. Temperature fields comparison

For that comparison, the time of interest is $t = 300dt_{CFL}^*$, chosen to easily discriminate the two schemes. Figure 3 displays the different temperature fields for the two schemes HEX8R (up) and TETRA2 (bottom), obtained with the explicit ($\theta = 0$) or implicit ($\theta = 1$) method (respectively left and right side). Explicit simulations have been performed using $dt = dt_{CFL}^*/2$ and implicit simulations with $dt = 2dt_{CFL}^*$. The HEX8R scheme induces important spurious oscillations, whether using an implicit or explicit method. Importantly, an overshoot is observed on the colorbar: the temperature field can be locally higher than the applied temperature $T^{left} = 100$ K. On the other hand, the TETRA2 scheme used with the implicit or explicit method, do not exhibit any sign of spurious oscillations or overshoot: the field is smooth within each material and the transition is sharp at the interface between the two materials. Finally, the quasi-identical fields obtained with the explicit method and the implicit method with 4 times larger step demonstrates the interest of the implicit method in that case (especially for quite large times of interest).

As a conclusion regarding the quality of the solution, using the TETRA2 scheme appears a relevant choice to avoid any artefact such as spurious oscillations and overshoot. This is a crucial point especially when considering the solution field as an input for another coupled physics (mechanics or phase transition, for example).

## 4.3. Performance

The purpose of this subsection is not to provide a deep analysis but rather a few elements of discussion on the choice of the step size, or between explicit or implicit method. Table 1 gathers the cumulated number of iterations up to the time of interest $t = 300dt_{CFL}^*$ with different time steps for the implicit method (all of them generate almost identical fields, as observed in Figure 3).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Temperature field distribution with the HEX8R scheme (up) and the TETRA2 scheme (bottom), using the explicit method ($\theta = 0$, $dt = dt_{CFL}/2$, left)), or the implicit method ($\theta = 1$, $dt = 2dt_{CFL}$, right)). Note the oscillations and overshoot (values higher than 100 K on the colorbar).</figcaption></figure>

Considering the choice between HEX8R and TETRA2, the number of iterations is a bit higher for TETRA2. In addition, each iteration of the TETRA2 scheme is a bit more complex as the flux is evaluated on two grids of centres (see Section 2). However, the price to pay, regarding the important enhancement of the solution quality, is rather low. Considering the choice of the step size for implicit algorithms: it is worth noting that while the step size is multiplied by 2 the number of iterations is only divided by a coefficient between 1.02 and 1.2. Actually, when increasing the step size, the number of fix-point iterations, required to solve each step, increases (especially at the beginning of the simulation). Finally, considering the choice between explicit and implicit algorithm, the difference in terms of iterations is rather important. In addition, an iteration of the explicit algorithm has a much lower cost as it consists of a direct evaluation of $\text{div}(q)$ in real space and does not require any back and force in Fourier space. From that point of view, and for that
given example (material, loading and time of interest), the explicit method appears much more performant than the implicit scheme.

**5. Example on a larger unit-cell**

The present section is just an example showing that, despite of a non-optimized *python* implementation for this proof of concept, the method applies straightforwardly to a larger microstructure that consists of 125 randomly distributed spheres with a volume fraction of 22%. The size of the unit-cell is $1 \times 1 \times 1 \text{ m}^3$, discretized in $128 \times 128 \times 128$ voxels. Microstructure is generated with periodicity condition and periodic conditions are considered in the transverse directions (wrt to the temperature gradient direction) and the same conditions are applied for the upper and lower face (100 K and 0 K homogeneous applied temperature). Initial temperature is 0 K. The mass density and specific heat capacity are assumed homogeneous (same values as previously), while the thermal conductivity is set to $10^4 \text{ W} \cdot \text{m}^{-1} \cdot \text{K}^{-1}$ for the inclusion and $100 \text{ W} \cdot \text{m}^{-1} \cdot \text{K}^{-1}$ for the matrix. The temperature field displayed on Figure 4 is evaluated after $5000d t_{\text{CFL}}^*$. A second simulation considers insulating inclusions with a thermal conductivity of $1 \text{ W} \cdot \text{m}^{-1} \cdot \text{K}^{-1}$. The temperature field is evaluated at the same time (corresponding to $50d t_{\text{CFL}}^*$ in that case). On this more realistic application, the temperatures fields are still very smooth and comparison clearly exhibits the influence of the inclusions, whether conducting or insulating, on the propagation of the temperature with the composite material.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Example exhibiting the influence of conducting (left) or insulating (right) spherical inclusions within a matrix on the temperature field after the same propagation time (the simulated unit-cell is shown at the centre).</figcaption></figure>

Important remark: it must be mentioned that, in this example, the inclusions are well separated from each others. Various tests, not reported in this short note, have revealed the existence of spurious oscillations in cases where inclusions are “touching” each other(s). If the amount of oscillation is still less than HEX8R scheme, this point must be kept in mind and needs further investigations.

**6. Conclusion and prospects**

As a conclusion of this letter, it is demonstrated that combining recent developments using DTTs, to apply non-periodic boundary conditions, with a recent Finite Difference scheme (noted
TETRA2), allows for taking advantage of the efficiency of FFT-based solvers for simulating tran-
sient diffusion in heterogeneous materials with realistic boundary conditions. As demonstrated,
the method is robust in time (using an implicit solver), accurate in space as the solution fields do
not exhibit any spurious oscillations (at least if the material consists of a matrix with well sepa-
rated inclusions), and accurate in term of “locality” as the TETRA2 finite difference scheme only
considers the 8 points of an hexahedra to approximate the divergence or the gradient. In that
sense, it is as local as the HEX8R scheme but do not suffer from spurious oscillations originat-
ing from the hourglass phenomenon. To the author’s opinion, the presentation of the TETRA2
scheme, proposed by Finel in his original paper [17], has been simplified, generalized (towards
other mixing of finite differences) and a bit extended (its usage is not limited to even grids). From
the numerical experiments in Section 4, using the TETRA2 scheme is a bit less efficient than the
HEX8R scheme, however the quality enhancement is so important that theTETRA2 choice should
be favoured (author’s opinion). Regarding the choice of explicit or implicit method, and step size
for the implicit method, the question has been addressed for a specific case of application (homo-
geneous ρc_p, high conductivity contrast, large time of interest) and the explicit method proved
more efficient than the implicit one. However, a strategy mixing Euler explicit and implicit to-
gether with an adaptive step size will probably be an efficient way to take benefit from the advan-
tages of the two methods.

As mentioned in Section 5, for matrix inclusion microstructures with inclusions touching each
other, spurious oscillations appear on the temperature field with the TETRA2 scheme. If less
important than with the HEX8R scheme, this point needs further investigations. The use of
composite voxels is generally an efficient way to smooth oscillations is such situations [21].

As a future prospect, the present FFT-based solver (θ-method in time and TETRA2 scheme in space) can be considered as a reliable and essential component for coupling transient diffusion with mechanics. Important applications in the context of materials science consider the diffusion of species strongly coupled with mechanics through a stress-free strain directly proportional to the concentration of diffusing atoms (e.g. [3–5]). Such couplings, generally performed with standard Finite Element methods, could now take benefit from the present efficient, robust and accurate FFT-based solver adapted to heterogeneous materials and realistic boundary conditions.

Finally, an implementation of this framework is in progress in the FFT-based code AMI-TEX [22] in order to take advantage of massively parallel computers and deal with large scale simulations.

**Declaration of interests**

The authors do not work for, advise, own shares in, or receive funds from any organization
that could benefit from this article, and have declared no affiliations other than their research
organizations.

**Appendix A. DTTs and boundary conditions**

Table 2 gathers the different discrete symmetry types with the corresponding types of DTT
(according to the names defined in the FFTW library [23]) and the associated types of boundary
conditions.

For the definition of the symmetry TXTY:

* T has the value W, if the symmetry axis are located on the last points, or H, if the axis is located a half element just after the last points. In the present paper, the temperature being located at voxels' corners, only W types of symmetry are considered,
* X has the value S or A, whether the signal is Symmetric or Anti-symmetric on the left point,

* Y has the value S or A, whether the signal is Symmetric or Anti-symmetric on the right point.

Finally, an Anti-symmetry is associated to a Dirichlet BC and a Symmetry to a Neumann BC.

<table><caption>**Table 2.** DTTs and boundary conditions</caption><thead><tr><td>Discrete symmetry type (TXTY)</td><td>Discrete Trigo. Transf. [23]</td><td>Boundary condition (left/right)</td></tr></thead><tbody><tr><td>WSWS</td><td>DCT1</td><td>Neumann/Neumann</td></tr><tr><td>HSHS</td><td>DCT2</td><td></td></tr><tr><td>WSWA</td><td>DCT3</td><td>Neumann/Dirichlet</td></tr><tr><td>HSHA</td><td>DCT4</td><td></td></tr><tr><td>WAWA</td><td>DST1</td><td>Dirichlet/Dirichlet</td></tr><tr><td>HAHA</td><td>DST2</td><td></td></tr><tr><td>WAWS</td><td>DST3</td><td>Dirichlet/Neumann</td></tr><tr><td>HAHS</td><td>DST4</td><td></td></tr></tbody></table>

## References

[1] F. Larché, J. W. Cahn, "A linear theory of thermochemical equilibrium of solids under stress", *Acta Metall.* **21** (1973), no. 8, p. 1051-1063.

[2] F. C. Larché, J. W. Cahn, "Overview no. 41 the interactions of composition and stress in crystalline solids", *Acta Metall.* **33** (1985), no. 3, p. 331-357.

[3] A. Villani, E. P. Busso, K. Ammar, S. Forest, M. G. D. Geers, "A fully coupled diffusional-mechanical formulation: numerical implementation, analytical validation, and effects of plasticity on equilibrium", *Arch. Appl. Mech.* **84** (2014), no. 9, p. 1647-1664.

[4] E. Martínez-Pañeda, A. Golahmar, C. F. Niordson, "A phase field formulation for hydrogen assisted cracking", *Comput. Methods Appl. Mech. Eng.* **342** (2018), p. 742-761.

[5] W. Jbara, V. Maurel, K. Ammar, S. Forest, "Effect of free surface, oxide, and coating layers on rafting in $\gamma - \gamma'$ superalloys", in *Superalloys 2024* (J. Cormier, I. Edmonds, S. Forsik, P. Kontis, C. O'Connell, T. Smith, A. Suzuki, S. Tin, J. Zhang, eds.), Springer Nature Switzerland, Cham, 2024, p. 694-704.

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

[7] L. Gélébart, "FFT-based simulations of heterogeneous conducting materials with combined non-uniform neumann, periodic and dirichlet boundary conditions", *Eur. J. Mech. - A/Solids* **105** (2024), article no. 105248.

[8] J. Paux, L. Morin, L. Gélébart, "A discrete sine-cosine transforms Galerkin method for the conductivity of heterogeneous materials with mixed Dirichlet/Neumann boundary conditions", *Int. J. Numer. Methods Eng.* **126** (2025), article no. e7615.

[9] L. Risthaus, M. Schneider, "Imposing different boundary conditions for thermal computational homogenization problems with FFT- and tensor-train-based Green's operator methods", *Int. J. Numer. Methods Eng.* **125** (2024), no. 7, article no. e7423.

[10] L. Gélébart, L. Morin, J. Paux, "FFT-based mechanical simulations of heterogeneous materials with combined non-uniform neumann, periodic and dirichlet boundary conditions", *Comput. Mech.* (2024), revision in progress.

[11] L. Risthaus, M. Schneider, "Imposing Dirichlet boundary conditions directly for FFT-based computational micromechanics", *Comput. Mech.* **74** (2024), p. 1089-1113.

[12] J. Paux, L. Morin, L. Gélébart, A. M. Amadou Sanoko, "A discrete sine-cosine based method for the elasticity of heterogeneous materials with arbitrary boundary conditions", *Comput. Methods Appl. Mech. Eng.* **433** (2025), article no. 117488.

[13] F. Willot, "Fourier-based schemes for computing the mechanical response of composites with accurate local fields", *C. R. Méc.* **343** (2015), no. 3, p. 232-245.

[14] M. Schneider, D. Merkert, M. Kabel, "FFT-based homogenization for microstructures discretized by linear hexahedral elements", *Int. J. Numer. Meth. Eng.* **109** (2017), no. 10, p. 1461-1489.

[15] T. Belytschko, J. S.-J. Ong, W. K. Liu, J. M. Kennedy, "Hourglass control in linear and nonlinear problems", *Comput. Methods Appl. Mech. Eng.* **43** (1984), no. 3, p. 251-276.
[16] M. Schneider, "Voxel-based finite elements with hourglass control in fast fourier transform-based computational homogenization", *Int. J. Numer. Methods Eng.* **123** (2022), no. 24, p. 6286-6313.

[17] A. Finel, "A fast and robust discrete FFT-based solver for computational homogenization", 2024, submitted, https://arxiv.org/pdf/2405.11168v1.

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

[19] D. G. Anderson, "Iterative procedures for nonlinear integral equations", *J. ACM (JACM)* **12** (1965), no. 4, p. 547-560.

[20] Y. Chen, L. Gélébart, C. Chateau, M. Bornert, C. Sauder, A. King, "Analysis of the damage initiation in a sic/sic composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography", *Int. J. Solids Struct.* **161** (2019), p. 111-126.

[21] R. Charière, A. Marano, L. Gélébart, "Use of composite voxels in FFT based elastic simulations of hollow glass microspheres/polypropylene composites", *Int. J. Solids Struct.* **182-183** (2020), p. 1-14.

[22] L. Gélébart, 2022, https://amitexfftp.github.io/amitex/index.html.

[23] M. Frigo, S. G. Johnson, "The design and implementation of FFTW3", *Proc. IEEE* **93** (2005), no. 2, p. 216-231.