# A FFT-based numerical scheme for the transient conductivity of heterogeneous materials with non-periodic boundary conditions

Abdoul Magid Amadou Sanoko, Simon Essongue, Lionel Gélébart, Lucas Lapostolle, Léo Morin, Joseph Paux

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

## ► To cite this version:

Abdoul Magid Amadou Sanoko, Simon Essongue, Lionel Gélébart, Lucas Lapostolle, Léo Morin, et al. A FFT-based numerical scheme for the transient conductivity of heterogeneous materials with non-periodic boundary conditions. European Journal of Mechanics - A/Solids, 2025, 113, pp.105680. 10.1016/j.euromechsol.2025.105680. hal-05067769

HAL Id: hal-05067769

https://hal.science/hal-05067769v1

Submitted on 14 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
Full length article

# A FFT-based numerical scheme for the transient conductivity of heterogeneous materials with non-periodic boundary conditions

Abdoul Magid Amadou Sanoko a,b, Simon Essongue a,b, Lionel Gélébart c, Lucas Lapostolle a,b,
Léo Morin a,b,*, Joseph Paux d

## ARTICLE INFO

### Keywords:

FFT-based solvers

Mixed Dirichlet/Neumann boundary conditions

Discrete sine-cosine transforms

Transient conductivity

## ABSTRACT

The aim of this work is to develop FFT-based solvers for transient diffusion in heterogeneous materials subjected to non-periodic (Dirichlet/Neumann) boundary conditions. We focus on a problem of thermal conductivity and derive a theta-method which includes an implicit solver for transient thermal conductivity in heterogeneous materials. The method is based on a fixed-point iterative solution of an auxiliary problem obtained by a Galerkin discretization using an approximation space based on mixed sine-cosine series. The solution field is decomposed as a known term verifying the boundary conditions and a fluctuation (unknown) term described by appropriate sine-cosine series. The solution of the auxiliary problem involves discrete sine-cosine transforms, of type I and III, which makes the solver rely on the computational complexity of fast Fourier transforms. The method is applied to the prediction of transient thermal fields in a composite material subjected to non-periodic boundary conditions.

## 1. Introduction

This work is concerned with the development of efficient FFT-based numerical schemes for diffusion equations in heterogeneous materials subjected to non-periodic Dirichlet/Neumann boundary conditions, as a first step towards the development of fast FFT-based multiphysics solvers. FFT-based methods for heterogeneous materials have been introduced in the seminal work of Moulinec and Suquet (1998), to compute the local and overall mechanical fields in composite materials subjected to periodic boundary conditions, mainly in the context of micromechanics. The advantages of this class of spectral solvers over standard finite element methods (FEM) are that (i) calculations are based on the computational complexity of fast Fourier transforms in $O(N \log N)$ (which is lower than in most of matrix-based methods), (ii) it is image-based (obtained from imagery techniques such as scanning electron microscopy (SEM) or tomography) so it does not require meshing operations and (iii) it can be suited for parallel implementation (see e.g. Gélébart (2022)). The initial fixed-point algorithm proposed by Moulinec and Suquet (1998) for the iterative procedure, i.e. the so-called basic scheme, has been improved through the use of accelerated iterative schemes (Eyre and Milton, 1999; Michel et al., 2000; Monchiet and Bonnet, 2012; Moulinec and Silva, 2014; Kabel et al., 2014), to

overcome some convergence issues of the method for highly-contrasted materials. FFT-based methods have therefore been successfully used to predict the local and macroscopic fields in numerical applications including elasticity (Moulinec and Suquet, 1998), J2-plasticity (Moulinec and Suquet, 1998), crystal (elasto)viscoplasticity (Lebensohn, 2001), dislocation-mediated plasticity (Brenner et al., 2014; Bertin et al., 2015), conductivity (Eyre and Milton, 1999), piezoelectricity (Brenner, 2009), porous ductile solids (Bilger et al., 2005; Paux et al., 2018), among others (see the review papers of Schneider (2021) and Lucarini et al. (2021)).

FFT-based methods are therefore adapted to coupled problems in heterogeneous microstructures due to the absence of resolution of a linear system, and have been successfully applied to the resolution of transient multiphysics problems, including electro-mechanics for the response of polycrystalline ferroelectric ceramics (Vidyasagar et al., 2017), vacancy diffusion-mediated crystal plasticity (Chakraborty et al., 2023), corrosion-driven damage in porous media (Pundir et al., 2023) and chemo-mechanical degradation of electrodes (Zarzoso et al., 2025).

However, the method proposed by Moulinec and Suquet (1998) is by essence restricted to periodic problems, as it relies on the use of Fourier series (and associated transforms), which makes it not suitable

<footnote>a Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France</footnote>

<footnote>b Arts et Metiers Institute of Technology, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France</footnote>

<footnote>c Université Paris-Saclay, CEA, SRMA, 91191, Gif/Yvette, France</footnote>

<footnote>d Université Sorbonne Paris Nord, Laboratoire des Sciences des Procédés et des Matériaux, LSPM, CNRS, UPR 3407, F-93430, Villetaneuse, France</footnote>

<footnote>* Corresponding author at: Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France.
E-mail address: leo.morin@u-bordeaux.fr (L. Morin).</footnote>
in several situations including transient (diffusion or elastodynamic) problems for which non-periodic boundary conditions are of prime importance (see e.g. Sancho et al. (2023) for elastodynamic calculations). Several methods have been proposed to extend Moulinec and Suquet (1998)'s method to non-periodic boundary conditions:

* Buffer-type approaches, introduced by Gélébart (2020), are based on the introduction of a buffer zone surrounding the domain of interest. The complete cell, composed of the domain of interest plus the buffer zone, can be periodic so this technique relies on standard FFT packages. It allows to apply Dirichlet or Neumann boundary conditions (Kaptchouang and Gélébart, 2022). This approach requires (i) a much larger unit cell than the microstructure studied because of the (elastic) buffer zone as well as (ii) additional internal iterations to find the eigendisplacement field at the boundary of the buffer zone. Mention has to be made of the recent work of Zecevic and Lebensohn (2025) in which non-periodic boundary conditions are applied by extending the domain of interest to satisfy the periodicity and adding a penalization with augmented lagrangian.

* Mirror-cell approaches, employed by Grimm-Strele and Kabel (2021) in elasticity problems and Monchiet and Bonnet (2024) in conductivity problems, take advantage of the equivalence between particular symmetries of the unit cell and boundary conditions on a portion of the unit cell, and allow to apply uniform Dirichlet and/or Neumann boundary conditions. These approaches use implicitly discrete sine-cosine transforms.

* Body force field methods, introduced in the work of Sancho et al. (2023), are based on the introduction of an artificial body force field emulating time-dependent Dirichlet boundary conditions. This approach has been applied to transient elasto-dynamic problems using a beta-Newmark approach.

* Discrete sine-cosine transforms based approaches use a description of the solution field through sine-cosine series (instead of Fourier series) which allows to take into account arbitrary (non-uniform) Dirichlet/Neumann boundary conditions on entire faces. These approaches have been investigated in conductivity (Risthaus and Schneider, 2024b; Gélébart, 2024; Morin and Paux, 2024; Paux et al., 2025a) and elasticity problems (Risthaus and Schneider, 2024c,a; Paux et al., 2025b). In this approach, the solution field of the elliptic problem is split into a known term verifying the Dirichlet/Neumann boundary conditions and an unknown fluctuation term described by an appropriate trigonometric series with no contribution on the boundary. With such description, the heterogeneous problem can be solved iteratively using discrete sine and cosine transforms which can be computed using FFT packages and have thus a computational complexity in $O(N \log N)$.

The aim of this work is to extend (Paux et al., 2025a)'s method for the stationary conductivity of heterogeneous materials subjected to non-periodic Dirichlet/Neumann boundary conditions to the transient case, using a Galerkin approach and discrete sine-cosine transforms. In transient problems, the non-periodic FFT-based approach is expected to be very numerically efficient (in comparison to standard FEM approaches) because (i) in the case of implicit schemes, it does not involve any resolution of a linear system at each time step, but only fast Fourier transforms and (ii) it requires less memory storage. Mention has to be made to the recent work of Gélébart (2025) which developed FFT-based solver for transient diffusion in heterogeneous materials in the context of discrete Green operators, i.e. using a finite difference scheme for the spatial derivatives. The main difference between the present work and Gélébart (2025)'s approach is the use of a Galerkin approach together with spectral derivative rules, which leads to a different Green operator solution of the auxiliary problem (as shown in the stationary case (Paux et al., 2025a)).

The paper is organized as follows. In Section 2, the framework of sine-cosine series (and the associated discrete transforms) is introduced. The problem transient conductivity in heterogeneous media subjected to non-periodic boundary conditions is then solved in Section 3, using an iterative solver of an auxiliary problem that is discretized using a Galerkin method in which the unknown field is described by sine-cosine series. In Section 4 the method is validated using an analytical solution in a homogeneous case. Finally the method is assessed in Section 5 through several numerical examples.

## 2. The framework of sine-cosine series and the associated discrete transforms

We recall in this section several important results pertaining to sine-cosine series and the associated discrete transforms, which will be necessary for the description of the solution fields of the thermal conductivity problem under mixed Dirichlet/Neumann boundary conditions. Indeed, sine-cosine series intrinsically describe functions than can verify null-Dirichlet and/or null-Neumann conditions on any boundaries and will thus be used to model the fluctuations solution of the transient conductivity equations considered hereafter. This framework is therefore necessary for the description of test functions used in the Galerkin discretization of the local problem.

### 2.1. The 1-d case

**Sine-cosine series.** We consider a function $f$ defined in the interval $[0, L_x]$ which verifies null Dirichlet and/or null Neumann boundary conditions. Such function can be described by half-range or quarter-range sine-cosine series as

$$f(x) = \sum_{i=0}^{+\infty} \alpha_i F_i g(k_i x), \quad \forall x \in [0, L_x], \quad (1)$$

where $k_i$ correspond to some “frequency” parameter, $g$ is a cosine or a sine function (depending on the type of boundary conditions) and $F_i$, $i \ge 0$ are the sine-cosine series coefficients given by

$$F_i = \frac{2}{L_x} \int_{0}^{L_x} f(x)g(k_i x)dx, \quad \forall i \ge 0. \quad (2)$$

The associated partial series of order $N$ (for every positive integer $N$) is classically by

$$f_N(x) = \sum_{i=0}^{N} F_i g(k_i x), \quad \forall x \in [0, L_x]. \quad (3)$$

A summary of all functions and coefficients is given in Table 1 to cover all types of boundary conditions on $x = 0$ and $x = L_x$ (Dirichlet-Dirichlet (DD), Neumann-Neumann (NN), Dirichlet-Neumann (DN) and Neumann-Dirichlet (ND)).

Derivative rules associated with (1) can be written as

$$\frac{df}{dx}(x) = \sum_{i=0}^{+\infty} \xi_i \alpha_i F_i \bar{g}(k_i x), \quad \forall x \in [0, L_x], \quad (4)$$

$$\frac{d^2 f}{dx^2}(x) = \sum_{i=0}^{+\infty} -(\xi_i)^2 \alpha_i F_i g(k_i x), \quad \forall x \in [0, L_x]. \quad (5)$$

In these expressions, $\bar{g}$ is the dual function of $g$ (i.e. a sine function if the initial function is a cosine function and vice versa) and $\xi_i$ is the coefficient arising from the derivation (see Table 1).

**Discrete sine-cosine transforms.** The computation of the sine-cosine coefficients (2) associated to a sine-cosine series (1) can be calculated numerically using the discrete sine transform (DST) and the discrete cosine transform (DCT), which are similar to the standard Discrete Fourier Transform (DFT) used for the computation of the coefficients of the Fourier series for periodic functions. An efficient computation of the coefficients of the DFT, DST and DCT is classically done using fast Fourier transforms (FFT) (Frigo and Johnson, 1998).
<table><caption>**Table 1**
Functions and coefficients defining the four mixed Dirichlet/Neumann boundary conditions on the faces $x=0$ and $x=L_x$.</caption><thead><tr><th></th><th>Boundary conditions</th><th>g and g&#773;</th><th>k<sub>i</sub></th><th>&alpha;<sub>i</sub></th><th>&xi;<sub>i</sub></th></tr></thead><tbody><tr><td rowspan="2">DD</td><td>f(x, y, z) = 0 for x = 0</td><td>g(x) = sin(x)</td><td rowspan="2">k<sub>i</sub> = <sup>i</sup>&frasl;<sub>L<sub>x</sub></sub>&pi;</td><td rowspan="2">&alpha;<sub>i</sub> = 1</td><td rowspan="2">&forall;i &ge; 0</td></tr><tr><td>f(x, y, z) = 0 for x = L<sub>x</sub></td><td>g&#773;(x) = cos(x)</td></tr><tr><td rowspan="2">NN</td><td><sup>&part;f</sup>&frasl;<sub>&part;x</sub>(x, y, z) = 0 for x = 0</td><td>g(x) = cos(x)</td><td rowspan="2">k<sub>i</sub> = <sup>i</sup>&frasl;<sub>L<sub>x</sub></sub>&pi;</td><td rowspan="2">&alpha;<sub>0</sub> = 1/2</td><td rowspan="2">&forall;i &ge; 1</td></tr><tr><td><sup>&part;f</sup>&frasl;<sub>&part;x</sub>(x, y, z) = 0 for x = L<sub>x</sub></td><td>g&#773;(x) = sin(x)</td></tr><tr><td rowspan="2">DN</td><td>f(x, y, z) = 0 for x = 0</td><td>g(x) = sin(x)</td><td rowspan="2">k<sub>i</sub> = <sup>2i + 1</sup>&frasl;<sub>2L<sub>x</sub></sub>&pi;</td><td rowspan="2">&alpha;<sub>i</sub> = 1</td><td rowspan="2">&forall;i &ge; 0</td></tr><tr><td><sup>&part;f</sup>&frasl;<sub>&part;x</sub>(x, y, z) = 0 for x = L<sub>x</sub></td><td>g&#773;(x) = cos(x)</td></tr><tr><td rowspan="2">ND</td><td><sup>&part;f</sup>&frasl;<sub>&part;x</sub>(x, y, z) = 0 for x = 0</td><td>g(x) = cos(x)</td><td rowspan="2">k<sub>i</sub> = <sup>2i + 1</sup>&frasl;<sub>2L<sub>x</sub></sub>&pi;</td><td rowspan="2">&alpha;<sub>i</sub> = 1</td><td rowspan="2">&forall;i &ge; 0</td></tr><tr><td>f(x, y, z) = 0 for x = L<sub>x</sub></td><td>g&#773;(x) = sin(x)</td></tr></tbody></table>

The 1-d domain of size $[0, L_x]$ is discretized uniformly with $N_x + 2$ points so the spatial scale associated with the uniform grid is $\Delta x = L_x/(N_x + 1)$ and we denote by $x_a = a\Delta x$ (for $a = 0, \dots, N+1$) the grid points. The grid points values of function $f$ are denoted by $f_a = f(x_a)$ (for $a = 0, \dots, N+1$). The discrete transform associated to each boundary conditions<fnref n="1" /> can be written under the generic form as follows

$$ \hat{F} = D_x(f), \quad (6) $$

where $f$ is the array of size $N+2$ containing the values of function $f$ at the grid points $x_a$ and $\hat{F}$ is the array of size $N+2$ containing the values of the associated discrete transform coefficients $\hat{F}_i$. The notation $D_x(\cdot)$ corresponds to the discrete transform operator (with $x = sI$ for the DST-I, $x = cI$ for the DCT-I, $x = sIII$ for the DST-III and $x = cIII$ for the DCT-III, depending on the boundary conditions, see footnote 3). The inverse discrete transform can be formally written as

$$ f = D_x^{-1}(\hat{F}), \quad (7) $$

where the notation $D_x^{-1}(\cdot)$ corresponds to the associated inverse discrete transform operator. The inverse transforms of type I (sine and cosine) discrete transforms are type I (sine and cosine) discrete transforms multiplied by a factor $2/(N+1)$, while the inverse transforms of type III (sine and cosine) discrete transforms are type II (sine and cosine) discrete transforms multiplied by a factor $2/(N+1)$.

The expressions of all discrete transforms is given in Table 2. Interestingly, one can notice that the coefficients $\hat{F}_i$ provide a numerical approximation of the coefficients $F_i$ (of the associated cosine or sine series), using a trapezoidal rule, through

$$ \hat{F}_i \approx \frac{N+1}{2} F_i, \quad \forall i \ge 0. \quad (8) $$

Finally, it is interesting to note that the four discrete transforms can be put under the following generic form

$$ \hat{F}_i = \sum_{a=0}^{N+1} \beta_a f_a g(k_i x_a), \quad \forall i \ge 0, \quad (9) $$

where $\beta_a$, $k_i$, $\alpha_i$ and $g(k_i x)$ depend on the boundary conditions. In addition, the inverse transforms can be written as

$$ f_a = \sum_{i=0}^{N+1} \gamma_i \hat{F}_i g(k_i x_i), \quad \forall a \ge 0. \quad (10) $$

## 2.2. The 3-d case

### 2.2.1. 3-d sine-cosine series

We consider a function $f$ defined in the prismatic domain $\Omega = [0, L_x] \times [0, L_y] \times [0, L_z]$ and verifying null Dirichlet boundary condition on faces $\partial\Omega_D$

$$ f(x, y, z) = 0, \quad \forall (x, y, z) \in \partial\Omega_D, \quad (11) $$

and null Neumann boundary condition on faces $\partial\Omega_N$

$$ \nabla f(x, y, z) \cdot n = 0, \quad \forall (x, y, z) \in \partial\Omega_N, \quad (12) $$

where $n$ is the outward normal to the boundary $\partial\Omega_N$ (with $\partial\Omega_N \cup \partial\Omega_D = \partial\Omega$). Function $f$ can be described by a generalized 3-d sine-cosine series as

$$ f(x, y, z) = \sum_{i=0}^{+\infty} \sum_{j=0}^{+\infty} \sum_{k=0}^{+\infty} \alpha_i^x \alpha_j^y \alpha_k^z F_{ijk} g^x(k_i^x x) g^y(k_j^y y) g^z(k_k^z z), \quad \forall (x, y, z) \in \Omega, \quad (13) $$

where the functions $g^x$, $g^y$, $g^z$, and the coefficients $\alpha_i^x$, $\alpha_j^y$, $\alpha_k^z$, $k_i^x$, $k_j^y$, $k_k^z$ depend on the type of boundary conditions on each couple of opposite faces<fnref n="2" /> (see Table 1), and $F_{ijk}$ are the sine-cosine series coefficients, given by

$$ F_{ijk} = \frac{8}{L_x L_y L_z} \int_{\Omega} f(x, y, z) g^x(k_i^x x) g^y(k_j^y y) g^z(k_k^z z) \, dxdydz, \quad \forall i, j, k \ge 0. \quad (14) $$

The partial derivatives of $f$ read

$$ \frac{\partial f}{\partial x}(x, y, z) = \sum_{i=0}^{+\infty} \sum_{j=0}^{+\infty} \sum_{k=0}^{+\infty} \epsilon_i^x \alpha_i^x \alpha_j^y \alpha_k^z F_{ijk} \overline{g^x(k_i^x x) g^y(k_j^y y) g^z(k_k^z z)}, \quad \forall (x, y, z) \in \Omega, \quad (15) $$

$$ \frac{\partial^2 f}{\partial x^2}(x, y, z) = \sum_{i=0}^{+\infty} \sum_{j=0}^{+\infty} \sum_{k=0}^{+\infty} -(\xi_i^x)^2 \alpha_i^x \alpha_j^y \alpha_k^z F_{ijk} g^x(k_i^x x) g^y(k_j^y y) g^z(k_k^z z), \quad \forall (x, y, z) \in \Omega, \quad (16) $$

where $\overline{g^x}$ is the dual function of $g^x$ (as defined in Table 1).

### 2.2.2. Discrete sine-cosine transforms in 3-d

The 3-d prismatic domain $\Omega$ is discretized uniformly with $(N_x+2) \times (N_y+2) \times (N_z+2)$ points so that the spatial scales are $\Delta x = L_x/(N_x+1)$, $\Delta y = L_y/(N_y+1)$ and $\Delta z = L_z/(N_z+1)$. The coordinates of the grid points are denoted by $x_a = a\Delta x$ (for $a = 0, \dots, N_x+1$), $y_b = b\Delta y$ (for $b = 0, \dots, N_y+1$) and $z_c = c\Delta z$ (for $c = 0, \dots, N_z+1$). The values of the function $f$ at the grid points are denoted by $f_{abc} = f(x_a, y_b, z_c)$ (for $a = 0, \dots, N_x+1$, $b = 0, \dots, N_y+1$ and $c = 0, \dots, N_z+1$).

Discrete sine-cosine transforms in 3-d are simply a composition of 1-d transforms in each spatial directions. The coefficients of the discrete sine-cosine transform $(\hat{F}_{ijk})$ of function $f$ are then given by (following the conventions of Table 2)

$$ \hat{F}_{ijk} = \sum_{a=0}^{N_x+1} \sum_{b=0}^{N_y+1} \sum_{c=0}^{N_z+1} \beta_a^x \beta_b^y \beta_c^z f_{abc} g^x(k_i^x x_a) g^y(k_j^y y_b) g^z(k_k^z z_c), \quad (17) $$

where the functions $g^x$, $g^y$, $g^z$, and the coefficients $k_i^x$, $k_j^y$, $k_k^z$, $\beta_i^x$, $\beta_j^y$ and $\beta_k^z$ depend on the type of boundary conditions. Let us denote by $\hat{F}$ the

<footnote n="1">¹ DST-I encodes Dirichlet–Dirichlet, DCT-I encodes Neumann–Neumann, DST-III encodes Dirichlet–Neumann and DCT-III encodes Neumann–Dirichlet.</footnote>

<footnote n="2">² There is a total of $2^6 = 64$ different sine-cosine series in the 3-d case (corresponding to the 2 types of boundary conditions on the 6 faces).</footnote>
<table><caption>**Table 2**
Computational details of the Discrete Sine-Cosine Transforms.</caption><thead><tr><th>BC</th><th>Transform</th><th>Î&#x0302;<sub>i</sub> (i = 0, ..., N + 1)</th><th>f<sub>a</sub> (a = 0, ..., N + 1)</th></tr></thead><tbody><tr><td rowspan="3">DD</td><td>DST-I</td><td>Î&#x0302;<sub>i</sub> = &sum;<sub>a=0</sub><sup>N+1</sup> &beta;<sub>a</sub>f<sub>a</sub> sin(&pi;a/N + 1)</td><td>f<sub>a</sub> = 2/(N + 1) &sum;<sub>i=0</sub><sup>N+1</sup> &gamma;<sub>i</sub>Î&#x0302;<sub>i</sub> sin(&pi;a/N + 1)</td></tr><tr><td>Î&#x0302; = D<sub>cl</sub>(f)</td><td>&beta;<sub>0</sub> = &beta;<sub>N+1</sub> = 1/2</td><td>&gamma;<sub>0</sub> = &gamma;<sub>N+1</sub> = 1/2</td></tr><tr><td>f = D<sup>-1</sup><sub>sl</sub>(Î&#x0302;)</td><td>&beta;<sub>a</sub> = 1 for a = 1, ..., N</td><td>&gamma;<sub>i</sub> = 1 for i = 1, ..., N</td></tr><tr><td rowspan="3">NN</td><td>DCT-I</td><td>Î&#x0302;<sub>i</sub> = &sum;<sub>a=0</sub><sup>N+1</sup> &beta;<sub>a</sub>f<sub>a</sub> cos(&pi;a/N + 1)</td><td>f<sub>a</sub> = 2/(N + 1) &sum;<sub>i=0</sub><sup>N+1</sup> &gamma;<sub>i</sub>Î&#x0302;<sub>i</sub> cos(&pi;a/N + 1)</td></tr><tr><td>Î&#x0302; = D<sub>cl</sub>(f)</td><td>&beta;<sub>0</sub> = &beta;<sub>N+1</sub> = 1/2</td><td>&gamma;<sub>0</sub> = &gamma;<sub>N+1</sub> = 1/2</td></tr><tr><td>f = D<sup>-1</sup><sub>cl</sub>(Î&#x0302;)</td><td>&beta;<sub>a</sub> = 1 for a = 1, ..., N</td><td>&gamma;<sub>i</sub> = 1 for i = 1, ..., N</td></tr><tr><td rowspan="3">DN</td><td>DST-III</td><td>Î&#x0302;<sub>i</sub> = &sum;<sub>a=0</sub><sup>N+1</sup> &beta;<sub>a</sub>f<sub>a</sub> sin(&pi;a(i + 1/2)/(N + 1))</td><td>f<sub>a</sub> = 2/(N + 1) &sum;<sub>i=0</sub><sup>N+1</sup> &gamma;<sub>i</sub>Î&#x0302;<sub>i</sub> sin(&pi;a(i + 1/2)/(N + 1))</td></tr><tr><td>Î&#x0302; = D<sub>slsl</sub>(f)</td><td>&beta;<sub>0</sub> = &beta;<sub>N+1</sub> = 1/2</td><td>&gamma;<sub>i</sub> = 1 for i = 0, ..., N</td></tr><tr><td>f = D<sup>-1</sup><sub>slsl</sub>(Î&#x0302;)</td><td>&beta;<sub>a</sub> = 1 for a = 1, ..., N</td><td>&gamma;<sub>N+1</sub> = 0</td></tr><tr><td rowspan="3">ND</td><td>DCT-III</td><td>Î&#x0302;<sub>i</sub> = &sum;<sub>a=0</sub><sup>N+1</sup> &beta;<sub>a</sub>f<sub>a</sub> cos(&pi;a(i + 1/2)/(N + 1))</td><td>f<sub>a</sub> = 2/(N + 1) &sum;<sub>i=0</sub><sup>N+1</sup> &gamma;<sub>i</sub>Î&#x0302;<sub>i</sub> cos(&pi;a(i + 1/2)/(N + 1))</td></tr><tr><td>Î&#x0302; = D<sub>slsl</sub>(f)</td><td>&beta;<sub>0</sub> = &beta;<sub>N+1</sub> = 1/2</td><td>&gamma;<sub>i</sub> = 1 for i = 0, ..., N</td></tr><tr><td>f = D<sup>-1</sup><sub>slsl</sub>(Î&#x0302;)</td><td>&beta;<sub>a</sub> = 1 for a = 1, ..., N</td><td>&gamma;<sub>N+1</sub> = 0</td></tr></tbody></table>

array of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the discrete sine-cosine coefficients and by $f$ the array of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of function $f$ at the grid points. The discrete transform (17) is then written as

$$ \hat{F} = D_{xyz}(f), \qquad (18) $$

where $x \in \{sI, cI, sIII, cIII\}, y \in \{sI, cI, sIII, cIII\}$ and $z \in \{sI, cI, sIII, cIII\}$, depending on the boundary conditions applied on the x-faces, y-faces and z-faces, respectively.

The associated inverse transform reads

$$ f_{abc} = \frac{8}{(N_x+1)(N_y+1)(N_z+1)} \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \gamma_i^x \gamma_j^y \gamma_k^z \hat{F}_{ijk} g^x(k_a^x x_i) \times g^y(k_j^x y_b)g^z(k_k^z z_c), \qquad (19) $$

and it can be written as

$$ f = D_{xyz}^{-1}(\hat{F}). \qquad (20) $$

A useful notation is finally introduced to account for dual mixed sine-cosine transforms. Assuming that the transform $D_{xyz}(\cdot)$ corresponds to Eq. (17), then we denote by

$$ \hat{F} = D_{\bar{z}yz}(f), \qquad (21) $$

the transform defined by

$$ \hat{F}_{ijk} = \sum_{a=0}^{N_x+1} \sum_{b=0}^{N_y+1} \sum_{c=0}^{N_z+1} \beta_a^x \beta_b^y \beta_c^z f_{abc} \bar{g}^x(k_i^x x_a) g^y(k_j^x y_b) g^z(k_k^z z_c). \qquad (22) $$

The notation $D_{\bar{z}yz}$ defines the dual transform of $D_{\bar{z}yz}$ in the x-direction. For example, if $D_{xyz}$ corresponds to the transform $D_{sIsIsI}$, then $D_{\bar{z}yz}$ corresponds to $D_{cIsIsI}$.

## 3. A fast numerical method for transient conductivity problems in heterogeneous media subjected to mixed Dirichlet/Neumann boundary conditions

### 3.1. Equations of conductivity

We are concerned in solving the equations of (transient) thermal conductivity in a heterogeneous finite medium. The finite cell is a cubic domain denoted by $\Omega = [0, L_x] \times [0, L_y] \times [0, L_z]$ in 3-d and a square denoted by $\Omega = [0, L_x] \times [0, L_y]$ in 2-d. Tensorial components refer to a system of Cartesian coordinates $(e_x; e_y; e_z)$ in 3-d (and $(e_x; e_y)$ in 2-d).

The problem of transient heat conduction (without source terms) consists in the computation of the heat flow $q(x, t)$, and temperature

$T(x, t)$, at each point $x$ in $\Omega$ and instants $t > 0$, for a given heterogeneous conductivity field $c(x)$, density $\rho(x)$ and specific heat capacity $c_p(x)$, solution of the local problem

$$ \left\{ \begin{array}{ll} \forall(x, t) \in \Omega \times [0, t_{\max}], & \text{div } q(x, t) = -\rho(x)c_p(x) \frac{\partial T}{\partial t}(x, t) \\ \forall(x, t) \in \Omega \times [0, t_{\max}], & q(x) = -c(x) : \nabla T(x, t). \end{array} \right. \qquad (23) $$

In the case of isotropic conductivity, the local second-order conductivity tensor $c(x)$ reduces to (in the case of a 3-d medium)

$$ c(x) = c(x)\mathbf{I}_3 \qquad (24) $$

with $c(x)$ the local (scalar) conductivity field and $\mathbf{I}_3$ the second-order identity tensor.

Mixed Dirichlet/Neumann boundary conditions are considered and are given by

$$ \left\{ \begin{array}{ll} \forall \mathbf{x} \in \partial\Omega_D, & T(\mathbf{x}, t) = T_{\partial\Omega_D}(\mathbf{x}, t) & \text{(Dirichlet)} \\ \forall \mathbf{x} \in \partial\Omega_N, & \nabla T(\mathbf{x}, t) \cdot \mathbf{n} = f_{\partial\Omega_N}(\mathbf{x}, t) & \text{(Neumann)}, \end{array} \right. \qquad (25) $$

where $T_{\partial\Omega_D}(\mathbf{x})$ is the prescribed temperature, on the boundary $\partial\Omega_D$, corresponding to the Dirichlet condition, $\mathbf{n}$ is the outward normal and $f_{\partial\Omega_N}(\mathbf{x})$ is the prescribed value for the normal derivative of $T$, on the boundary $\partial\Omega_N$, corresponding to the Neumann condition (and can be related to the normal heat flux). It should be noted that the boundary conditions are applied on entire faces (with of course $\partial\Omega_D \cup \partial\Omega_N = \partial\Omega$).

In addition, we consider the following initial condition:

$$ \forall \mathbf{x} \in \Omega, \quad T(\mathbf{x}, 0) = T_{\text{init}}(\mathbf{x}), \qquad (26) $$

where $T_{\text{init}}(\mathbf{x})$ is the initial temperature field at $t = 0$.

In the particular case of a homogeneous and isotropic medium, i.e. the scalar conductivity field $c$, the density $\rho$ and the specific heat capacity $c_p$ are homogeneous, the equations to be solved reduce to the standard heat equation

$$ \rho c_p \frac{\partial T}{\partial t}(\mathbf{x}, t) = c \Delta T(\mathbf{x}, t), \qquad (27) $$

where $\Delta$ is the Laplacian operator.

### 3.2. Principle of resolution

The problem is solved incrementally by introducing a finite number of time steps $t^n$ defined by

$$ t^n = n\Delta t, \qquad (28) $$

where $\Delta t$ is the time step size. The quantities at time $t^n$ are denoted by

$$ \mathbf{q}^n(\mathbf{x}) = \mathbf{q}(\mathbf{x}, t^n), \quad T^n(\mathbf{x}) = T(\mathbf{x}, t^n). \qquad (29) $$
The resolution will consist in finding the thermal state $T^{n+1}(\mathbf{x})$ at time $t^{n+1}$, knowing the previous thermal state $T^n(\mathbf{x})$ at time $t^n$.

First, we consider a theta-method to discretize problem (23) in time. The thermal problem therefore reads

$$ \left\{ \begin{array}{ll} \forall \mathbf{x} \in \Omega, & \theta \operatorname{div} \mathbf{q}^{n+1}(\mathbf{x}) + (1-\theta)\operatorname{div} \mathbf{q}^n(\mathbf{x}) = -\rho(\mathbf{x})c_p(\mathbf{x}) \frac{T^{n+1}(\mathbf{x}) - T^n(\mathbf{x})}{\Delta t} \\ \forall \mathbf{x} \in \Omega, & \mathbf{q}^{n+1}(\mathbf{x}) = -\mathbf{c}(\mathbf{x}) : \nabla T^{n+1}(\mathbf{x}) \end{array} \right. \qquad (30) $$

subjected to the boundary conditions

$$ \left\{ \begin{array}{ll} \forall \mathbf{x} \in \partial\Omega_D, & T^{n+1}(\mathbf{x}) = T_{\partial\Omega_D}(\mathbf{x}, t^{n+1}) \\ \forall \mathbf{x} \in \partial\Omega_N, & \nabla T^{n+1}(\mathbf{x}) \cdot \mathbf{n} = f_{\partial\Omega_N}(\mathbf{x}, t^{n+1}). \end{array} \right. \qquad (31) $$

The values $\theta = 0$, $\theta = 0.5$ and $\theta = 1$ respectively correspond to a forward Euler explicit scheme, Crank–Nicolson scheme and backward Euler implicit scheme (which will be simply called implicit scheme hereafter).

In order to solve (30), we extend (Paux et al., 2025a)'s approach (see also Gélébart (2024)) to the transient case using the following procedure:

1. Write the problem (30) as an auxiliary problem involving a reference homogeneous material and a polarization term (as done by Moulinec and Suquet (1998) in the periodic setting);

2. Split the temperature field in two contributions, a known term verifying the boundary conditions (31) and an unknown fluctuation term, defined such that it is null on the boundary $\partial\Omega_D$ and it has normal derivatives null on the boundary $\partial\Omega_N$, to be determined;

3. Write the weak (Galerkin) formulation by (i) describing the fluctuation term by a partial 3-d mixed sine–cosine series (depending on the type of boundary conditions) and (ii) using trial functions based on the same 3-d mixed sine–cosine series;

4. Solve approximately the problem using discrete sine and cosine transforms (to approximate the integrals defined in the weak formulation);

5. Find the polarization term solution of the problem (30) using an iterative scheme (as done by Moulinec and Suquet (1998) in the periodic setting with a fixed-point method).

The main difference between the stationary and the transient cases is the heterogeneity of material properties in Eq. (30)$_1$ ($\rho(\mathbf{x})$ and $c_p(\mathbf{x})$; this will require the introduction of two polarization terms in the auxiliary problem.

### 3.3. The auxiliary problem

Following Moulinec and Suquet (1998) (see also Paux et al. (2025a) and Gélébart (2024)), the local conductivity problem (subjected to the boundary conditions (31)) is written alternatively in the Lippmann–Schwinger form by introducing a homogeneous isotropic reference material

$$ \left\{ \begin{array}{ll} \forall \mathbf{x} \in \Omega, & \theta \operatorname{div} \mathbf{q}^{n+1}(\mathbf{x}) + \rho_0 c_{p0} \frac{T^{n+1}(\mathbf{x})}{\Delta t} + \tau_1(\mathbf{x}) \\ & = \rho(\mathbf{x}) c_p(\mathbf{x}) \frac{T^n(\mathbf{x})}{\Delta t} + (\theta - 1) \operatorname{div} \mathbf{q}^n(\mathbf{x}) \\ \forall \mathbf{x} \in \Omega, & \mathbf{q}^{n+1}(\mathbf{x}) = -c_0 \nabla T^{n+1}(\mathbf{x}) + \tau_2(\mathbf{x}) \end{array} \right. \qquad (32) $$

where $c_0$, $\rho_0$ and $c_{p0}$ are respectively the conductivity, density and specific heat capacity of some homogeneous isotropic reference material and

$$ \tau_1(\mathbf{x}) = (\rho(\mathbf{x})c_p(\mathbf{x}) - \rho_0 c_{p0}) \frac{T^{n+1}(\mathbf{x})}{\Delta t}, \quad \tau_2(\mathbf{x}) = -(\mathbf{c}(\mathbf{x}) - c_0 \mathbf{I}_3) : \nabla T^{n+1}(\mathbf{x}), \qquad (33) $$

are the so-called polarization tensors. The polarization tensors are introduced in order to write an auxiliary problem that leads to a diagonal system to be solved after spatial discretization.

Assuming momentarily that the polarization tensors $\tau_1(\mathbf{x})$ and $\tau_2(\mathbf{x})$ are known, the problem reduces to the calculation of the current temperature $T^{n+1}(\mathbf{x})$ solution of

$$ -\theta c_0 \nabla^2 T^{n+1}(\mathbf{x}) + \rho_0 c_{p0} \frac{T^{n+1}(\mathbf{x})}{\Delta t} = \rho(\mathbf{x}) c_p(\mathbf{x}) \frac{T^n(\mathbf{x})}{\Delta t} - \tau_1(\mathbf{x}) + \operatorname{div} \left( (\theta - 1) \mathbf{q}^n(\mathbf{x}) - \theta \tau_2(\mathbf{x}) \right), \qquad (34) $$

where the notation $\nabla^2$ corresponds to the Laplacian operator.

Following the ideas of Gélébart (2024) and Paux et al. (2025a), the boundary conditions are taken into account by splitting the temperature field $T^{n+1}(\mathbf{x})$ solution of the problem in two terms

$$ T^{n+1}(\mathbf{x}) = T_0^{n+1}(\mathbf{x}) + \tilde{T}^{n+1}(\mathbf{x}), \qquad (35) $$

where $T_0^{n+1}(\mathbf{x})$ is a known analytical $c^1$ function verifying the applied boundary conditions (25) and $\tilde{T}^{n+1}(\mathbf{x})$ is an unknown fluctuation function verifying the properties

$$ \left\{ \begin{array}{ll} \forall \mathbf{x} \in \partial\Omega_D, & \tilde{T}^{n+1}(\mathbf{x}) = 0, \\ \forall \mathbf{x} \in \partial\Omega_N, & \nabla \tilde{T}^{n+1}(\mathbf{x}) \cdot \mathbf{n} = 0. \end{array} \right. \qquad (36) $$

The known function $T_0^{n+1}(\mathbf{x})$ can be easily reconstructed from the knowledge of the prescribed temperature $T_{\partial\Omega_D}$ and heat fluxes $f_{\partial\Omega_N}$ at $t^{n+1}$ by interpolation within the cell (see e.g. Paux et al. (2025a)).

Assuming that $\tau_1(\mathbf{x})$, $\tau_2(\mathbf{x})$, $T_0^{n+1}(\mathbf{x})$, $T^n(\mathbf{x})$ and $\mathbf{q}^n(\mathbf{x})$ are known, the problem reduces to the determination of the fluctuation $\tilde{T}^{n+1}(\mathbf{x})$ (verifying the null boundary conditions (36)) solution of

$$ -\theta c_0 \nabla^2 \tilde{T}^{n+1}(\mathbf{x}) + \frac{\rho_0 c_{p0}}{\Delta t} \tilde{T}^{n+1}(\mathbf{x}) = A(\mathbf{x}) + \operatorname{div} \mathbf{B}(\mathbf{x}), \qquad (37) $$

where $A(\mathbf{x})$ and $\mathbf{B}(\mathbf{x})$ are (known) terms given by

$$ \begin{aligned} A(\mathbf{x}) &= \rho(\mathbf{x})c_p(\mathbf{x}) \frac{T^n(\mathbf{x})}{\Delta t} - \tau_1(\mathbf{x}) + \theta c_0 \nabla^2 T_0^{n+1}(\mathbf{x}) - \frac{\rho_0 c_{p0}}{\Delta t} T_0^{n+1}(\mathbf{x}), \\ \mathbf{B}(\mathbf{x}) &= (\theta - 1)\mathbf{q}^n(\mathbf{x}) - \theta \tau_2(\mathbf{x}). \end{aligned} \qquad (38) $$

The determination of $\tau_1(\mathbf{x})$ and $\tau_2(\mathbf{x})$ solution of the problem is then done using an iterative scheme.

### 3.4. Weak formulation

The 3-d prismatic domain is discretized uniformly with $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ points and the spatial scales associated with the uniform grid are $\Delta x = L_x/(N_x + 1)$, $\Delta y = L_y/(N_y + 1)$ and $\Delta z = L_z/(N_z + 1)$. The coordinates of the grid points are denoted by $x_a = a\Delta x$ (for $a = 0, ..., N_x + 1$), $y_b = b\Delta y$ (for $b = 0, ..., N_y + 1$) and $z_c = c\Delta z$ (for $c = 0, ..., N_z + 1$).

The fluctuation term $\tilde{T}^{n+1}(\mathbf{x})$ is defined by construction so as to verify Eq. (36). Therefore, it will be described by a partial 3-d mixed sine–cosine series, chosen to ensure “null” boundary conditions. We consider that $\tilde{T}^{n+1}(\mathbf{x})$ is of the form

$$ \tilde{T}^{n+1}(\mathbf{x}) = \sum_{i=0}^{N_x} \sum_{j=0}^{N_y} \sum_{k=0}^{N_z} \alpha_i^x \alpha_j^y \alpha_k^z \tilde{T}_{ijk} g^x(k_i^x x) g^y(k_j^y y) g^z(k_k^z z), \qquad (39) $$

where $\tilde{T}_{ijk}$ are the 3-d sine–cosine series coefficients of function $\tilde{T}^{n+1}$, which constitute the unknowns of our problem. The choice of the functions $g^x$, $g^y$, $g^z$ and coefficients $k_i^x$, $k_j^y$, $k_k^z$, $\alpha_i^x$, $\alpha_j^y$, $\alpha_k^z$ then simply depends on the type of boundary conditions (see Section 2). As explained by Paux et al. (2025a), the number of “modes” for the sine–cosine series is equal to $(N_x + 1) \times (N_y + 1) \times (N_z + 1)$ and does not correspond to the number of grid points, because the last coefficients $\tilde{T}_{ijk}$ may be poorly estimated for the high frequencies (i.e. for $i = N_x + 1$ or $j = N_y + 1$ or $k = N_z + 1$). For homogeneity purposes, Eq. (39) can be alternatively
written

$$ \tilde{T}^{n+1}(\mathbf{x}) = \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \alpha_i^x \alpha_j^y \alpha_k^z \tilde{T}_{ijk} g^x(k_i^x) g^y(k_j^y) g^z(k_k^z), \quad (40) $$

where we enforce that

$$ \tilde{T}_{N_x+1,jk} = \tilde{T}_{iN_y+1,k} = \tilde{T}_{ijN_z+1} = 0. \quad (41) $$

The weak formulation of Lippmann-Schwinger equation (37) reads

$$ \int_{\Omega} \left( -\theta c_0 \nabla^2 \tilde{T}^{n+1}(\mathbf{x}) + \frac{\rho_0 c_p \theta}{\Delta t} \tilde{T}^{n+1}(\mathbf{x}) \right) v(\mathbf{x}) d\Omega \\ = \int_{\Omega} (A(\mathbf{x}) + \text{div } \mathbf{B}(\mathbf{x})) v(\mathbf{x}) d\Omega, \quad (42) $$

where $v(\mathbf{x})$ is a testing differentiable function. Following (39), these trial functions are the interpolation basis of $\tilde{T}(\mathbf{x})$, i.e.,

$$ v_{lmn}(\mathbf{x}) = g^x(k_i^x)g^y(k_j^y)g^z(k_n^z), \quad (43) $$

with $l = 0, \dots, N_x$, $m = 0, \dots, N_y$ and $n = 0, \dots, N_z$. The integrals involved in the weak formulation (42), that need to be calculated, are denoted by

$$ \left\{ \begin{aligned} I_{lmn} &= \int_{\Omega} \left( -\theta c_0 \nabla^2 \tilde{T}^{n+1}(\mathbf{x}) + \frac{\rho_0 c_p \theta}{\Delta t} \tilde{T}^{n+1}(\mathbf{x}) \right) v_{lmn}(\mathbf{x}) d\Omega \\ J_{lmn} &= \int_{\Omega} \operatorname{div} \mathbf{B}(\mathbf{x}) v_{lmn}(\mathbf{x}) d\Omega \\ K_{lmn} &= \int_{\Omega} A(\mathbf{x}) v_{lmn}(\mathbf{x}) d\Omega \end{aligned} \right. \quad (44) $$

After calculation of some elementary integrals (see Appendix for detailed calculations), one obtains

$$ \left\{ \begin{aligned} I_{lmn} &= \frac{L_x L_y L_z}{8} \left( \theta c_0 \left( (\xi_l^x)^2 + (\xi_m^y)^2 + (\xi_n^z)^2 \right) + \frac{\rho_0 c_p \theta}{\Delta t} \right) \tilde{T}_{lmn} \\ J_{lmn} &= \frac{-L_x L_y L_z}{(N_x + 1)(N_y + 1)(N_z + 1)} \left[ \xi_l^x (D_{\bar{x}\bar{y}z}(\tilde{\mathbf{B}}_x))_{lmn} + \xi_m^y (D_{\bar{x}\bar{y}z}(\tilde{\mathbf{B}}_y))_{lmn} + \xi_n^z (D_{\bar{x}\bar{y}z}(\tilde{\mathbf{B}}_z))_{lmn} + S_{lmn} \right] \\ K_{lmn} &= \frac{L_x L_y L_z}{(N_x + 1)(N_y + 1)(N_z + 1)} (D_{\bar{x}\bar{y}z}(\mathbf{A}))_{lmn} \end{aligned} \right. \quad (45) $$

where $\mathbf{A}, \tilde{\mathbf{B}}_x, \tilde{\mathbf{B}}_y$, and $\tilde{\mathbf{B}}_z$ respectively denote the arrays of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of functions $\mathbf{A}, \tilde{\mathbf{B}}_x, \tilde{\mathbf{B}}_y$, and $\tilde{\mathbf{B}}_z$ at the grid points, where the functions $\tilde{\mathbf{B}}_x, \tilde{\mathbf{B}}_y, \tilde{\mathbf{B}}_z$ are given by Eq. (A.10). The quantity $S_{lmn}$ is given by (A.15).

Let us denote by $\tilde{T}$ the array of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of the sine-cosine series coefficients $\tilde{T}_{ijk}$ of the unknown fluctuation $\tilde{T}^{n+1}$ defined by Eq. (39). Combining Eqs. (42) and (45), the coefficients $\tilde{T}$ are given by

$$ \tilde{T} = C \odot \left[ D_{xyz}(\mathbf{A}) - \frac{1}{c_0} \left( \xi^x \odot D_{\bar{x}\bar{y}z}(\tilde{\mathbf{B}}_x) + \xi^y \odot D_{\bar{x}\bar{y}z}(\tilde{\mathbf{B}}_y) + \xi^z \odot D_{\bar{x}\bar{y}z}(\tilde{\mathbf{B}}_z) + S \right) \right], \quad (46) $$

where $C$ is the $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ array given by

$$ (\mathbf{C})_{ijk} = \frac{1}{\theta c_0 \left( (\xi_l^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2 \right) + \frac{\rho_0 c_p \theta}{\Delta t}} \times \frac{8}{(N_x + 1)(N_y + 1)(N_z + 1)}, \quad (47) $$

$\xi^x, \xi^y, \xi^z$ are the $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ arrays given by

$$ (\xi^x)_{ijk} = \xi_i^x, \quad (\xi^y)_{ijk} = \xi_j^y, \quad (\xi^z)_{ijk} = \xi_k^z, \quad (48) $$

$\mathbf{S}$ is the $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ array containing the values $S_{lmn}$ and $\odot$ denotes Hadamard product (pointwise product) given by

$$ (\mathbf{a} \odot \mathbf{b})_{ijk} = (\mathbf{a})_{ijk}(\mathbf{b})_{ijk} \quad (\text{no summation on the indices}). \quad (49) $$

It must be noted that Eq. (46) allows the determination of $\tilde{T}_{lmn}$ for all frequencies except for $l = N_x + 1, m = N_y + 1$ and $n = N_z + 1$ where Eq. (41) holds. Once the array $\tilde{T}$ is known, the fluctuation term $\tilde{T}^{n+1}(\mathbf{x})$ is known using Eq. (39).

### 3.5. The iterative scheme

The last step is to determine the polarization tensors $\tau_1$ and $\tau_2$ solution of the problem (23) using an iterative scheme. We consider, following Moulinec and Suquet (1998), a fixed-point iterative scheme in order to calculate the thermal state $T^{n+1}(\mathbf{x})$ at time $t^{n+1}$, knowing the previous thermal state $T^n(\mathbf{x})$ at time $t^n$:

$$ \left\{ \begin{array}{ll} \text{Initialization} & \text{Thermal state at time } t^n \text{ is known: } T^n(\mathbf{x}) \\ & (T^{n+1}(\mathbf{x}))^0 = T^n(\mathbf{x}) \\ & (\mathbf{q}^{n+1}(\mathbf{x}))^0 = -c(\mathbf{x}) : \nabla (T^{n+1}(\mathbf{x}))^0 \\ & (T^{n+1}(\mathbf{x}))^i \text{ and } (\mathbf{q}^{n+1}(\mathbf{x}))^i \text{ being known} \\ \text{Iterate } i+1 & (\tau_1(\mathbf{x}))^i = \frac{\rho(\mathbf{x})c_p(\mathbf{x}) - \rho_0 c_p \theta}{\Delta t} (T^{n+1}(\mathbf{x}))^i \\ & (\tau_2(\mathbf{x}))^i = -(c(\mathbf{x}) - c_0 \mathbf{I}_3) : \nabla (T^{n+1}(\mathbf{x}))^i, \\ & (A(\mathbf{x}))^i = \rho(\mathbf{x})c_p(\mathbf{x}) \frac{T^n(\mathbf{x})}{\Delta t} - (\tau_1(\mathbf{x}))^i + \theta c_0 \nabla^2 T_0^{n+1}(\mathbf{x}) \\ & - \frac{\rho_0 c_p \theta}{\Delta t} T_0^{n+1}(\mathbf{x}) \\ & (\mathbf{B}(\mathbf{x}))^i = (\theta - 1) \mathbf{q}^n(\mathbf{x}) - \theta (\tau_2(\mathbf{x}))^i \\ & (b) \quad \tilde{T}^{i+1} = C \odot \left[ D_{xyz}(\mathbf{A}^i) - \frac{1}{c_0} \left( \xi^x \odot D_{\bar{x}\bar{y}z}(\tilde{\mathbf{B}}_x^i) + \xi^y \odot D_{\bar{x}\bar{y}z}(\tilde{\mathbf{B}}_y^i) + \xi^z \odot D_{\bar{x}\bar{y}z}(\tilde{\mathbf{B}}_z^i) + S^i \right) \right] \\ & (c) \quad \left( \tilde{T}^{n+1}(\mathbf{x}) \right)^{i+1} = \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \alpha_i^x \alpha_j^y \alpha_k^z \tilde{T}_{ijk}^{i+1} g^x(k_i^x) \\ & \times g^y(k_j^y) g^z(k_k^z) \\ & (d) \quad (\mathbf{q}^{n+1}(\mathbf{x}))^{i+1} = -c(\mathbf{x}) : \left( \nabla T_0^{n+1}(\mathbf{x}) + \nabla (\tilde{T}^{n+1}(\mathbf{x}))^{i+1} \right) \\ & (e) \quad \text{Convergence test.} \end{array} \right. \quad (50) $$

The convergence test consists in verifying the conservation Eq. (30)$_1$ in the weak sense. Let us introduce the error function $\epsilon^{n+1}(\mathbf{x})$

$$ \epsilon^{n+1}(\mathbf{x}) = \theta \operatorname{div} \mathbf{q}^{n+1}(\mathbf{x}) + (1 - \theta) \operatorname{div} \mathbf{q}^n(\mathbf{x}) + \rho(\mathbf{x})c_p(\mathbf{x}) \frac{T^{n+1}(\mathbf{x}) - T^n(\mathbf{x})}{\Delta t}. \quad (51) $$

Following Parseval's identity, the $L^2$ norm of the error function reads

$$ \| \epsilon^{n+1}(\mathbf{x}) \|_{L^2}^2 = \frac{1}{L_x L_y L_z} \int_{\Omega} (\epsilon^{n+1}(\mathbf{x}))^2 d\Omega = \sum_{i=0}^{N_x} \sum_{j=0}^{N_y} \sum_{k=0}^{N_z} (E_{ijk}^{n+1})^2 \quad (52) $$

where $E_{ijk}^{n+1}$ is given by

$$ E_{ijk}^{n+1} = \frac{8}{L_x L_y L_z} \int_{\Omega} \epsilon^{n+1}(\mathbf{x}) v_{ijk}(\mathbf{x}) d\Omega. \quad (53) $$

These integrals are calculated using discrete sine-cosine transforms as described in Appendix; convergence is considered achieved when

$$ \Delta x \sqrt{\sum_{i=0}^{N_x} \sum_{j=0}^{N_y} \sum_{k=0}^{N_z} (E_{ijk}^{n+1})^2} \le \text{tol.} \quad (54) $$

It should be noted that, introducing the iterative expressions of $(\tau_1)^i$, $(\tau_2)^i$, $A^i$ and $\mathbf{B}^i$ (Eq. (50)a) in the expression of $\tilde{T}^{i+1}$ (Eq. (50)b), the expression $\tilde{T}^{i+1}$ can be expressed thanks to the error functions (51) at iteration $i$. Formally, it can be written as

$$ \tilde{T}^{i+1} = \tilde{T}^i - C \odot E^i, \quad (55) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1. Assessment of the numerical method in a homogeneous reference problem. (a) Distribution of the temperature in half of the domain and (b) Evolution of the temperature at the center of the cell.</figcaption></figure>

with $E^i$ the array containing the coefficients of $E_{ijk}$ of the error function at iteration $i$. Introducing the arrays $dq^n$ and $dq^{n+1}$ corresponding to the divergences of $q^n$ and $q^{n+1}$, respectively, and the array $dT$ corresponding to the term $\rho(x)c_p(x)^{T^{n+1}(x)-T^n(x)}$ in (51), and following the computation method described in Section 3.4, it leads to the following modified iterative scheme,

$$
\left\{
\begin{array}{ll}
\text{Initialization} & \text{Thermal state at time } t^n \text{ is known: } T^n(x), \\
& E^0 = 0, \quad \tilde{T}^0 = 0 \\
\text{Iterate } i+1 & dq^n = \xi^x \odot D_{\bar{x}\bar{y}z}(\tilde{q}_x^n) + \xi^y \odot D_{\bar{x}\bar{y}z}(\tilde{q}_y^n) + \xi^z \odot D_{\bar{x}\bar{y}\bar{z}}(\tilde{q}_z^n) \\
& \tilde{T}^i \text{ and } E^i \text{ being known} \\
\text{(a)} & \tilde{T}^{i+1} = \tilde{T}^i - C \odot E^i \\
\text{(b)} & (\tilde{T}^{n+1}(x))^{i+1} = \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \alpha_i^x \alpha_j^y \alpha_k^z \tilde{T}^{i+1} g^x(k_i^x x) \\
& \times g^y(k_j^y y) g^z(k_z^z z) \\
\text{(c)} & (dT(x))^{i+1} = \rho(x)c_p(x) \frac{(T^{n+1}(x))^{i+1}-T^n(x)}{\Delta t} \\
\text{(d)} & (q^{n+1}(x))^{i+1} = -c(x) : \left( \nabla T_0^{n+1}(x) + \nabla (\tilde{T}^{n+1}(x))^{i+1} \right) \\
\text{(e)} & (dq^{n+1})^{i+1} = \xi^x \odot D_{\bar{x}\bar{y}z}((\tilde{q}_x^{n+1})^{i+1}) \\
& + \xi^y \odot D_{\bar{x}\bar{y}z}((\tilde{q}_y^{n+1})^{i+1}) + \xi^z \odot D_{\bar{x}\bar{y}\bar{z}}((\tilde{q}_z^{n+1})^{i+1}) \\
\text{(f)} & E^{i+1} = \theta dq^{n+1})^{i+1} + (1-\theta) dq^n + D_{\bar{x}\bar{y}z}(dT^{i+1}) \\
\text{(g)} & \text{Convergence test}
\end{array}
\right.
\quad (56)
$$

This modified iterative scheme avoids redundancies in the computations and has the advantage to ensure consistency between the stabilization of the numerical fixed point and the convergence test. Therefore, this scheme will be used for the forthcoming applications.

## 4. Assessment of the method in a homogeneous case

The numerical method developed in Section 3 is assessed using a 1d problem with homogeneous properties for which an analytical solution is known. In the domain $[0, L]$, we consider the following diffusion equation

$$ \frac{\partial T^2}{\partial x^2}(x, t) = \frac{\rho c_p}{c} \frac{\partial T}{\partial t}(x, t) \quad (57) $$

with null Dirichlet conditions

$$ T(0, t) = T(L, t) = 0 \quad (58) $$

and initial condition

$$ T(x, t = 0) = T_0(x). \quad (59) $$

The analytical solution of this equation can be easily calculated using separation of variables and reads

$$ T(x, t) = \sum_{k=1}^{\infty} b_k \exp \left( - \left( \frac{\rho c_p}{c} \frac{k\pi}{L} \right)^2 t \right) \sin \left( \frac{k\pi x}{L} \right), \quad (60) $$

with $b_k$ given by

$$ b_k = \frac{2}{L} \int_0^L T_0(x) \sin \left( \frac{k\pi x}{L} \right) dx. \quad (61) $$

In order to assess the method, we consider the following test case

$$ T(x, t = 0) = T_0(x) = 4T_{\max} \frac{xL - x^2}{L^2}, \quad (62) $$

for which the analytical solution is given by Eq. (60), where the coefficients $b_k$ read

$$ b_k = \frac{16LT_{\max}}{k^3 \pi^3} (1 - (-1)^k). \quad (63) $$

The following parameters are taken for the assessment:

* Initial maximal temperature: $T_{\max} = 1$ °C;

* Material properties: $c = 1$ W m$^{-1}$ K$^{-1}$, $c_p = 1$ J kg$^{-1}$ K$^{-1}$, $\rho_1 = 1$ kg m$^{-3}$.

* Domain size $L = 1$ m.

For the numerical simulations, we consider 129 grid points and two time steps ($\Delta t = 0.01$ s and $\Delta t = 0.001$ s). Both the implicit scheme ($\theta = 1$) and the Crank–Nicolson scheme ($\theta = 0.5$) are considered.

The distribution of the temperature is represented in Fig. 1a at several time steps and the evolution of the temperature at the center of the cell ($x = 0.5$ m) is shown in 1b. The numerical solution converges towards the analytical solution by decreasing the time step increment $\Delta t$; the results obtained using the Crank–Nicolson scheme appears to be more precise than the standard implicit scheme, for a given time step increment. This analytical solution allows a validation of the numerical scheme.

## 5. Application of the method to a heterogeneous case

### 5.1. Description of the simulations

The method is applied to the prediction of transient thermal conduction in a composite laminate material made of two phases of properties $c_1$, $\rho_1$, $c_{p1}$ in phase 1 and $c_2$, $\rho_2$, $c_{p2}$ in phase 2. The microstructure is 2-d square domain $[0, L] \times [0, L]$ composed of 13 layers of phase 1, 12 layers of phase 2 of thickness $5/128L$ and two exterior layers of phase 2 of
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Description of the problem considered. (a) Sketch of the microstructure, (b) Function G and (c) Function Amp.</figcaption></figure>

thickness 3/256L, as shown in Fig. 2. This square domain is discretized with 129 × 129 grid points, corresponding to 128 × 128 pixels.

The cell is subjected to the following boundary and initial conditions:

$$ 
\left\{ 
\begin{array}{l}
T(\mathbf{x}, 0) = T_0 \\
T_{\text{top}}(\mathbf{x}, t) = T_{\text{bottom}}(\mathbf{x}, t) = T_{\text{right}}(y, t) = T_0 \\
T_{\text{left}}(y, t) = T_0 + (T_{\text{max}} - T_0)G(y) \times \text{Amp}(t),
\end{array}
\right. 
\quad (64) $$

where G(y) and Amp(t) are non-dimensional (respectively spatial and temporal) functions represented in Fig. 2 and defined by

$$ G(y) = \begin{cases} 0 & \text{if } y \le L/4 \\ 0 & \text{if } y \ge 3L/4 \\ \frac{1}{2} \left( 1 + \cos \left[ \frac{4\pi}{L} \left(y - \frac{L}{2}\right) \right] \right) & \text{if } L/4 < y < 3L/4 \end{cases} \quad (65) $$

and

$$ \text{Amp}(t) = \begin{cases} \sin \left( \pi \frac{t}{t_{\text{max}}} \right) & \text{if } t \le t_{\text{max}}/2 \\ 1 & \text{if } t \ge t_{\text{max}}/2. \end{cases} \quad (66) $$

These boundary conditions correspond to a full Dirichlet case so the fluctuation function $\tilde{T}$ is described by a sine-sine series.

The parameters taken for the simulation are as follows:

* Spatial and time domain: $L_x = L_y = L = 1$ m, $\Delta x = \Delta y = 0.0079$ m, $t_{\text{max}} = 1$ s;

* Temperature conditions: $T_0 = 20$ °C, $T_{\text{max}} = 200$ °C;

* Phase 1: $c_1 = 1$ W m⁻¹ K⁻¹, $c_{p1} = 1$ J kg⁻¹ K⁻¹, $\rho_1 = 1$ kg m⁻³;

* Phase 2: $c_2 = 0.1$ W m⁻¹ K⁻¹, $c_{p2} = 0.1$ J kg⁻¹ K⁻¹, $\rho_2 = 1$ kg m⁻³.

The values taken for the phases' properties do not correspond to any real material but are taken to highlight the influence of material contrast (i.e. the microstructure heterogeneity) on the evolution of the local fields; in the present case, the contrast is $c_2/c_1 = (\rho_2 c_{p2})/(\rho_1 c_{p1}) = 0.1$.

Simulations are performed using an Euler implicit scheme ($\theta = 1$) and a Crank–Nicolson scheme ($\theta = 0.5$) as these cases represent the most interesting situations for the method. Several time step sizes are considered to investigate the predictions: $\Delta t = 0.01$ s and $\Delta t = 0.001$ s. It should be noted that the time steps considered for this application are a few order of magnitude higher than the minimal value $\Delta t_{\text{explicit}} = \frac{\Delta x^2}{2} \times \min(\frac{\rho c_p}{c}) \approx 3.12 \times 10^{-5}$ s of the explicit scheme (in order to verify the CFL condition).

In all calculations, the tolerance for the iterative procedure is taken as $\text{tol} = 10^{-6}$.

## 5.2. Numerical results

Implicit scheme ($\theta = 1$). We begin with the Euler implicit scheme. The distribution of the temperature field T and the component $q_x$ of the
<figure><img src="image_6.png" /><figcaption>Fig. 3. Distribution of the temperature in the case $\theta = 1$ (implicit) and $\Delta t = 0.01$ s.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Distribution of the component $q_x$ of the heat flux in the case $\theta = 1$ (implicit) and $\Delta t = 0.01$ s.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Distribution of the temperature in the case $\theta = 0.5$ (Crank–Nicolson) and $\Delta t = 0.01$ s.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Distribution of the component $q_x$ of the heat flux in the case $\theta = 0.5$ (Crank–Nicolson) and $\Delta t = 0.01$ s.</figcaption></figure>

heat flux are represented, in the case $\Delta t = 0.01$ s, at several snapshots ($t = [0.2, 0.3, 0.4, 1]$ s) respectively in Figs. 3 and 4. The distribution of both quantities highlight the effect of the (laminate) microstructure. Interestingly, no numerical artefacts (such as oscillations) are observed on the temperature field and very moderate oscillations can be observed on the heat flux, at some interfaces, which is standard using spectral differentiation rules (see e.g. Morin et al. (2021)); it should be noted that these oscillations do not accumulate during the simulation. The effect of the time step is not visible on the distributions and will be investigated hereafter on more quantitative plots. The number of iterations to convergence is rather constant during each simulation and

it does not depend much on the time step discretization; in both cases $\Delta t = 0.01$ s and $\Delta t = 0.01$ s, the number of iterations to convergence is of about 45.

Crank–Nicolson scheme ($\theta = 0.5$). We continue with the Crank–Nicolson scheme. The distribution of the temperature field $T$ and the component $q_x$ of the heat flux are again represented, in the case $\Delta t = 0.01$ s, at several snapshots ($t = [0.2, 0.3, 0.4, 1]$ s) respectively in Figs. 5 and 6. Only the results for the time step value $\Delta t = 0.01$ s are represented because the effect of $\Delta t$ is negligible in that case. As in the Euler implicit scheme, the number of iterations to convergence is rather constant during each simulation and it does not depend much
<figure><img src="image_10.png" /><figcaption>Fig. 7. Influence of the time step on (a) The evolution of the temperature at the center of the cell and (b) The evolution of the component $q_x$ of the heat flux at the center of the cell.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Comparison with FEM (reference) results. (a) Evolution of the temperature at the center of the cell and (b) Evolution of the component $q_x$ of the heat flux at the center of the cell.</figcaption></figure>

on the time step discretization; in both cases $\Delta t = 0.01$ s and $\Delta t = 0.01$ s, the number of iterations to convergence is of about 40 for the Crank–Nicolson scheme.

Effect of the time discretization $\Delta t$. The effect of the time discretization is now investigated. The evolution of the temperature and the component $q_x$ of the heat flux are represented at the center of the cell (i.e. for $x = y = L/2$) in Fig. 7. The effect of the time step is moderate using the implicit scheme and almost negligible using the Crank–Nicolson scheme, i.e. coarse time steps using the Crank–Nicolson scheme provide an accurate solution. This is expected as the Crank–Nicolson scheme is a second-order method in time while the implicit scheme is a first-order method in time.

### 5.3. Comparison with FEM calculations

The (heterogeneous) test case examined in this section lacks an analytical solution. Therefore, to assess the validity of the results obtained with the newly developed FFT procedure, a numerical reference solution was calculated; for this purpose, we relied on the well-established finite element method (FEM), using a finite element code validated in challenging scenarios such as crack propagation (Essongue et al., 2023) and thermal modeling of additive manufacturing (Essongue et al., 2024). To ensure the FEM solution's accuracy, we utilized a refined mesh with 0.5 mm elements (fifteen times finer than the FFT grid)

and a time step size of 0.1 ms (ten times smaller than that used for FFT calculations). Robust, parameter-free choices were made, including the use of bilinear quadrilateral elements, the imposition of boundary conditions through Lagrange multipliers, a backward Euler time discretization scheme, and a direct solver to handle the resulting equations. As shown in Fig. 8, the FFT results align perfectly with the reference solution, thereby validating the newly developed FFT-based solver.

## 6. Conclusion

The aim of this work was to develop efficient FFT-based solvers for transient diffusion in heterogeneous materials subjected to non-periodic (Dirichlet/Neumann) boundary conditions. We focus on a problem of thermal conductivity, which constitutes the prototype equation of any diffusion process. The numerical method proposed in this work relies on the use of (i) theta-method (which includes Euler implicit, Euler explicit and Crank–Nicolson schemes) for time discretization and (ii) a Galerkin-based discretization using an approximation space based on mixed sine–cosine series for the spatial discretization. The method is based on a fixed-point iterative solution of an auxiliary problem involving homogeneous material properties, following the pioneering work of Moulinec and Suquet (1998). The solution field is decomposed into a known term verifying the non-periodic boundary conditions and a fluctuation (unknown) term described by appropriate sine–cosine
series. The solution of this auxiliary problem involves discrete sine-cosine transforms, of type I and III, depending on the type of boundary conditions. The system arising in this formulation is diagonal and therefore no linear system needs to be solved; the solver relies on the computational complexity of fast Fourier transforms. The method has finally been applied to the prediction of transient thermal fields in a composite material subjected to non periodic boundary conditions.

The method developed in this work investigated, in a Galerkin-based framework, the problem considered by Gélébart (2025) using a finite difference scheme in space. Both methods rely on discrete sine-cosine transforms and but mainly differ by the Green operator arising in the solution of the auxiliary problem. In both works, the classical oscillations of the FFT-based framework are observed using spectral derivation or hexahedral-type schemes (Willot et al., 2014; Morin et al., 2021; Gélébart, 2024). One might remark that Gélébart (2025) also used a tetrahedral stencil, recently proposed by Finel (2024), which appears to suppress most of the spurious oscillations. Additionally, the method developed in this work is based on the calculation of the temperature (or equivalently the displacement field in elasticity problems) in contrast with standard FFT-based method which operates in the gradient fields (see e.g. Vondřejc et al. (2014)); therefore accelerated schemes (Eyre and Milton, 1999; Michel et al., 2000; Monchiet and Bonnet, 2012; Moulinec and Silva, 2014; Kabel et al., 2014) cannot be used as is and need to be adapted to the displacement-based formulation.

Future developments of this work will concern the coupling between a diffusion and mechanics, to handle thermomechanical problems as well as multiphysics problems such as, e.g., hydrogen-assisted deformation of polycrystalline materials (Singh et al., 2022), stress-dependent polymer ageing (Konica and Sain, 2021) or electro-chemical mechanical coupling in lithium-ion batteries (Zhao et al., 2019). In addition, the methods developed in this work could benefit the simulation of processes such additive manufacturing for which efficient numerical strategies are required (see e.g. Essongue et al. (2024)).

### CRediT authorship contribution statement

**Abdoul Magid Amadou Sanoko:** Writing – review & editing, Writing – original draft, Software, Investigation. **Simon Essongue:** Writing – review & editing, Writing – original draft, Software, Investigation. **Lionel Gélébart:** Writing – review & editing, Writing – original draft, Software, Investigation. **Lucas Lapostolle:** Writing – review & editing, Writing – original draft, Software, Investigation. **Léo Morin:** Writing – review & editing, Writing – original draft, Software, Investigation. **Joseph Paux:** Writing – review & editing, Writing – original draft, Software, Investigation.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Appendix. Elementary integrals

The integrals $I_{lmn}$, $J_{lmn}$ and $K_{lmn}$ arising in the weak formulation (42) of the transient conductivity problem are calculated as follows.

* The integral $I_{lmn}$ reads

$$I_{lmn} = \sum_{i=0}^{N_x+1} \sum_{j=0}^{N_y+1} \sum_{k=0}^{N_z+1} \left[ \left( \theta c_0 \left( (\xi_j^x)^2 + (\xi_j^y)^2 + (\xi_k^z)^2 \right) + \frac{\rho_0 c_p 0}{\Delta t} \right) \alpha_i^x \alpha_j^y \alpha_k^z \tilde{\tau}_{lmn} \times \int_{x=0}^{x=L_x} g^x(k_i^x) g^x(k_j^x) dx \times \int_{y=0}^{y=L_y} g^y(k_j^y) g^y(k_m^y) dy \times \int_{z=0}^{z=L_z} g^z(k_k^z) g^z(k_n^z) dz \right]. \quad (A.1)$$

After calculation of the elementary integrals in Eq. (A.1), it can be shown that

$$I_{lmn} = \frac{L_x L_y L_z}{8} \left( \theta c_0 \left( (\xi_j^x)^2 + (\xi_m^y)^2 + (\xi_n^z)^2 \right) + \frac{\rho_0 c_p 0}{\Delta t} \right) \tilde{\tau}_{lmn}, \quad (A.2)$$

since

$$\int_{x=0}^{x=L_x} g^x(k_i^x) g^x(k_j^x) dx = \frac{\delta_{il} L_x}{2 \alpha_i^x}, \quad \int_{y=0}^{y=L_y} g^y(k_j^y) g^y(k_m^y) dy = \frac{\delta_{jm} L_y}{2 \alpha_m^y}, \quad (A.3)$$

$$\int_{z=0}^{z=L_z} g^z(k_k^z) g^z(k_n^z) dz = \frac{\delta_{kn} L_z}{2 \alpha_n^z}, \quad (A.4)$$

with $\delta$ the Kronecker symbol.

* The integral $J_{lmn}$ reads

$$J_{lmn} = J_{lmn}^x + J_{lmn}^y + J_{lmn}^z, \quad (A.5)$$

with

$$J_{lmn}^x = \int_{\Omega} \left( \frac{\partial B_x}{\partial x} \right) g^x(k_i^x) g^y(k_m^y) g^z(k_n^z) dz \Omega, \quad (A.6)$$

$$J_{lmn}^y = \int_{\Omega} \left( \frac{\partial B_y}{\partial y} \right) g^x(k_i^x) g^y(k_m^y) g^z(k_n^z) dz \Omega, \quad (A.7)$$

$$J_{lmn}^z = \int_{\Omega} \left( \frac{\partial B_z}{\partial z} \right) g^x(k_i^x) g^y(k_m^y) g^z(k_n^z) dz \Omega, \quad (A.8)$$

where $B_x$, $B_y$ and $B_z$ are given by Eq. (38). In the definition of the auxiliary problem, the known fields are $B_x$, $B_y$ and $B_z$, and not the partial derivatives involved in Eqs. (A.6), (A.7) and (A.8). Consequently, as in the stationary case (Paux et al., 2025a), the integral $J_{lmn}$ can be calculated using integration by parts in order to avoid calculation of the aforementioned partial derivatives. In the case of $J_{lmn}^x$ (the other integrals can be calculated similarly), one has

$$J_{lmn}^x = \int_{y=0}^{y=L_y} \int_{z=0}^{z=L_z} \left( g^x(k_i^x L_x) B_x(L_x, y, z) - g^x(0) B_x(0, y, z) \right) \times g^y(k_m^y) g^z(k_n^z) dz dy - \int_{\Omega} \frac{\delta^x}{\delta z} \int_{z=0}^{z=L_z} B_x(x, y, z) g^x(k_i^x) g^y(k_m^y) g^z(k_n^z) dz \Omega. \quad (A.9)$$

Let us introduce the function $\tilde{B}_x$, whose expression depends on the type of boundary conditions on the x-faces:

$$\left\{ \begin{aligned} \tilde{B}_x(x, y, z) &= B_x(x, y, z) & (\text{DD}) \\ \tilde{B}_x(x, y, z) &= B_x(x, y, z) - B_x(0, y, z) - \frac{B_x(L_x, y, z) - B_x(0, y, z)}{L_x} & (\text{NN}) \\ \tilde{B}_x(x, y, z) &= B_x(x, y, z) - B_x(L_x, y, z) & (\text{DN}) \\ \tilde{B}_x(x, y, z) &= B_x(x, y, z) - B_x(0, y, z) & (\text{ND}) \end{aligned} \right. \quad (A.10)$$

The calculation of $J_{lmn}^x$ then reduces to

$$J_{lmn}^x = - \frac{\delta^x}{\delta z} \int_{\Omega} \tilde{B}_x(x, y, z) \overline{g^x(k_i^x) g^y(k_m^y) g^z(k_n^z)} dz \Omega + s_{lmn}^x, \quad (A.11)$$

where $s_{lmn}^x$ is related to the constant term arising in the derivative in the NN case and will be detailed hereafter. For arbitrary functions $\tilde{B}_x$, $\tilde{B}_y$ and $\tilde{B}_z$, the elementary integral defined in Eq. (A.11) cannot be calculated analytically; therefore those integrals are calculated approximately, using discrete sine-cosine transforms:

$$J_{lmn}^x = \frac{-L_x L_y L_z}{(N_x + 1)(N_y + 1)(N_z + 1)} \left[ \frac{\xi_l^x}{D_{xyz}(\tilde{B}_x)} \right]_{lmn} + S_{lmn}^x, \quad (A.12)$$

where $\tilde{B}_x$ denotes the array of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of the function $\tilde{B}_x$ at the grid points. The term $S_{lmn}^x$ is always null, expect in the case $l = 0$ of the NN boundary conditions, where it is given by

$$S_{0mn}^x = - \frac{(N_x + 1)}{L_x} (D_{yz}(\tilde{B}_x^0 - \tilde{B}_x^{L_x}))_{mn}. \quad (A.13)$$
The term $J_{lmn}$ finally reads

$$J_{lmn} = \frac{-L_x L_y L_z}{(N_x + 1)(N_y + 1)(N_z + 1)} \left[ \begin{aligned} & \xi_l^x (D_{\bar{x}\bar{y}}(\dot{\mathbf{B}}_x))_{lmn} \\ & + \xi_m^y (D_{\bar{x}\bar{y}}(\dot{\mathbf{B}}_y))_{lmn} \\ & + \xi_n^z (D_{\bar{x}\bar{z}}(\dot{\mathbf{B}}_z))_{lmn} + S_{lmn} \end{aligned} \right] \quad (A.14)$$

where $\dot{\mathbf{B}}_x, \dot{\mathbf{B}}_y$ and $\dot{\mathbf{B}}_z$ respectively denote the arrays of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of functions $\dot{\mathbf{B}}_x, \dot{\mathbf{B}}_y$ and $\dot{\mathbf{B}}_z$ at the grid points, and $S_{lmn}$ is given by

$$S_{lmn} = 0 \qquad \text{if } l \neq 0 \text{ and } m \neq 0 \text{ and } n \neq 0 \quad (A.15)$$

$$\begin{aligned} S_{0mn} &= -\frac{(N_x + 1)}{L_x} (D_{yz}(\mathbf{B}_x^0 - \mathbf{B}_x^{L_x}))_{mn} & \text{if NN on } x \text{-faces} \\ S_{l0n} &= -\frac{(N_y + 1)}{L_y} (D_{xz}(\mathbf{B}_y^0 - \mathbf{B}_y^{L_y}))_{ln} & \text{if NN on } y \text{-faces} \\ S_{lm0} &= -\frac{(N_z + 1)}{L_z} (D_{xy}(\mathbf{B}_z^0 - \mathbf{B}_z^{L_z}))_{lm} & \text{if NN on } z \text{-faces,} \end{aligned} \quad (A.16)$$

where $\mathbf{B}_x^0, \mathbf{B}_x^{L_x}, \mathbf{B}_y^0, \mathbf{B}_y^{L_y}, \mathbf{B}_z^0$ and $\mathbf{B}_z^{L_z}$ denote the arrays of sizes $(N_x + 2) \times (N_y + 2)$ for $\mathbf{B}_x^0$ and $\mathbf{B}_x^{L_x}$, $(N_x + 2) \times (N_z + 2)$ for $\mathbf{B}_y^0$ and $\mathbf{B}_y^{L_y}$, and $(N_x + 2) \times (N_y + 2)$ for $\mathbf{B}_z^0$ and $\mathbf{B}_z^{L_z}$, containing the values of functions $B_x, B_y$ and $B_z$ on the faces.

• Finally, the integral $K_{lmn}$ reads

$$K_{lmn} = \int_{\Omega} A(x) g^x(k_l^x x) g^y(k_m^y y) g^z(k_n^z z) d\Omega, \quad (A.17)$$

where A is given by Eq. (38). The integral (A.17) is then calculated approximately using discrete sine-cosine transforms as

$$K_{lmn} = \frac{L_x L_y L_z}{(N_x + 1)(N_y + 1)(N_z + 1)} (D_{xyz}(A))_{lmn}, \quad (A.18)$$

where A is the array of size $(N_x + 2) \times (N_y + 2) \times (N_z + 2)$ containing the values of the function A at the grid points.

## Data availability

Data will be made available on request.

## References

Bertin, N., Upadhyay, M.V., Pradalier, C., Capolungo, L., 2015. A FFT-based formulation for efficient mechanical fields computation in isotropic and anisotropic periodic discrete dislocation dynamics. Modelling Simul. Mater. Sci. Eng. 23, 065009.

Bilger, N., Auslender, F., Bornert, M., Michel, J.C., Moulinec, H., Suquet, P., Zaoui, A., 2005. Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis. Int. J. Solids Struct. 42, 517–538.

Brenner, R., 2009. Numerical computation of the response of piezoelectric composites using Fourier transform. Phys. Rev. B 79, 184106.

Brenner, R., Beaudoin, A.J., Suquet, P., Acharya, A., 2014. Numerical implementation of static field dislocation mechanics theory for periodic media. Phil. Mag. 94, 1764–1787.

Chakraborty, A., Lebensohn, R.A., Capolungo, L., 2023. Coupled chemo-mechanical modeling of point-defect diffusion in a crystal plasticity fast fourier transform framework. J. Mech. Phys. Solids 173, 105190.

Essongue, S., Couégnat, G., Martin, E., 2023. On the use of drilling degrees of freedom to stabilise the augmented finite element method. Appl. Mech. 4, 1140–1171.

Essongue, S., Nain, V., Carin, M., 2024. Efficient thermal modeling of laser directed energy deposition using the forward euler scheme: Methodology, merits and limitations. Finite Elem. Anal. Des. 242, 104270.

Eyre, D.J., Milton, G.W., 1999. A fast numerical scheme for computing the response of composites using grid refinement. Eur. Phys. J. - Appl. Phys. 6, 41–47.

Finel, A., 2024. A fast and robust discrete FFT-based solver for computational homogenization. arXiv preprint arXiv:2405.11168.

Frigo, M., Johnson, S.G., 1998. FFTW: An adaptive software architecture for the FFT. In: Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP'98 (Cat. No. 98CH36181). IEEE, pp. 1381–1384.

Gélébart, L., 2020. A modified FFT-based solver for the mechanical simulation of heterogeneous materials with Dirichlet boundary conditions. C. R. Mec. 348, 693–704.

Gélébart, L., 2022. Amitex. https://amitexftft.github.io/AMITEX/index.html.

Gélébart, L., 2024. FFT-based simulations of heterogeneous conducting materials with combined non-uniform Neumann, periodic and Dirichlet boundary conditions. Eur. J. Mech. A Solids 105248.

Gélébart, L., 2025. An accurate and robust fft-based solver for transient diffusion in heterogeneous materials. C. R. Mec. 353, 113–125.

Grimm-Strele, H., Kabel, M., 2021. FFT-based homogenization with mixed uniform boundary conditions. Internat. J. Numer. Methods Engrg. 122, 7241–7265.

Kabel, M., Böhlke, T., Schneider, M., 2014. Efficient fixed point and Newton–krylov solvers for FFT-based homogenization of elasticity at large deformations. Comput. Mech. 54, 1497–1514.

Kaptchouang, N.B.N., Gélébart, L., 2022. Multiscale coupling of FFT-based simulations with the Idc approach. Comput. Methods Appl. Mech. Engrg. 394, 114921.

Konica, S., Sain, T., 2021. A reaction-driven evolving network theory coupled with phase-field fracture to model polymer oxidative aging. J. Mech. Phys. Solids 150, 104347.

Lebensohn, R., 2001. N-site modeling of a 3D viscoplastic polycrystal using fast Fourier transform. Acta Mater. 49, 2723–2737.

Lucarini, S., Upadhyay, M.V., Segurado, J., 2021. FFT based approaches in micromechanics: fundamentals, methods and applications. Modelling Simul. Mater. Sci. Eng. 30, 023002.

Michel, J.C., Moulinec, H., Suquet, P., 2000. A computational method based on augmented Lagrangians and fast Fourier transforms for composites with high contrast. Comput. Model. Eng. Sci. 1, 79–88.

Monchiet, V., Bonnet, G., 2012. A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast. Internat. J. Numer. Methods Engrg. 89, 1419–1436.

Monchiet, V., Bonnet, G., 2024. FFT based iterative schemes for composite conductors with uniform boundary conditions. Eur. J. Mech. A Solids 103, 105146.

Morin, L., Brenner, R., Derrien, K., Dorhmi, K., 2021. Periodic smoothing splines for FFT-based solvers. Comput. Methods Appl. Mech. Engrg. 373, 113549.

Morin, L., Paux, J., 2024. A fast numerical method for the conductivity of heterogeneous media with Dirichlet boundary conditions based on discrete sine–cosine transforms. Comput. Methods Appl. Mech. Engrg. 421, 116772.

Moulinec, H., Silva, F., 2014. Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials. Internat. J. Numer. Methods Engrg. 97, 960–985.

Moulinec, H., Suquet, P., 1998. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput. Methods Appl. Mech. Engrg. 157, 69–94.

Paux, J., Brenner, R., Kondo, D., 2018. Plastic yield criterion and hardening of porous single crystals. Int. J. Solids Struct. 132–133, 80–95.

Paux, J., Morin, L., Gélébart, L., 2025a. A discrete sine-cosine transforms galerkin method for the conductivity of heterogeneous materials with mixed dirichlet/neumann boundary conditions. Internat. J. Numer. Methods Engrg. 126, e7615.

Paux, J., Morin, L., Gélébart, L., Sanoko, A.M.A., 2025b. A discrete sine–cosine based method for the elasticity of heterogeneous materials with arbitrary boundary conditions. Comput. Methods Appl. Mech. Engrg. 433, 117488.

Pundir, M., Kammer, D.S., Angst, U., 2023. An FFT-based framework for predicting corrosion-driven damage in fractal porous media. J. Mech. Phys. Solids 179, 105388.

Risthaus, L., Schneider, M., 2024a. FFT-based computational micromechanics with dirichlet boundary conditions on the rotated staggered grid. Internat. J. Numer. Methods Engrg. 125, e7569.

Risthaus, L., Schneider, M., 2024b. Imposing different boundary conditions for thermal computational homogenization problems with FFT- and tensor-train-based green's operator methods. Internat. J. Numer. Methods Engrg. 125, e7423.

Risthaus, L., Schneider, M., 2024c. Imposing Dirichlet boundary conditions directly for FFT-based computational micromechanics. Comput. Mech. 74, 1089–1113.

Sancho, R., Rey-de Pedraza, V., Lafourcade, P., Lebensohn, R.A., Segurado, J., 2023. An implicit FFT-based method for wave propagation in elastic heterogeneous media. Comput. Methods Appl. Mech. Engrg. 404, 115772.

Schneider, M., 2021. A review of nonlinear FFT-based computational homogenization methods. Acta Mech. 232, 2051–2100.

Singh, V., Kumar, R., Charles, Y., Mahajan, D.K., 2022. Coupled diffusion-mechanics framework for simulating hydrogen assisted deformation and failure behavior of metals. Int. J. Plast. 157, 103392.

Vidyasagar, A., Tan, W.L., Kochmann, D.M., 2017. Predicting the effective response of bulk polycrystalline ferroelectric ceramics via improved spectral phase field methods. J. Mech. Phys. Solids 106, 133–151.

Vondřejc, J., Zeman, J., Marek, I., 2014. An FFT-based galerkin method for homogenization of periodic media. Comput. Math. Appl. 68, 156–173.
Willot, F., Abdallah, B., Pellegrini, Y.P., 2014. Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields. Internat. J. Numer. Methods Engrg. 98, 518–533.

Zarzoso, G., Roque, E., Montero-Chacón, F., Segurado, J., 2025. An FFT based chemo-mechanical framework with fracture: Application to mesoscopic electrode degradation. Mech. Mater. 201, 105211.

Zecevic, M., Lebensohn, R.A., 2025. Extended FFT-based micromechanical formulation to consider general non-periodic boundary conditions. Int. J. Solids Struct. 311, 113225.

Zhao, Y., Stein, P., Bai, Y., Al-Siraj, M., Yang, Y., Xu, B.X., 2019. A review on modeling of electro-chemo-mechanics in lithium-ion batteries. J. Power Sources 413, 259–283.