# A simple extension of FFT-based methods to strain gradient loadings - Application to the homogenization of beams and plates with linear and non-linear behaviors

Lionel Gélébart

► To cite this version:

Lionel Gélébart. A simple extension of FFT-based methods to strain gradient loadings - Application to the homogenization of beams and plates with linear and non-linear behaviors. Journal of Theoretical, Computational and Applied Mechanics, 2022, pp.1-27. 10.46298/jtcam.6790. hal-02942202v3

HAL Id: hal-02942202

https://hal.science/hal-02942202v3

Submitted on 1 Jul 2022

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.6790

OAI hal-02942202v3

**History**

Received Sep 18, 2020

Accepted May 19, 2021

Published Jun 16, 2022

**Associate Editor**

Julien RÉTHORÉ

**Reviewers**

Anonymous
Sébastien BRISARD

**Open Review**

OAI hal-03667777

**Licence**

CC BY 4.0
©The Authors

# A simple extension of FFT-based methods to strain gradient loadings – Application to the homogenization of beams and plates with linear and non-linear behaviors

Lionel GÉLÉBART

Université Paris-Saclay, CEA, ISAS/DMN/SRMA, 91191, Gif/Yvette, France

Because of their simplicity, efficiency and ability for parallelism, FFT-based methods are very attractive in the context of numerical periodic homogenization, especially when compared to standard FE codes used in the same context. They allow applying to a unit-cell a uniform average strain with a periodic strain fluctuation that is an unknown quantity. Solving the problem allows to evaluate the complete stress-strain fields. The present work extends straightforwardly the method from uniform loadings (i.e. uniform applied strain) to strain gradient loadings (i.e. strain fields with a uniform strain gradient) and proposes an application to the homogenization of beams and plates, while keeping the simplicity of the original method. It is observed that, due to the use of periodic boundary conditions, among the 18 strain gradient components, only 9 can be really prescribed. However, they are sufficient to consider the homogenization of beams and plates, taking into account stress-free boundary conditions through enlarged unit-cells with null stiffness. A first application validates the approach and compares it to another FFT-based method dedicated to the homogenization of plates. The second application concerns the homogenization of beams, for the first time considered (to author's knowledge) with an FFT-based solver. The method applies to different beam cross-sections and the proposition of using composite voxels drastically improves the numerical solution when the beam cross-section is not conform with the spatial discretization, especially for torsion loading. Finally, the method is applied to non-linear behaviors, considering a beam unit-cell with an elastic inclusion within an elasto-plastic matrix, perfectly plastic. As a result, the massively parallel AMITEX_FFTP code has been slightly modified and now offers a new functionality, allowing the users to prescribe torsions and flexions to beam or plate heterogeneous unit-cells, with both linear and non-linear behaviors.

**Keywords:** numerical homogenization; FFT; strain gradients; beams; plates; flexion; torsion; composite voxels; elasto-plasticity

## 1 Introduction

The present paper lies in the overall framework of numerical homogenization of heterogeneous materials for which the mechanical behavior at the upper scale derives from spatial averages of stress-strain fields simulated on heterogeneous unit-cells. The first developments and most common applications focus on first order homogenization, the homogenized behavior relating the displacement gradient tensor to a stress tensor, or more specifically, under the small perturbations assumption, the linearized strain to the Cauchy stress. The type of boundary conditions applied to the unit-cell is quite important in that context. First, they must satisfy the Hill-Mandel condition so that the average microscopic strain energy is equal to the macroscopic one. Then, various choices satisfying this condition are available: Kinematic Uniform Boundary Conditions (BC), Stress Uniform BC, periodic BC or even normal-mixed BC (Gélébart et al. 2009). If the choice of periodic BC is obvious for periodic microstructures, it also provides the best estimate for random microstructures (Kanit et al. 2003; Chateau et al. 2015) compared to KUBC or SUBC, that
provide upper and lower bounds, at least for linear material (Huet 1990). If applying PBC is not straightforward in Finite Element codes and can deteriorate the computation performance, it is natural for the FFT-based methods proposed in the 1990's (Moulinec and Suquet 1998).

These methods have been extensively used and improved during the last decade, reducing the initial drawbacks such as convergence issues for highly contrasted materials, spurious oscillations or sensitivity of the convergence to the reference material behavior (an algorithm parameter). Among significant improvements: modified discrete Green operators reducing spurious oscillations and improving the convergence (Brisard and Dormieux 2010; Willot 2015; Schneider et al. 2017), algorithms replacing the initial fixed-point algorithm to improve convergence (Zeman et al. 2010; Gélébart and Mondon-Cancel 2013; Kabel et al. 2014; Chen et al. 2019b) or composite voxels accounting for multi-phase materials, for voxels crossed by an interface, to improve the numerical solution and reduce spurious oscillations (Brisard and Dormieux 2010; Gélébart and Ouaki 2015; Kabel et al. 2015). The list is not exhaustive and numerous theoretical works reinforce the mathematical basis of the method.

Due to their high numerical efficiency compared to finite element codes used in the same context (periodic boundary conditions), various works intend to extend the method to various applications. Among these applications, we can cite the extension to non-local behaviors such as strain gradient plasticity (Lebensohn and Needleman 2016) or damage phase field (Chen et al. 2019b), the coupling with metallurgical phase fields such as martensitic phase transformation (Kochmann et al. 2015) or other physics such as magneto-electricity (Brenner and Bravo-Castillero 2010), the coupling with Dislocation Dynamics (DD) codes (Bertin and Capolungo 2018). Once again, this short list is far from exhaustive but gives an idea of the increasing development of the method towards new research fields.

Another kind of extension concerns higher order numerical homogenization (Tran et al. 2012; Dietrich et al. 2019). In that case, the homogenized behavior is not only dependent of the first gradient of the displacement but also of higher order terms. Following this idea of going beyond the classical application of an average strain loading (i.e. first order displacement gradient) the present paper lies in the same general context.

The first part of the paper proposes a straightforward extension of the classical FFT-based method in order to account for strain gradient loadings. It is observed that, due to the use of periodic boundary conditions, among the 18 strain gradient components, only 9 can be applied. The second part of the paper demonstrates its ability to solve the question of the homogenization of beams and plates, taking into account the traction free boundary conditions through enlarged unit-cells with null stiffness.

Actually, the questions of beams and plates homogenization can be regarded as second order term homogenizations (Geers et al. 2007) for which strain gradients are associated to flexion and torsion loadings. For plates, if the classical numerical implementation relies on the finite element method (Geers et al. 2007; Helfen and Diebels 2014), Nguyen et al. (2008) also proposes an FFT-based method devoted to this application. It relies on a specific algorithm combined with a dedicated Green operator to account for the traction free boundary conditions applied at the plates' free surfaces. On the contrary, if the finite element method has also been used is also used for beams homogenization (Cartraud and Messager 2006), to the best of the author's knowledge, FFT-based methods have never been proposed in that context. One purpose of the paper is then to demonstrate that the same slightly modified FFT-based code, used with enlarged unit-cells with null stiffness, can be used for both classical (i.e. first gradient term), as well as beams and plates, numerical periodic homogenization. This is made possible and quite efficient since recent FFT-based methods (see beginning of the introduction) allow for the simulation of unit-cells enlarged with void voxels (i.e. infinite contrast), to account for traction free boundary conditions.

To summarize, the first part of the paper (Section 2) introduces the simple extension of the FFT-based method to a reduced set of strain gradient loadings and the second part (Section 3) applies the method to the homogenization of beams and plates. In this part, the first sub-section presents the application to plates homogenization and focuses on the comparison with Nguyen's method using exactly the same cases (Nguyen et al. 2008). Then, the second sub-section presents the new FFT-based application to beams homogenization with an additional emphasis on the use of composite voxels when the beam cross-section is not conform with the regular grid
discretization. Finally, a last section is proposed to demonstrate that the method can be applied to non-linear behaviors: a beam unit-cell with an elastic inclusion within an elasto-plastic matrix, perfectly-plastic, is submitted to bending and torsion.

The method proposed in (Nguyen et al. 2008), based on a modified Green operator, is limited to plates with plane free surfaces and its application was proposed in the context of linear elasticity. However, the present work makes use of unit-cells enlarged with void voxels and thus accounts for both plates and beams with any arbitrary cross-section. For plates, the free surfaces do not have to be plane and corrugated plates could be simulated as well. The FFT-based simulation of beams, which is itself a novelty, is enriched with the use of Voigt composite voxels to improve the quality of the simulations (here, for beams with disk cross-sections). In addition, the algorithm is extended to non-linear behaviors and validated by an application to a heterogeneous beam with an elasto-plastic matrix. Finally, the discussion on the strain gradient components that can be prescribed when using periodic boundary conditions (i.e. 9 components among 18) could also be of interest for higher order numerical homogenization.

## 2 FFT-based method extension to strain gradient loadings

### 2.1 Periodicity condition

The description of the periodicity condition prescribed on a unit-cell $\Omega$ can be done equivalently from two different points of view. In the first case, fields are defined on $\Omega$ and boundary conditions have to be applied on its boundary $\partial\Omega$. In the second case, fields are defined on an infinite space and $\Omega$-periodicities are assumed for the fields. The first description is more appropriate for Finite-Element solvers, the second is well-suited for FFT-based solvers. The two equivalent descriptions of the equilibrium condition, submitted to a periodicity condition, are given as

$$ \left\{ \begin{array}{l} \mathrm{div}(\sigma(\mathbf{x})) = 0 \\ \sigma \cdot \mathbf{n} \text{ antiperiodic on } \partial\Omega \end{array} \right. \quad \Leftrightarrow \quad \left\{ \begin{array}{l} \mathrm{div}(\sigma(\mathbf{x})) = 0 \\ \sigma \cdot \Omega \text{-periodic} \end{array} \right. \qquad (1) $$

where “$\sigma \cdot n$ antiperiodic on $\partial\Omega$” means that the traction vector $\sigma \cdot n$ on opposite points of opposite faces of $\Omega$ are of opposite sign. The periodicity condition of a displacement field $\tilde{u}$ in

$$ \tilde{u} \text{ periodic on } \partial\Omega \quad \Leftrightarrow \quad \tilde{u} \cdot \Omega\text{-periodic} \quad \Leftrightarrow \quad \left\{ \begin{array}{l} \tilde{\varepsilon} \text{ compatible} \\ \tilde{\varepsilon} \cdot \Omega\text{-periodic} \\ \langle \tilde{\varepsilon} \rangle = 0 \end{array} \right. \qquad (2) $$

is also equivalently formulated on $\Omega$ in terms of strain field $\tilde{\varepsilon}$, adding compatibility and null average conditions in that case. Below, “$\tilde{u}$ periodic on $\partial\Omega$” means that the displacements $\tilde{u}$ evaluated on opposite points of opposite faces are equal.

### 2.2 FFT-based algorithm

When considering classical (i.e. first order term) periodic homogenization under the small perturbations assumption, the problem to solve consists of a heterogeneous unit cell $\Omega$ described by its stiffness tensor field $c$, submitted to an average displacement (first order) gradient $\nabla U$ with periodic boundary conditions

$$ \left\{ \begin{array}{l} \mathrm{div}(\sigma(\mathbf{x})) = 0 \\ \sigma(\mathbf{x}) = c(\mathbf{x}) : \varepsilon(\mathbf{x}) \\ \varepsilon(\mathbf{x}) = (\nabla \mathbf{u})^{\mathrm{sym}}(\mathbf{x}) \\ \mathbf{u}(\mathbf{x}) = \nabla U \cdot \mathbf{x} + \tilde{\mathbf{u}}(\mathbf{x}) \end{array} \right. \quad \text{with} \quad \left\{ \begin{array}{l} \tilde{\mathbf{u}} \text{ periodic on } \partial\Omega \\ \sigma \cdot \mathbf{n} \text{ antiperiodic on } \partial\Omega \end{array} \right. \qquad (3) $$

The displacement field is written as the sum of the prescribed linear term and the unknown fluctuation field $\tilde{u}$. Using Equation (1) and Equation (2), Problem (3) written on $\Omega$, is equivalent to

$$ \left\{ \begin{array}{l} \mathrm{div}(\sigma) = 0 \\ \sigma = c : \varepsilon \quad \text{with} \quad \left\{ \begin{array}{l} \tilde{\varepsilon} \text{ compatible} \\ \tilde{\varepsilon} \cdot \Omega\text{-periodic} \\ \tilde{\sigma} \cdot \Omega\text{-periodic} \end{array} \right. \\ \varepsilon = E + \tilde{\varepsilon} \end{array} \right. \qquad (4) $$
where the x dependence is removed for the sake of concision, the average strain E being equal to the symmetrized average gradient E = ∇U<sup>sym</sup>. Also, ̃ε being defined as a fluctuation term, its average vanishes, 〈 ̃ε 〉 = 0. Introducing a homogeneous reference stiffness tensor c<sub>0</sub> allows to equivalently define

$$
\left\{
\begin{array}{l}
\mathbf{div}(c_0 : \varepsilon + \tau) = 0 \\
\tau = c : \varepsilon - c_0 : \tilde{\varepsilon} \\
\varepsilon = E + \tilde{\varepsilon}
\end{array}
\right.
\quad \text{with} \quad
\left\{
\begin{array}{l}
\tilde{\varepsilon} \text{ compatible} \\
\tilde{\varepsilon} \Omega\text{-periodic} \\
\sigma = c_0 : \tilde{\varepsilon} + \tau \Omega\text{-periodic}
\end{array}
\right.
\qquad (5)
$$

Temporarily assuming the so-called polarization τ as a known Ω-periodic field yields the auxiliary problem

$$
\tilde{\epsilon} = -\Gamma_0 * \tau \quad (\leftrightarrow \text{Fourier}) \quad \Leftrightarrow \quad \left\{
\begin{array}{l}
\mathbf{div}(c_0 : \tilde{\epsilon} + \tau) = 0 \\
\tilde{\epsilon} \text{ compatible} \\
\tilde{\epsilon} \Omega\text{-periodic} \\
\sigma = c_0 : \tilde{\epsilon} + \tau \Omega\text{-periodic}
\end{array}
\right.
\qquad (6)
$$

which can be solved by the application of the Green operator Γ₀, straightforward in Fourier space as the stiffness c₀ is homogeneous (see the Appendix of (Moulinec and Suquet 1998) for a detailed description).

Finally, the initial Problem (4), equivalent to Problem (5), reduces to Problem (6) along with

$$
\tau = c : \varepsilon - c_0 : \tilde{\varepsilon}, \tag{7a}
$$

$$
\tilde{\epsilon} = -\Gamma_0 * \tau \quad (\text{Fourier}), \tag{7b}
$$

$$
\epsilon = E + \tilde{\epsilon}. \tag{7c}
$$

Note that periodicity and compatibility conditions are gathered in Problem (6) and automatically
fulfilled when applying the Green operator in Fourier space, so that Problem (7) gathers a set
of equations, which, followed step by step, defines the simple fixed-point algorithm proposed
initially by Moulinec and Suquet (1998). Initializing the strain field ε (with ε = E) allows to
evaluate the polarization field τ in Problem (7a). Then, the Green operator applies (in Fourier
space) on the polarization to evaluate the fluctuation strain field ̃ε, Problem (7b), which, added to
the average (applied) strain field provides the strain field ε, Problem (7c) (In practice, it can be
done easily in Fourier space on the null frequency). This strain field allows to begin a second
iterate and so on until convergence.

2.3 A straightforward extension to strain gradient loadings

The purpose of this section is to add a compatible non-uniform applied strain ε* to the homoge-
neous applied strain E, modifying the initial Problem (4) in Problem (8):

$$
\left\{
\begin{array}{l}
\mathbf{div}(\sigma) = 0 \\
\sigma = c : \varepsilon \\
\varepsilon = E + \varepsilon^* + \tilde{\varepsilon}
\end{array}
\right.
\quad \text{with} \quad
\left\{
\begin{array}{l}
\tilde{\varepsilon} \text{ compatible} \\
\tilde{\varepsilon} \Omega\text{-periodic} \\
\tilde{\sigma} \Omega\text{-periodic}
\end{array}
\right.
\qquad (8)
$$

Following the different steps introduced in Section 2.2, Problem (8) is equivalent to

$$
\tau = c : \varepsilon - c_0 : \tilde{\varepsilon} \tag{9a}
$$

$$
\tilde{\epsilon} = -\Gamma_0 * \tau \quad (\text{Fourier}) \tag{9b}
$$

$$
\epsilon = E + \epsilon^* + \tilde{\epsilon} \tag{9c}
$$

It depicts a fixed-point algorithm (succession of steps (a) → (b) → (c) → (a) → (b) ...). The only simple difference with Equation (7) is step (c) with the addition of the applied strain field ε*.

Among the compatible non-uniform applied strain fields ε*, the simplest ones are linear fields. Actually, the compatibility condition that involves double derivatives of the strain components is satisfied as the double derivatives automatically vanish for linear fields. In the following, the heterogeneous applied strain fields are limited to such linear fields characterized by their strain gradients G as ε*(x) = G · x or equivalently ε*<sub>ij</sub> = G<sub>ijk</sub>x<sub>k</sub>, assuming x = 0 is the center of the unit-cell. Note that G is symmetric with respect to its first two indices (due to the symmetry of the strain tensor) so that it consists of 6 × 3 = 18 components.
## 2.4 Restrictions on the applied strain gradient

It must be emphasized that among the 18 components of the applied strain gradient only 9 can really be applied in this context of periodic boundary conditions. Mathematical foundations can be found in Appendix A. For any $G_{ijj}$ applied strain gradient, a periodic and compatible fluctuation displacement component $\tilde{u}_i$ can be found that modifies its effect. For example, a fluctuation displacement $\tilde{u}_1 = \tilde{G}_{122}x_2^2$ can modify the effect of the applied strain gradient $G_{122}$. Actually, $\tilde{u}_1$ satisfies the periodicity condition (i.e $\tilde{u}_1(h_2/2) = \tilde{u}_1(-h_2/2)$, with $h_2$ the cell size in direction 2) and it is associated to a fluctuation strain field $\tilde{\epsilon}_{12} = \tilde{G}_{122}x_2$, so that the average strain gradient seen by the unit-cell is not $G_{122}$ but $G_{122} + \tilde{G}_{122}$. Hence, these strain gradient components $ijj$ cannot be prescribed. Note that this property has been validated in our FFT-based implementation, see Equation (9), on a simple homogeneous unit-cell: the resulting stress and strain fields remain homogeneous even when introducing the strain gradient components $G_{ijj}$. The nine remaining components are associated to six bending loadings (components $G_{ijj}$ with $i \neq j$) and three torsion loading (components $G_{ijk}$ with $i \neq j, i \neq k$, and $j \neq k$).

## 2.5 Implementation considerations

Considering an existing FFT-based code able to solve Problem (7), transforming this code to solve Problem (9) is obvious. It consists in adding the applied strain gradient field when adding the homogeneous applied strain field to the periodic strain fluctuation. This is done, first when initializing the strain field (with a null fluctuation) and then, at the end of each iteration. The impact on the memory footprint is null, as no additional field needs to be allocated, and the additional computational cost is low compared to the complete cost of an iteration.

The implementation is detailed in Algorithm 1. Note that in the loop, $\tilde{\epsilon}$ stores the strain

**Algorithm 1** Standard FFT-based algorithm with strain gradient: linear behaviour.

1 Initialization: $\epsilon \leftarrow E + \epsilon^*$, $\tilde{\epsilon} \leftarrow 0$

2 **for** *i* **do**

3 $\sigma \leftarrow c : \epsilon$

4 $\hat{\sigma} \leftarrow \text{FFT}(\sigma)$

5 $\hat{\tilde{\epsilon}} \leftarrow -\Gamma_0 : (\hat{\sigma} - c_0 : \hat{\tilde{\epsilon}})$

6 $\tilde{\epsilon}(0) \leftarrow E$

7 $\tilde{\epsilon} \leftarrow \text{iFFT}(\tilde{\epsilon})$

8 $\epsilon \leftarrow \tilde{\epsilon} + \epsilon^*$

9 **end for**

fluctuation added to the average strain $E$ using $\hat{\tilde{\epsilon}}(0) = E$. In addition to this fixed-point algorithm, a convergence acceleration procedure is implemented in the code AMITEX_FFTP (Gélébart 2022). Its description can be found in (Chen et al. 2019a). It consists of storing $N$ couples of solutions $\tilde{\epsilon}^i$ and residual $\tilde{\epsilon}^i - \tilde{\epsilon}^{i-1}$ associated to the last four iterations, and, every $n$ iterations, suggests a solution deduced from a combination of the stored solutions (in AMITEX_FFTP, $N = 4$ and the default value for $n$ is 3).

## 2.6 Extension to non-linear behaviors

The extension to non-linear behaviors is rather straightforward. The evolution of the loading as a function of time, in terms of average strain $E$ and average strain gradient $G$, is discretized into time steps. The non-linear behavior of materials is classically introduced through internal variables $\alpha$ and the stress and internal variables between time steps $i$ and $i + 1$, associated to an applied strain increment $\tilde{\epsilon}^{i+1} - \tilde{\epsilon}^i$, is formally given by

$$[\sigma^{i+1}, \alpha^{i+1}] = \mathcal{F}([\sigma^i, \alpha^i, \epsilon^i, \epsilon^{i+1}]) \quad (10)$$

This generic description allows for a standard implementation of behaviors law. For example, these implementations are introduced in the Finite Element code ABAQUS through umat procedures. A similar standard is used in our FFT-based code AMITEX_FFTP. This standard allows for using the code generator MFRONT (Helfer et al. 2015) which greatly simplifies the numerical implementation, sometimes tedious, of Equation (10) that generally consists of the resolution of a non-linear system of differential equations.
A slight modification of Algorithm 1, dedicated to linear behaviours, is used to determine the fields $\sigma^{i+1}$, $\alpha^{i+1}$, $\varepsilon^{i+1}$ from $\sigma^i$, $\alpha^i$, $\varepsilon^i$ and the average applied load $E^{i+1}$, $\varepsilon^{*i+1}$ in Algorithm 2. The two minor adaptations to non-linearity are: (1) the linear behaviour $\sigma = c : \varepsilon$ is simply replaced by the non-linear behaviour in Equation (10), (2) the initialization step is replaced by a linear time extrapolation from the knowledge of the current and previous steps $i$ and $i-1$: $\varepsilon^{i+1} \leftarrow \varepsilon^i + (\varepsilon^i - \varepsilon^{i-1})(t^{i+1} - t^i)/(t^i - t^{i-1})$.

**Algorithm 2** Standard FFT-based algorithm with strain gradient: non-linear behaviour.

1 Initialization: $\varepsilon^{i+1} \leftarrow \varepsilon^i + (\varepsilon^i - \varepsilon^{i-1})(t^{i+1} - t^i)/(t^i - t^{i-1})$, $\hat{\varepsilon} \leftarrow 0$
2 **for** $i$ **do**
3 $[\sigma^{i+1}, \alpha^{i+1}] \leftarrow \mathcal{F}([\sigma^i, \alpha^i, \varepsilon^i, \varepsilon^{i+1}])$
4 $\hat{\sigma}^{i+1} \leftarrow \text{FFT}(\sigma^{i+1})$
5 $\hat{\varepsilon} \leftarrow -\Gamma_0 : (\hat{\sigma}^{i+1} - c_0 : \hat{\varepsilon})$
6 $\tilde{\varepsilon}(0) \leftarrow E^{i+1}$
7 $\tilde{\varepsilon} \leftarrow \text{iFFT}(\tilde{\varepsilon})$
8 $\varepsilon^{i+1} \leftarrow \tilde{\varepsilon} + \varepsilon^{*i+1}$
9 **end for**

# 3 Application to the homogenization of beams and plates

The purpose of this section is to demonstrate that the slight modification proposed in the previous section, applied to an existing FFT-based code that is able to simulate microstructures with null elastic properties, is enough to consider the problem of thin plates or elongated beams homogenization, considering enlarged unit-cells with zero stiffness material to account for traction-free surfaces. Nguyen et al. (2008) has initially addressed the question of thin plate homogenization with an FFT-based solver, but his approach relied on a dedicated Green operator, limited to plates with plane free surfaces, whereas our approach relies on the standard Green operator and could be applied as well to corrugated plates. The approach requires an FFT-based solver able to account for null elastic properties. For the sake of validation and comparison, Section 3.1 reproduces all the simulation cases proposed in (Nguyen et al. 2008). Then, Section 3.2 extends the approach to the homogenization of beams, treated for the first time, to the author's knowledge, with an FFT-based solver. Finally, Section 3.3 extends the applications to non-linear behaviors, considering a heterogeneous beam with an elasto-plastic matrix.

## 3.1 Homogenization of plates

The purpose of this section is to demonstrate that a general FFT-based code, equipped with the minor modification proposed in this study, extends its application domain to the homogenization of plates. Validation comes from comparisons with analytical and numerical results obtained by Nguyen et al. (2008). Note that an extensive and refined comparison between our simulations and Nguyen's simulations is not easy because of the differences between the iterative algorithms, the discrete Green operators and the convergence criterion. Especially, the convergence criterion used by Nguyen is a 'macroscopic' criterion based on the 'relative' difference of the strain energy between two iterations whereas the criterion used in AMITEX_FFTP relies on the equilibrium condition of the 'local' stress field at the current iteration.

### 3.1.1 Method

**General setting** We consider a thin heterogeneous plate described by a unit cell $\Omega$ periodically repeated in the plane $(e_1, e_2)$, with direction $e_3$ normal to the plate. Assuming an applied strain field consistent with the Kirchhoff-Love modeling for plates, the problem to solve reads, according to Nguyen et al. (2008) citing (Caillerie and Nedelec 1984):

$$ \left\{ \begin{array}{l} \text{div}(\sigma) = 0 \\ \sigma = c : \varepsilon \\ \varepsilon = E + \chi x_3 + (\nabla \tilde{u})^{\text{sym}} \end{array} \right. \quad \text{with} \quad \left\{ \begin{array}{l} \tilde{u} \text{ periodic on } \partial_{12}\Omega \\ \sigma \cdot n \text{ antiperiodic on } \partial_{12}\Omega \\ \sigma \cdot n = 0 \text{ on } \partial_3\Omega \end{array} \right. \qquad (11) $$

where $\partial_{12}\Omega$ denotes the four faces of normal directions $e_1$ and $e_2$, submitted to periodic boundary conditions and $\partial_3\Omega$ the two faces submitted to a traction-free boundary condition. The prescribed
components of both the average strain E and strain gradient χ are limited to the in-plane components (see Appendix A for a detailed description of prescribed and 'measured' components) indicated in boldface:

$$E = \begin{pmatrix} E_{11} & E_{12} & E_{13} \\ E_{12} & E_{22} & E_{23} \\ E_{13} & E_{23} & E_{33} \end{pmatrix} \quad \text{and} \quad \chi = \begin{pmatrix} \chi_{11} & \chi_{12} & \chi_{13} \\ \chi_{12} & \chi_{22} & \chi_{23} \\ \chi_{13} & \chi_{23} & \chi_{33} \end{pmatrix} = \begin{pmatrix} G_{113} & G_{123} & G_{133} \\ G_{123} & G_{223} & G_{233} \\ G_{133} & G_{133} & G_{333} \end{pmatrix}. \quad (12)$$

Note that, due to mechanical couplings between directions, the out-of-plane components do not systematically vanish (i.e. even in the simplest case of an elastic, homogeneous and isotropic material, the components $E_{33}$ and $\chi_{33}$ arise from the classical Poisson effect). These out-of-plane components, that are not prescribed, can be evaluated as a post-treatment of the simulation. Note that Problem (11) is exactly the same as the one proposed in (Nguyen et al. 2008, Eqs 32 and 33) where it is mentioned that the out-of-plane components of E and χ are enforced to 0. The author's opinion is that prescribing these quantities is not consistent with the stress free boundary condition applied on the upper and lower plate surfaces (i.e. $\sigma \cdot n = 0$ on $\partial_3 \Omega$). For example, the transverse strain $E_{33}$ (but also $\chi_{33}$) arising when applying a uniaxial strain $E_{11}$ depends on the Poisson coefficient and cannot be set to 0. However, as observed in Section 3.2.1, the numerical results obtained with the present approach are consistent with the analytical and numerical results given in (Nguyen et al. 2008) for heterogeneous plates.

**FFT-based resolution on an enlarged unit-cell** In order to solve Problem (11), Nguyen et al. (2008) follow exactly the reasoning described in Section 2. However, by keeping the traction-free condition on $\partial_3 \Omega$, they have to modify the Green operator to account for this condition in the auxiliary problem.

In the present paper, instead of solving Problem (11) with a modified Green operator, we solve Problem (9) with the classical periodic Green operator, built with full periodic conditions, available in every FFT-based codes. In order to satisfy the traction-free boundary condition that appears in Problem (11) but not in Problem (9), the fully periodic problem is not solved on the unit-cell $\Omega$ but on an enlarged unit-cell $\Omega^*$ adding to $\Omega$ one layer with null elastic properties on each side of $\partial_3 \Omega$, see Figure 1. In practice the added layers are one voxel thickness in the

<figure><img src="image_2.png" /><figcaption>**Figure 1** Unit-cells used in simulations of heterogeneous plates: (a) laminate plate and (b) composite plate with cylinder (square base) inclusions. Additional void layers used in the simulations in orange.</figcaption></figure>

discretized problem and using one layer on each side allows to keep the same middle surface for $\Omega$ and $\Omega^*$. Due to the added void layers, the in-plane components of $E^*$ and $\chi^*$ applied to $\Omega^*$ are the same as applied to the sub-domain $\Omega$ but the out-of-plane components of $E^*$ and $\chi^*$, can be chosen arbitrarily. For the out-of-plane components of $\chi^*$, i.e. $G_{i33}^*$, it has been demonstrated in Section 2.4 and Appendix A that these components cannot be prescribed with periodic boundary conditions. Their value can be arbitrarily set to 0:

$$E^* = \begin{pmatrix} E_{11} & E_{12} & 0 \\ E_{12} & E_{22} & 0 \\ 0 & 0 & 0 \end{pmatrix} \quad \text{and} \quad \chi^* = \begin{pmatrix} \chi_{11} & \chi_{12} & 0 \\ \chi_{12} & \chi_{22} & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} G_{113} & G_{123} & 0 \\ G_{123} & G_{223} & 0 \\ 0 & 0 & 0 \end{pmatrix}. \quad (13)$$

In order to clarify this aspect, let us consider the in-plane component 22 and the out-of-plane component 33 successively. Let $\tilde{\epsilon}_{22}^*(x_1, x_3)$ be the linear average strain (component 22) in direction
2 so that $\forall \mathbf{x} \in \Omega^*$, $\bar{\epsilon}_{22}^*(x_1, x_3) = E_{22}^* + G_{223}^*x_3$ (the linear average contribution of strain fluctuation $\bar{\epsilon}_{22}$ vanishes). As, the linear average in direction 2 is performed in a single domain (whether $\Omega$ or in the void) then $\bar{\epsilon}_{22}(x_1, x_3)$, the linear average strain in domain $\Omega$, follows the same relation $\bar{\epsilon}_{22}(x_1, x_3) = \bar{\epsilon}_{22}^*(x_1, x_3) = E_{22} + G_{223}x_3$ with $E_{22} = E_{22}^*$ and $G_{223} = G_{223}^*$. In other words, the macroscopic components $E_{22}$ and $G_{223}$ applied to $\Omega$ are the same as the macroscopic components $E_{22}^*$ and $G_{223}^*$ applied to $\Omega^*$. The same idea can be followed with components 11 and 12.

Now, $\bar{\epsilon}_{33}^*(x_1, x_2)$ is the linear average strain (component 33) in direction 3 so that $\forall \mathbf{x} \in \Omega^*$, $\bar{\epsilon}_{33}^*(x_1, x_2) = E_{33}^*$. In that case, $\bar{\epsilon}_{33}(x_1, x_2)$, the linear average strain in domain $\Omega$, is different from $\bar{\epsilon}_{33}^*(x_1, x_2)$ as $\bar{\epsilon}_{33}^*(x_1, x_2) = f\bar{\epsilon}_{33}(x_1, x_2) + (1-f)\bar{\epsilon}_{33}^{\text{void}}(x_1, x_2)$, where $\bar{\epsilon}_{33}^{\text{void}}(x_1, x_2)$ represents the linear average strain in the void layers and $f$ the volume fraction of $\Omega$ in $\Omega^*$. As a consequence, $\bar{\epsilon}_{33}(x_1, x_2) \neq E_{33}^*$ and then $E_{33} \neq E_{33}^*$. In other words, the additional void layers are free to deform so that the average $E_{33}^*$ can be chosen arbitrarily. The same idea can be followed with components 13 and 23.

Finally, as the out-of-plane components of $\chi^*$ cannot be prescribed on $\Omega^*$ (see Section 2.4 and Appendix A), they cannot be prescribed on $\Omega$, a subset of $\Omega^*$. To make it short, in-plane components of $E^*$ and $\chi^*$, prescribed on $\Omega^*$, are also prescribed on $\Omega$. On the other hand, the out-of-plane components of $E$ and $\chi$ cannot be prescribed on $\Omega$ but can be evaluated as post-treatments.

**Homogenization of plates** Once, the problem solved at the micro-scale (on the heterogeneous unit-cell), the last step consists in defining appropriate averages for the macroscopic quantities used for modeling plates. For an homogenized material, considering $\sigma_{\alpha\beta}(\alpha)$, $\beta = 1$ or 2) as a function of $x_3$, the two macroscopic ‘stress’ tensors are the linear force tensor $N$, expressed in [N/m], and the linear moment tensor $M$, in [Nm/m], defined by:

$$N_{\alpha\beta} = \int_{-h}^{h} \sigma_{\alpha\beta} \, dx_3 \quad \text{and} \quad M_{\alpha\beta} = \int_{-h}^{h} x_3 \sigma_{\alpha\beta} \, dx_3. \qquad (14)$$

Now, considering a heterogeneous unit-cell (and $\sigma_{\alpha\beta}$ as a function of $x_1, x_2$ and $x_3$), the following averaging procedure is used, with $S$ the area of $S_{12}$ the cross-sections of $\Omega$ in plane (1, 2):

$$\begin{aligned} N_{\alpha\beta} &= \int_{-h}^{h} \left( \frac{1}{S} \int_{S_{12}}^{n} \sigma_{\alpha\beta} \, dS \right) dx_3 = \frac{1}{S} \int_{\Omega} \sigma_{\alpha\beta} \, dV \\ M_{\alpha\beta} &= \int_{-h}^{h} \left( \frac{1}{S} \int_{S_{12}}^{n} x_3 \sigma_{\alpha\beta} \, dS \right) dx_3 = \frac{1}{S} \int_{\Omega} x_3 \sigma_{\alpha\beta} \, dV \end{aligned} \qquad (15)$$

From the in-plane periodicity of $\tilde{u}$ and the equilibrium condition ($\int_{S_{12}} \sigma_{i3} \, dS = 0$ for all sections $S_{12}(x_3)$), it can be shown that the average of the microscopic energy is equal to the macroscopic energy. For plates, the energy surface density is considered and the so-called macro-homogeneity condition of Hill-Mandel reads

$$\frac{1}{S} \int_{\Omega} \sigma : \varepsilon \, dV = \frac{1}{S} \int_{\Omega} \sigma : (E + \chi x_3 + \tilde{\varepsilon}) \, dV = N : E^p + M : \chi^p. \qquad (16)$$

Finally, the homogenized elastic behavior relates the kinematic tensors ($E^p, \chi^p$), restrictions of the tensors ($E, \chi$) to their in-plane components, to the macroscopic ‘stress’ tensors ($N, M$) as

$$\begin{bmatrix} N \\ M \end{bmatrix} = \begin{bmatrix} A & B \\ B & D \end{bmatrix} \begin{bmatrix} E^p \\ \chi^p \end{bmatrix}. \qquad (17)$$

Superscript “P” is associated to in-plane components. Tensor $A$ describes the membrane behavior, $D$ the bending and torsion behaviors and $B$ the possible coupling arising, for example, with non-symmetric unit-cells with respect to the middle surface. From this definition, the energy surface density reads:

$$\frac{1}{S} \int_{\Omega} \sigma : \varepsilon \, dV = E^p \cdot A \cdot E^p + 2E^p \cdot B \cdot \chi^p + \chi^p \cdot D \cdot \chi^p. \qquad (18)$$
In practice, the identification of the homogenized behavior results from six simulations corresponding to six independent loadings ($E^P, \chi^P$) applied to the unit-cell. The mechanical approach post-treats the stress fields to evaluate the corresponding resultant stress tensors ($N, M$) with Equation (15) in order to identify the coefficients of Equation (17). The energetic approach post-treats the stress-strain fields to evaluate the corresponding energies in order to identify the coefficients of Equation (18). As soon as the macro-homogeneity Condition (16) holds, the energetic approach followed by Nguyen et al. (2008) is strictly equivalent to the mechanical approach followed in the present paper.

### 3.1.2 Laminate plate

Following the validation test used by Nguyen et al. (2008), we consider the case of a two phases laminate plate for which a closed form solution exists. The unit-cell $\Omega$ represented in Figure 1, is symmetric so that the coupling tensor $B$ vanishes. The volume fractions of phases are equal (50%), and their elastic properties (Young's modulus, Poisson's coefficient) are respectively (46 GPa, 0.3) and (10 GPa, 0.3) for phases 1 and 2. According to Nguyen et al. (2008), the spatial resolution used for the simulations are $2^n \times 2^n$ ($\times 1$) voxels (the code is three-dimensional) for the unit cell $\Omega$ and two additional layers enlarge the unit-cell to simulate traction-free boundary condition in $\Omega^*$, see Section 2. The non-vanishing components ($E_{22}, \chi_{22}$) of $E^P$ and $\chi^P$ define the two applied loadings: (1, 0) for the tensile loading and (0, 1) for the bending loading. They directly allow to evaluate $A_{22} = N_{22}/E_{22}$ and $D_{22} = M_{22}/\chi_{22}$. The corresponding loadings applied to the enlarged unit-cell $\Omega^*$ are:

$$ \text{tension: } E^* = \begin{pmatrix} 0 & 0 & 0 \\ 0 & E_{22} & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \chi^* = 0 \quad \text{and bending: } E^* = 0, \quad \chi^* = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \chi_{22} & 0 \\ 0 & 0 & 0 \end{pmatrix}. \quad (19) $$

The choice made for the reference material's Lamé coefficients ($\lambda_0, \mu_0$) is $\lambda_0 = (\lambda_{min} + \lambda_{max})/2$ with $\lambda_{min} = 0$ (for the void layers) and $\mu_0 = (\mu_{min} + \mu_{max})/2$ with $\mu_{min} = 0$ (for the void layers).

For the tensile loading, the algorithm reaches convergence after 6 iterations with a convergence criterion on the local equilibrium of $10^{-15}$ far below the threshold of $10^{-4}$ used in the simulation. It seems that the solver reaches the exact solution, up to the double precision round-off error.

Actually, the exact solution being constant per layer, it can be observed that the FFT-based method is able to find it if the discretization grid is parallel to the layers. From the macroscopic point of view, Table 1 compares the analytical results and numerical results from (Nguyen et al. 2008). The comparison with analytical result is exact up to the 8th digit. The comparison

<table><caption>**Table 1** Axial stiffness $A_{22}$ for laminate case: comparison between analytical results and two FFT-based methods.</caption><thead><tr><th>Resolution</th><th>Analytical</th><th>Present work</th><th>(Nguyen et al. 2008)</th></tr></thead><tbody><tr><td>8 × 8</td><td>3.07692307</td><td>3.07692303</td><td>3.077</td></tr><tr><td>4 × 4</td><td>3.07692307</td><td>3.07692303</td><td></td></tr></tbody></table>

with (Nguyen et al. 2008) is perfect up to the number of available digits. However, the converged value is obtained after ten iterations to converge on a macroscopic criterion ($10^{-4}$ on the average energy between two iterations), compared to 6 for this study on a more stringent criterion (local equilibrium). Finally, as the solution is constant per layer, it does not depend on the mesh resolution and exactly the same result is obtained with a resolution $4 \times 4$ (i.e. with one voxel per layer in the thickness).

For the bending loading, the algorithm reaches convergence after 15 iterations with, once again, a convergence criterion on the local equilibrium of $10^{-10}$ far below the threshold of $10^{-4}$. In the case of bending, the stress field is not constant per phase so that the stress field and consequently the resultant 'stress' tensor $M$, depends on the mesh size. Table 2 compares the analytical results with the numerical results in (Nguyen et al. 2008) and the present study, for various spatial resolutions. Both results are in excellent agreement with the analytical case and the error becomes very low as soon as resolution reaches $32 \times 32$. For resolution $16 \times 16$, results are still in good agreement (relative error ~ 0.26 %) and Nguyen's method provides a similar precision (a bit less but with a lower number of iterations, i.e. 10).
<table><caption>**Table 2** Bending stiffness $D_{22}$ for laminate case: comparison between analytical results and two FFT-based methods.</caption><thead><tr><th>Resolution</th><th>Analytical [GPa/m]</th><th>Present study</th><th>Rel. error</th><th>(Nguyen et al. 2008)</th><th>Rel. error</th></tr></thead><tbody><tr><td>16 × 16</td><td>0.0038003663</td><td>0.0037903502</td><td>26 × 10<sup>-4</sup></td><td>0.00382</td><td>52 × 10<sup>-4</sup></td></tr><tr><td>32 × 32</td><td>0.0038003663</td><td>0.0037978623</td><td>7 × 10<sup>-4</sup></td><td>0.003805</td><td>13 × 10<sup>-4</sup></td></tr><tr><td>64 × 64</td><td>0.0038003663</td><td>0.0037997403</td><td>2 × 10<sup>-4</sup></td><td>0.003802</td><td>5 × 10<sup>-4</sup></td></tr><tr><td>128 × 128</td><td>0.0038003663</td><td>0.0038002098</td><td>4 × 10<sup>-5</sup></td><td></td><td></td></tr><tr><td>256 × 256</td><td>0.0038003663</td><td>0.0038003272</td><td>1 × 10<sup>-5</sup></td><td></td><td></td></tr><tr><td>512 × 512</td><td>0.0038003663</td><td>0.0038003566</td><td>2 × 10<sup>-6</sup></td><td></td><td></td></tr><tr><td>1024 × 1024</td><td>0.0038003663</td><td>0.0038003638</td><td>7 × 10<sup>-7</sup></td><td></td><td></td></tr></tbody></table>

An additional example demonstrating that our FFT-based approach converges towards analytical results can be found in Appendix B with a Poisson ratio of $10^{-4}$. To conclude, the minor modification applied to our general FFT-based solver, allows reproducing the analytical solution obtained for a laminate plate, at least as precisely and as efficiently as an FFT-based solver dedicated to the homogenization of plates.

### 3.1.3 Plate with periodic fibers

This second validation test proposed in (Nguyen et al. 2008) considers a panel with a periodic distribution of parallelepiped fibers. Figure 1 describes the periodic unit-cell. The matrix and the inclusion have the same Poisson ratio 0.3, and the Young moduli of inclusion/matrix (in [GPa]) are (0/1), (10/1) and (1/1000). The AMITEX_FFTP code being able to apply whether the average stress or the average strain $E^*$ for each of the six components, applying a null average stress for the components 33, 13 and 23, was found to reduce the number of iterations at convergence of about 30% for the tensile loading. Hence, the loading used in this case is identical to the loading introduced for the laminate case in Equation (19), except for these three components, now set to null stress instead of null strain. The bending loading remains unchanged. The choice made for Lamé's coefficients follows the same rule as proposed in the laminate case. The spatial resolution used for the simulations are $2^n \times 2^n$ (×1) voxels for the unit cell $\Omega$ and two additional layers enlarge the unit-cell to simulate traction-free boundary condition in $\Omega^*$, see Section 2.

All the simulations reported in (Nguyen et al. 2008) have been reproduced and all the results are gathered in tables of Appendix C, in which the membrane stiffness $A_{22}$, bending stiffness $D_{22}$ and number of iterations are given as a function of the spatial resolution for the three elastic contrasts. The convergence criterion being very different between AMITEX_FFTP and (Nguyen et al. 2008), the comparison between iteration numbers at convergence is not significant. Hence, the evaluation of strain energy has been implemented in AMITEX_FFTP and the relative difference between two iterations is evaluated as a post-treatment to determine the number of iterations at convergence if using this criterion with a threshold of $10^{-4}$ used in (Nguyen et al. 2008).

If not useful to discuss all the tables one by one, their analysis allows drawing a few main conclusions. First, even for the coarsest resolution 16 × 16, with relative errors below 1%, results obtained with both FFT-based methods are very close to the reference solution obtained from finite element simulations in (Nguyen et al. 2008). Second, the results are also quite comparable to (Nguyen et al. 2008) with a relative error that decreases quite monotonously as a function of the spatial resolution. Third, the number of iterations obtained with our approach (associated to a field equilibrium criterion) remains reasonable with values comprised between 20 and 76. It is a bit lower for bending loadings than for tensile loadings, whereas the type of loading has no influence in (Nguyen et al. 2008). Finally, using the strain energy based criterion, the number of iterations at convergence are rather similar to (Nguyen et al. 2008).

Finally, additional results obtained with realistic material properties are also given in Appendix C. The material corresponds to a glass fiber/epoxy composite with the couple (Young's modulus, Poisson's coefficient) being respectively (73.1 GPa, 0.22) and (3.45 GPa, 0.35) for glass and matrix. The elastic contrast is around 20 and here, the Poisson coefficients are different. The homogenized coefficient exhibit a monotonous convergence and the number of iteration of the algorithm is a bit higher than the previous ones but remains below 100.
To conclude, the minor modification applied to our general FFT-based solver, allows repro-
ducing the numerical results obtained for a matrix/inclusion plate, almost as precisely and as
efficiently as a solver dedicated to the homogenization of plates.

3.2 Homogenization of beams

If the extension of the approach proposed for plates in (Nguyen et al. 2008) to the homogenization of beams is not straightforward especially with arbitrary cross-sections, it is direct with the present approach. After a description of the method proposed for the homogenization of beams, it is applied to different loadings with various geometries and elastic contrasts. As uniaxial tension is quite straightforward and has already been applied in previous studies (Chen et al. 2019a), the paper focuses on bending and torsion loadings.

This section considers two different beams with a square and a disk cross-section. In the first case the discretization is conform to the geometry while it is not in the second case. Cubic inclusions, with conform discretization are considered in both cases, so that the effect of a conforming discretization is limited to the shape of the cross-section which is of a prime interest in the present context.

Figure 2 represents the geometry of the unit-cell with a cube diameter of 50 mm, a unit-cell length of 100 mm, a square and disk diameter of 100 mm for the cross-section. Void volumes in blue complement the unit-cells Ω, with a minimal thickness of two voxels, to form the enlarged unit-cells Ω*. The center of the inclusion lies on the axis of the unit-cells. The spatial

<figure><img src="image_3.png" /><figcaption>**Figure 2** Rendering of the two unit-cells: on the left with a square cross-section, on the right with a disk cross-section. The cubic inclusion is in red, the added void volume is in blue and the matrix in transparency.</figcaption></figure>

resolutions, defined by the number of voxels on 100 mm, are 40, 80 and 160 (so that the unit-cells
are 40 × 44 × 44, 80 × 84 × 84 and 160 × 164 × 164). With these resolutions, the discretization of
the inclusion is conform to its geometry.

The matrix Young's modulus is 10 MPa and Poisson's ratio, 0.3 with three elastic behaviors for the inclusion: the matrix behavior, so that the beam is homogeneous, a void behavior leading to an infinite contrast and a behavior stiffer than the matrix with a high contrast of 1000.

In addition, the use of composite voxels (Gélébart and Ouaki 2015; Kabel et al. 2015; Kabel et al. 2017; Charière et al. 2020) is tested to improve the description of the geometry when using a non-conforming discretization, as observed with the disk cross-section. The behavior of the composite voxels, crossed by the interface between the matrix and the additional void volume, follows a homogenization rule. Among the three classical rules, the Reuss and laminate homogenization rules are not well-suited with void behavior: the first one simply replaces the voxel by a void voxel and the second leads to an infinitely anisotropic behavior which generates convergence issues in the FFT-based algorithm (observed in torsion loading). In addition, the Voigt homogenization rule, which assumes a homogeneous strain over the phases, is extremely simple to implement, even with non-linear behaviors. Voigt is used in the following. For the sake of simplicity, an approximate definition of volume fractions is used and reported in Appendix D.

The choice made for the reference material's Lamé coefficients ($\lambda_0, \mu_0$) is $X_0 = (X_{\min}+X_{\max})/2$, with $X$ equal to $\lambda$ or $\mu$. Note that the min and max are evaluated over the beam material, excluding the additional void volume in $\Omega^*$: for contrasts 0 and 1000, the differences are negligible, and for contrast 1, the reference material corresponds to the behavior of the homogeneous beam.
### 3.2.1 Method

**General settings** The unit-cell Ω with an arbitrary cross-section is periodic in direction **e**₁. Assuming an applied strain field consistent with the Bernoulli modeling for beams, the problem to solve reads:

$$ 
\left\{ 
\begin{array}{l} 
\mathrm{div}(\sigma) = 0 \\ 
\sigma = c : \varepsilon \\ 
\varepsilon = E + \chi_3 x_3 + \chi_2 x_2 + (\nabla \tilde{u})^{\mathrm{sym}} 
\end{array} 
\right. 
\quad \text{with} 
\quad 
\left\{ 
\begin{array}{l} 
\tilde{u} \text{ periodic on } \partial_1 \Omega \\ 
\sigma \cdot \mathbf{n} \text{ antiperiodic on } \partial_1 \Omega \\ 
\sigma \cdot \mathbf{n} = 0 \text{ on } \partial_{23} \Omega 
\end{array} 
\right. 
\qquad (20) 
$$

where $\partial_1 \Omega$ denotes the two faces of outward normal $\mathbf{e}_1$, subject to periodic boundary conditions, and $\partial_{23} \Omega$ the external contour of the beam subject to a traction-free boundary condition. The macroscopic applied strain is limited to components in boldface in

$$ E = \begin{pmatrix} E_{11} & E_{12} & E_{13} \\ E_{12} & E_{22} & E_{23} \\ E_{13} & E_{23} & E_{33} \end{pmatrix}, \quad \chi_3 = \begin{pmatrix} G_{113} & -\alpha & G_{133} \\ -\alpha & G_{223} & G_{233} \\ G_{133} & G_{233} & G_{333} \end{pmatrix} \quad \text{and} \quad \chi_2 = \begin{pmatrix} G_{112} & G_{122} & \alpha \\ G_{122} & G_{222} & G_{232} \\ \alpha & G_{232} & G_{332} \end{pmatrix}. \quad (21) $$

The component $E_{11}$ prescribes a tensile loading, the components $G_{112}$ and $G_{113}$ prescribe two bending loadings respectively, and $\alpha = G_{132} = -G_{123}$ prescribes the torsion loading. According to the arguments proposed in Appendix A.1 based on the vanishing components of E and G imposed by the periodicity condition in direction 1, only the component 11 of the E, $\chi_2$ and $\chi_3$ tensors can be prescribed. However, the components $G_{132} = -G_{123} = \alpha$ were successfully prescribed for the torsion loading in the simulations below. Arguments regarding the symmetry of the problem are given in Appendix A.1.

**FFT-based resolution on an enlarged unit-cell** In order to satisfy the traction-free boundary condition, the full periodic problem is solved on the enlarged unit-cell $\Omega^*$ with a void volume surrounding $\Omega$ in order to fulfill the traction free boundary condition on $\partial_{23}\Omega$. Due to the added void layers, the out-of-plane components prescribed to $\Omega^*$ are the same as prescribed to $\Omega$ (in bold) but the in plane components can be chosen arbitrarily. They are set to 0:

$$ E^* = \begin{pmatrix} E_{11} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \chi_3^* = \begin{pmatrix} G_{113} & -\alpha & 0 \\ -\alpha & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \quad \text{and} \quad \chi_2^* = \begin{pmatrix} G_{112} & 0 & \alpha \\ 0 & 0 & 0 \\ \alpha & 0 & 0 \end{pmatrix}. \qquad (22) $$

The discussion on prescribed components and components that can be chosen arbitrarily is similar to the discussion proposed for plates. It has been checked on a homogeneous beam with $E_{ij}^* = 1$ and $G_{223}^* = G_{332}^* = 1$, that the strain field is homogeneous with an average on $\Omega$ corresponding to $E_{11} = 1$, $E_{22} = E_{33} = -vE_{11}$ (and null shear strain components). The arbitrary choice of the components $G_{ij}^*$ is general, see Appendix A.1, and was already checked in Section 2.4. Finally, it has been verified that the torsion loading applied to $\Omega^*$ is actually prescribed on $\Omega$, as expected from similar arguments than those used for plates.

**Homogenization of beams** Once, the problem solved at the micro-scale (on the heterogeneous unit-cell), the last step consists in defining appropriate spatial averages for the macroscopic quantities used for modeling beams. Classically, integrations of the traction vector $\sigma \cdot \mathbf{e}_1$ over the cross-section provide the macroscopic force $\mathbf{N}$ (expressed in [N]) and moments $\mathbf{M}$ (in [N m]):

$$ \mathbf{N} = \int_{S_{23}} \sigma \cdot \mathbf{e}_1 \, d\mathbf{S} = \int_{S_{23}} \begin{pmatrix} \sigma_{11} \\ \sigma_{12} \\ \sigma_{13} \end{pmatrix} \, d\mathbf{S}, \quad \mathbf{M} = \int_{S_{23}} \begin{pmatrix} 0 \\ x_2 \\ x_3 \end{pmatrix} \wedge (\sigma \cdot \mathbf{e}_1) \, d\mathbf{S} = \int_{S_{23}} \begin{pmatrix} \sigma_{13}x_2 - \sigma_{12}x_3 \\ \sigma_{11}x_3 \\ -\sigma_{11}x_2 \end{pmatrix} \, d\mathbf{S}. \quad (23) $$

Component $N_1$ of the force vector is the axial force and the other components are shear forces. Component $M_1$ of the moment vector is the torsion moment and the other components are bending moments. For homogenization purpose, these quantities evaluated on a cross-section are averaged along $\mathbf{e}_1$ of the unit-cell $\Omega$ so that they can be evaluated from volume integrations over
the complete unit-cell $\Omega$:

$$ 
\begin{aligned}
N &= \frac{1}{\ell_1} \int_{L_1} \int_{S_{23}} \sigma \cdot e_1 \, dS \, dx_1 = \frac{1}{\ell_1} \int_{\Omega} \sigma \cdot e_1 \, dV, \\
M &= \frac{1}{\ell_1} \int_{L_1} \int_{S_{23}} \begin{pmatrix} 0 \\ x_2 \\ x_3 \end{pmatrix} \wedge \sigma \cdot e_1 \, dS \, dx_1 = \frac{1}{\ell_1} \int_{\Omega} \begin{pmatrix} 0 \\ x_2 \\ x_3 \end{pmatrix} \wedge \sigma \cdot e_1 \, dV.
\end{aligned}
\quad (24) $$

Using the axial periodicity of $\tilde{u}$ and the equilibrium condition

$$ \int_{S_{12}(x_3)} \sigma_{i3} \, dS = \int_{S_{13}(x_2)} \sigma_{i2} \, dS = 0 \quad (25) $$

for all longitudinal sections $S_{12}(x_3)$ and $S_{13}(x_2)$, it can be demonstrated that the average of the microscopic energy is equal to the macroscopic energy. In the case of beams, the energy linear density is considered and the so-called macro-homogeneity condition of Hill-Mandel reads

$$ \frac{1}{\ell_1} \int_{\Omega} \sigma : \varepsilon \, dV = \frac{1}{\ell_1} \int_{\Omega} \sigma : (E + \chi_2 x_2 + \chi_3 x_3 + \tilde{\varepsilon}) \, dV = N_1 E_{11} + M_2 G_{113} - M_3 G_{112} + 2M_1 \alpha. \quad (26) $$

The different terms are respectively associated to the uniaxial, the two bendings and the torsion loadings. The minus sign of the third term comes from the minus sign of $M_3$ in Equation (23). Finally, the homogenized elastic behavior relates the kinematic quantities $U = (E_{11}, \alpha, G_{112}, G_{113})$ to the forces and moments $F = (N_1, M_1, M_2, M_3)$ with the linear relationship $F = K \cdot U$.

As mentioned for plates, the Hill-Mandel condition allows using indifferently the energetic or the mechanical approach. The last one is used in the following: a single kinematic loading is applied and the stress field is post-treated to evaluate forces and moments and identify the stiffness components of $K$.

### 3.2.2 Beam bending

The bending loading is applied to the unit-cell $\Omega^*$ with $G_{112}^* = G_{112}$. Components $G_{222}^*$ and $G_{332}^*$ arbitrarily set to 0 in Equation (22) are now set to $-vG_{112}^*$. Even if components 222 and 332 are not prescribed to the unit-cell $\Omega$, they provide a better initialization of the strain field and reduce the number of iterations at convergence. For the homogeneous strain loading, three different propositions have been tested: all the average strain components equal to 0, all the average stress components equal to 0, average strain components $(1j)_{j=1,2,3}$ and average stress components $(ij)_{i=1,2,j=1,2}$ equal to zero. All the propositions provide the same bending moment $M_3$ (up to a relative difference of $\sim 10^{-4}$) and none of them outperforms the others regarding the number of iterations at convergence. Results below are presented with all the average strain components equal to 0. In practice, $G_{112}$ is set to 1 in the simulations so that the bending moment is directly the bending stiffness.

**Homogeneous beam bending** The case of a homogeneous beam is interesting, at first as a validation case for which the analytical solution is known, but also to focus on the question of the non-conforming discretization of the beam (for the disk cross-section) and the potential interest of using composite voxels (with Voigt homogenization rule in that case). This specific case is also interesting because the initialization of the strain field in the FFT-based algorithm directly provides the solution of the problem and convergence should achieve at the first iterate. Actually, the initial strain field is a null homogeneous strain field added to a linearly varying strain field associated to $G_{112}$ and $G_{222}^* = G_{332}^* = -vG_{112}$ (see introduction of Section 3.2.2). Hence, the initialization of the strain field imposes $\varepsilon_{12} = \varepsilon_{13} = \varepsilon_{23} = 0$ for shear strains, $\varepsilon_{11} = G_{112}x_2$ for the axial strain and $\varepsilon_{22} = \varepsilon_{33} = -v\varepsilon_{11}$ for the transverse strains. This strain field initialized on the whole unit-cell $\Omega^*$ also applies on the sub-volume $\Omega$ so that it is directly the solution of the bending problem. As a validation part, the implementation in the AMITEX_FFTP code confirms this property: convergence is reached at the first iterate.

The bending moments derived from the Euler-Bernoulli beam theory are:

$$ M_3 = E \frac{\pi D^4}{64} G_{112} \quad (\text{disk cross-section}) \quad M_3 = E \frac{D^4}{12} G_{112} \quad (\text{square cross-section}) \quad (27) $$
Figure 3 represents the evolution of the relative error of the bending moments as a function of the spatial resolution. For the square cross-section, with a conform discretization the error

<figure><img src="image_4.png" /><figcaption>**Figure 3** Relative error on the evaluation of the bending moment for homogeneous beams with square and disk cross-sections, and with the use of composite voxels for disk cross-sections.</figcaption></figure>

is very low even for the lowest resolution $6 \times 10^{-4}$ and decreases monotonously to around $10^{-5}$. On the contrary, for the disk cross-section, the error is much higher and decreases with the spatial resolution from 0.01 to reach a value around $10^{-4}$ for the highest resolution. This difference clearly exhibits, on the macroscopic property, the difference between conforming and non-conforming discretization. For the square cross-section, the error is purely associated to the numerical integration of a quadratic field ($\sigma_{11}x_2$) on a rectilinear grid, whereas for the disk cross-section, an additional and predominant error due to the non-conforming discretization of the cross-section is added. In that case, the composite voxels demonstrate their efficiency to improve the results: the relative error is between 5 and 100 times lower than without composite voxels. As expected, the factor for resolution 40 (factor 37) is higher than for resolution 160 (factor 5). The factor 100 for resolution 80 was not really expected. However, it must be kept in mind that the convergence criterion (normalized norm of the stress equilibrium gap) is $10^{-4}$ in the simulations, so that the errors below $10^{-4}$ should be considered with caution. Hence, the increase from $10^{-5}$ to $3 \times 10^{-5}$ for a resolution increasing from 80 to 160 should be considered as a stagnation. From the stress field point of view, as observed in Figure 4, the composite voxels only modify the stress within the voxels crossed by the interface. In that case, their effect is purely local (it will be different for torsion loading).

<figure><img src="image_5.png" /><figcaption>**Figure 4** Homogeneous beam: axial stress field obtained for the disk cross-section with (left) and without (right) composite voxels for the lowest spatial resolution 40. For visualization purposes, a quarter of the unit-cell is removed.</figcaption></figure>

**Heterogeneous beam bending** Figure 5 (right) represents the evolution of the bending moments (or bending stiffness as $G_{112}=1$) for the three different elastic contrasts. It is worth noting that in spite of high contrasts, the evolution of the bending stiffness is moderate. Different observations can explain this result. At first, the volume fraction of inclusions is not so high (12.5 % and 15.9 %, respectively for the square and disk cross-sections). Then, and probably more importantly, the inclusions are located at the center of the beam, which is the less stressed region of the beam. It is also noticeable that the difference on the moments between square and disk cross-sections is almost independent of the contrast, respectively (3.42, 3.42, 3.45) $\times 10^7$ for the
contrasts 0, 1, and 1000. In other words, the stiffening (or softening) induced by these stiff (or soft) inclusions, with respect to the homogeneous beam, corresponds to an additional positive (or negative) bending stiffness that is rather independent of the beam cross-section.

<figure><img src="image_6.png" /><figcaption>**Figure 5** Bending moments: (a) normalized by their converged values obtained for the highest resolution 160 with solid, dotted and dashed lines for contrasts 1, 0 and 1000, respectively; (b) converged values obtained at highest resolution: the arrows symbolize null contrasts in log axis. Note that ‘disk’ and ‘disk+composite voxels’ symbols are superimposed, confirming that the moments obtained with the ‘disk’ (without composite voxel) at resolution 160 are converged values.</figcaption></figure>

Regarding the convergence analysis with respect to the spatial resolution, Figure 5(a) displays the evolution of the bending moments, normalized by their converged value (obtained for resolution 160). The same conclusions as obtained for the homogeneous beam still hold for heterogeneous beams: convergence is very fast for the square cross-section, unlike the case of the disk cross-section without composite voxels, and using composite voxels significantly improves the convergence. The effect of the contrast, moderate on the converged values of the moments, in Figure 5(b), is also quite moderate on the convergence curves, in Figure 5(a).

Regarding the convergence analysis of the FFT-based algorithm (in our case, a fixed-point algorithm combined with a convergence acceleration procedure), presented in Table 3, four main conclusions can be drawn. First, the convergence is very efficient for void inclusions. Second,

<table><caption>**Table 3** Bending loading case: number of iterations at convergence for the various beams, contrasts, resolutions and use of composite voxels.</caption><thead><tr><th rowspan="2">Resolution</th><th colspan="3">Square</th><th colspan="3">Disk</th><th colspan="3">Disk with Composite Voxel</th></tr><tr><th>C = 0</th><th>C = 1</th><th>C = 1000</th><th>C = 0</th><th>C = 1</th><th>C = 1000</th><th>C = 0</th><th>C = 1</th><th>C = 1000</th></tr></thead><tbody><tr><td>40</td><td>16</td><td>1</td><td>456</td><td>16</td><td>1</td><td>423</td><td>25</td><td>1</td><td>699</td></tr><tr><td>80</td><td>21</td><td>1</td><td>600</td><td>15</td><td>1</td><td>336</td><td>24</td><td>1</td><td>666</td></tr><tr><td>160</td><td>28</td><td>1</td><td>759</td><td>15</td><td>1</td><td>366</td><td>24</td><td>1</td><td>763</td></tr></tbody></table>

stiff inclusions drastically deteriorate the convergence. Note that, in that case, the unit-cell Ω* exhibits two types of large elastic contrast: between the matrix and the infinitely soft voids and between the matrix and the 1000 times stiffer inclusions. Such high and opposite contrasts are not really well suited for FFT-based methods. However, convergence is achieved and the number of iteration remains reasonable. Third, the use of composite voxels deteriorates the convergence by a factor of around 2. Fourth, for the square cross-section, increasing the spatial resolution deteriorates the convergence, but it tends to stabilize in the worst case with the contrast of 1000. The understanding of these last two unexpected results is still under investigation.

The analysis of the local stress fields focuses on the simulations with an elastic contrast of 1000 that induces the highest contribution to the bending stiffness, compared to the homogeneous beam in Figure 5(b), and exhibits the worst algorithm convergence, see Table 3. The worse convergence of the algorithm is probably to be correlated with the spurious oscillations arising in and around the stiff inclusion. Comparing the disk and square cross-sections, the stress fields in and around the inclusion are almost identical, independent of the cross-section. This is consistent
with the previous observation, at the macroscopic scale, of an additional contribution to the bending stiffness independent of the cross-section. Finally, comparing the disk cross-section with and without composite voxels, the conclusion reported for homogeneous beam still holds in the heterogeneous cases: they essentially affect the voxels crossed by the interface, keeping unchanged the rest of the volume (it is not true for torsion loading).

Finally, Figure 6 displays the distribution of the axial stress $\sigma_{11}$ within the beam. The effect of composite voxels used for the disk cross-section is similar to the homogeneous case, see Figure 4. In addition, the stress distribution in the inclusion seems rather independent of the beam cross-section (disk or square).

<figure><img src="image_7.png" /><figcaption>**Figure 6** Axial stress field $\sigma_{11}$ for the square cross-section (bottom), and the disk cross-section with (top) and without (center) composite voxels, for the lowest spatial resolution 40 and a contrast of 1000 for the inclusion. For visualization purposes, a quarter of the unit-cell is removed. Right handside figures are zooms on the inclusion.</figcaption></figure>

### 3.2.3 Beam torsion

The torsion loading applies to the unit-cell $\Omega^*$ with $\alpha = G_{132} = -G_{123}$ and the other strain gradient components set to 0. For the homogeneous strain loading, three different propositions have been tested: all the average strain components equal to 0, all the average stress components equal to 0, average strain components $(1j)_{j=1,2,3}$ and average stress components $(ij)_{i=1,2j=1,2}$ equal to zero. All the propositions provide the same torsion moment $M_1$ (up to a relative difference of ~ $10^{-4}$) and none of the proposition outperforms the others regarding the number of iterations at convergence. Results below are presented with all the average strain components equal to 0. In
practice, $G_{132}$ is set to 1 in the simulations so that the torsion moment is directly the torsion stiffness.

**Homogeneous beam torsion** The torsion of homogeneous beams, submitted to bending in the previous section, is also interesting for validation purpose. However, when the initialization of the strain field solution directly provides the solution in the bending case, the torsion case requires iterating the algorithm for balancing the stress field. The torsion moments derived from the Euler-Bernoulli beam theory are

$$M_1 = \mu \frac{\pi D^4}{32} \alpha \text{ (disk cross-section)} \qquad M_1 = \mu D^4 \alpha \times 0.1406 \text{ (square cross-section).} \qquad (28)$$

For the square cross-section, the expression is deduced from numerical simulation or from Fourier series evaluated numerically (Bruhns 2003).

Figure 7 reports the evolution of the relative error of the torsion moment (or torsion stiffness) as a function of the spatial resolution. For the square cross-section, the error is very small

<figure><img src="image_8.png" /><figcaption>**Figure 7** Relative error on the evaluation of the torsion moment for homogeneous beams with square and disk cross-sections, and with the use of composite voxels for disk cross-sections.</figcaption></figure>

$5 \times 10^{-4}$ for the lowest resolution 40 and decreases up to $2 \times 10^{-4}$ for the highest resolution 160. The shear stress field $\sigma_{13}$, observed in Figure 8 (down), is not a simple linear function of $x_2$ but also depends on $x_3$, emphasizing that the solution, obtained after 9 iterations, is quite different from the initial stress field (linear function of $x_2$). In spite of an infinite contrast (with the surrounding void volume), the solution is perfectly smooth without any spurious oscillations.

On the other hand, for the disk cross-section, the relative error is much larger: decreasing from 1.4 % to 0.6 % between resolutions 40 and 160. The stress field, in Figure 8 (up), is almost a linear function of $x_2$ but exhibits important spurious oscillations in the neighborhood of the beam interface to accommodate the traction free boundary condition approximated on a regular grid. These important fluctuations are suspected to maintain a quite high error (0.6 %) for the

<figure><img src="image_9.png" /><figcaption>**Figure 8** Homogeneous beam under torsion loading: shear stress field $\sigma_{13}$ for disk cross-section without (left) and with composite voxels (center) and square cross-section (right). Spatial resolution 40. For visualization purposes, a quarter of the unit-cell is removed.</figcaption></figure>

highest resolution, whereas in the bending case, without any spurious oscillations, the error reduces to $10^{-4}$. Finally, the use of composite voxels is, once again, of a great help for improving the solution. Macroscopically, the relative error is very low, compared to the solution without
composite voxels: decreasing from $12 \times 10^{-4}$ up to $4 \times 10^{-4}$ between resolution 40 and 160. The effect is also very important on the stress field, observed in Figure 8: it drastically reduces the amount of spurious oscillations. Compared to the bending loading case for which the effect is purely local, here, composite voxels also affect the solution in their neighborhood. Surprisingly, in spite of a smoother solution, the number of iterations at convergence is higher with composite voxels, see Table 4. This unexpected result is still not well understood.

**Heterogeneous beam torsion** Figure 9(b) represents the evolution of the torsion moments (or torsion stiffness as $\alpha = 1$) for the three different contrasts. Similarly to the bending loading cases,

<figure><img src="image_10.png" /><figcaption>**Figure 9** Torsion moments: (a) normalized by their converged values for the highest resolution 160: solid, dotted and dashed lines respectively for contrasts 1, 0 and 1000; (b) converged values obtained at highest resolution: arrows symbolize infinite contrasts). Symbols ‘disk’ and ‘disk+composite voxels’ are superimposed confirming that the moments obtained with the ‘disk’ (without composite voxel) at resolution 160 are converged.</figcaption></figure>

the effect of inclusions on the stiffness is quite moderate, even if a bit more pronounced for the highest contrast 1000 with an increase from $1.08 \times 10^8$ N m to $1.41 \times 10^8$ N m. In addition, and in agreement with the bending loading case: the stiffening (or softening) induced by the stiff (or soft) inclusions, with respect to the homogeneous beam, corresponds to an additional positive (or negative) torsion stiffness that is rather independent of the beam cross-section.

Regarding the convergence analysis with respect to the spatial resolution shown in Figure 9(a), observations are similar to the bending case: convergence is better for square cross-section than for disk cross-section, and in that case, convergence is better with composite voxels.

The convergence analysis of the FFT-based algorithm (fixed-point algorithm combined with a convergence acceleration procedure), presented in Table 4, draws conclusions very similar to the bending loading case. For the sake of conciseness, the reader is invited to refer to the analysis of

<table><caption>**Table 4** Torsion loading case: number of iterations at convergence for the different beams, contrasts, resolutions and use of composite voxels.</caption><thead><tr><th rowspan="2">Resolution</th><th colspan="3">Square</th><th colspan="3">Disk</th><th colspan="3">Disk with Composite Voxel</th></tr><tr><th>C = 0</th><th>C = 1</th><th>C = 1000</th><th>C = 0</th><th>C = 1</th><th>C = 1000</th><th>C = 0</th><th>C = 1</th><th>C = 1000</th></tr></thead><tbody><tr><td>40</td><td>18</td><td>9</td><td>462</td><td>19</td><td>12</td><td>475</td><td>57</td><td>48</td><td>1236</td></tr><tr><td>80</td><td>19</td><td>12</td><td>453</td><td>21</td><td>12</td><td>648</td><td>57</td><td>51</td><td>1437</td></tr><tr><td>160</td><td>25</td><td>16</td><td>378</td><td>21</td><td>10</td><td>744</td><td>51</td><td>36</td><td>1350</td></tr></tbody></table>

Table 3 for bending. The main difference arises from the fourth conclusion: in the torsion case, for the square cross-section, increasing the spatial resolution do not deteriorates the convergence.

Finally, the analysis of the local stress fields focuses on the simulations with an elastic contrast of 1000 that induces the highest contribution to the torsion stiffness, compared to the homogeneous beam, see Figure 9(b), and exhibits the worst algorithm convergence (see Table 4). As for the bending loading case, the stress field in and around the inclusion (zooms in Figure 10) does not seem to be affected by the beam cross-section (square or disk), observation which is also consistent with the macroscopic observation of an additional contribution on the torsion stiffness
independent of the cross-section. If spurious oscillations almost disappear in torsion, these smoother solutions are not associated to a better convergence of the algorithm when comparing Table 3 with Table 4.

<figure><img src="image_11.png" /><figcaption>**Figure 10** Shear stress $\sigma_{13}$ obtained for the square cross-section (bottom), and the disk cross-section with (middle) and without (top) composite voxels, for the lowest spatial resolution (40) and a contrast of 1000 for the inclusion. For visualization purposes, a quarter of the unit-cell is removed. Right handside figures are zooms in the inclusion.</figcaption></figure>

### 3.3 Non-linear beam application

The aim of this final subsection is to demonstrate the ability of the method to account for non-linear behaviors. For that purpose, a beam with a square cross-section ($D \times D$ with $D = 100$ mm) and a cube inclusion ($d \times d \times d$ with $d = 75$ mm) is used. Note that, in order to emphasize the effect of the inclusion on the beam behavior, the inclusion size $d$ has been increased. Two spatial resolution are used, 40 and 80.

As a validation test, the material is considered homogeneous with an elasto-plastic behavior with perfect plasticity. Then, while the matrix remains elasto-plastic, the inclusion is now purely elastic. The elastic properties of the matrix are set to (10 MPa, 0.3) for the couple (Young's modulus, Poisson's coefficient). Time independent perfect plasticity is associated to a Von Mises criterion with a yield stress $R_0$, set to 0.1 MPa. The Poisson coefficient of the inclusion is 0.3 and three different Young's moduli are tested: 0 MPa (elastic contrast $C = 0$ for a void inclusion), 10 MPa (no elastic contrast with the matrix, $C = 1$), 200 MPa (elastic contrast $C = 20$ as observed
in glass/epoxy composites for example).

For bending tests, the applied load corresponds to a null average strain E and the average strain gradient component $G_{112}^* = G_{112}$ is set to $0.001 \text{ mm}^{-1}$ at the end of the loading. If the gradients $G_{222}$ and $G_{332}$, cannot be prescribed (see Appendix A), the choice of $G_{222}^*$ and $G_{332}^*$ influences the convergence of the iterative algorithm. The choice $G_{222}^* = G_{332}^* = -vG_{112}$ made in the previous section with linear behaviors is kept. Other strain gradient components $G_{ijk}^*$ are set to zero. For torsion tests, the non-zero components of the strain gradient are $G_{132}^* = G_{132}$ and $G_{123}^* = G_{123}$ with $G_{132} = -G_{123} = 0.001 \text{ mm}^{-1}$ at the end of the loading. Finally, in both tests, the load is discretized in 20 steps.

### 3.3.1 Homogeneous beam with perfect plasticity

For bending test, the bending moment can be evaluated analytically. It is the sum of two contributions. The first comes from the region located around the center of the beam (i.e $x_2 \in [-y_0, y_0]$) where the material behaves elastically. The second comes from the rest of the beam where the material flows plastically with $\sigma_{11} = R_0$. The coordinate $y_0$ that separates the two regions corresponds to the coordinate for which the stress in the elastic region reaches $R_0 = E\epsilon_0 = EG_{112}y_0$. If the loading is too low, the whole beam remains elastic and $y_0 = D/2$, so that the expression can be summarized by $y_0 = \min(D/2, R_0/(EG_{112}))$. Then, the evaluation of the bending moment from the two contributions reads

$$M = M_1 + M_2 = E \frac{D(2y_0)^3}{12} G_{112} + D R_0 \left(\frac{D^2}{4} - y_0^2\right) \quad \text{with} \quad y_0 = \min\left(\frac{D}{2}, \frac{R_0}{EG_{112}}\right). \quad (29)$$

A comparison of numerical and analytical results is given in Figure 11. The comparison exhibits an excellent agreement even with the lowest resolution.

A similar analysis can be conducted for the torsion loading applied to a beam with a disk cross-section. The analytical evaluation of the torsion moment reads

$$M = M_1 + M_2 = \mu G_{123} \pi r_0^4 + \frac{2\pi R_0}{\sqrt{3}} \left(\frac{D^3}{24} - \frac{r_0^3}{3}\right) \quad \text{with} \quad r_0 = \min\left(\frac{D}{2}, \frac{R_0/\sqrt{3}}{2\mu G_{112}}\right). \quad (30)$$

The comparison between analytical and numerical results displayed in Figure 11 reveals that without composite voxels, increasing the resolution from 40 to 80 improves the agreement and, as reported in the linear case, using composite voxels leads to a very nice agreement even with the lowest resolution.

<figure><img src="image_12.png" /><figcaption>Figure 11 (a) Bending moment for a square cross-section. (b) Torsion moment for a disk cross-section.</figcaption></figure>

For the bending test, the total number of iterations required by the algorithm was 129 and 164 for resolutions 40 and 80, respectively, or 6.5 and 8.2 iterations per step in average, which is quite low. In spite of the highly non-linear behavior of perfect plasticity, the algorithm remains very efficient. For the torsion test, the required iterations are 356 and 386 without composite voxels and 400 and 441 with composite voxels. The use of Voigt composite voxels appears very interesting as the increase of the number of iterations is balanced by an improved precision.
### 3.3.2 Heterogeneous beam with perfectly plastic matrix

Now the inclusion is elastic while the matrix remains elasto-plastic with perfect plasticity and the elastic contrasts are $C = 0$ (void inclusion), $C = 1$ (no elastic contrast) and $C = 20$ (similar to a glass/epoxy elastic contrast). The square cross-section is adapted to the regular mesh and composite voxels are not required (as in the previous section in the similar case).

Figure 12 represents the evolution of the bending and torsion moments as a function of the applied strain gradients. It is worth noting that the curves associated to the two resolutions (40

<figure><img src="image_13.png" /><figcaption>**Figure 12** Moments for a square cross-section and two resolutions *r*. Black curves correspond to homogeneous beams (no inclusion) obtained analytically for bending and numerically for torsion. Legend applies to both figures.</figcaption></figure>

and 80) are almost superimposed so that the numerical simulation is considered as converged with respect to the spatial resolution. Then, it is noticeable that bending and torsion results exhibit three similar tendencies:

* compared to the homogeneous beam, the inclusion with $C = 1$ gives the same elastic behavior, the same plastic yield (moment at the onset of plasticity), but a higher plastic limit (asymptotic moment),

* compared to the heterogeneous beam with $C = 1$, the stiff inclusion with $C = 20$ stiffens the elastic behavior, increases the yield stress but does not change the plastic limit,

* compared to the heterogeneous beam with $C = 1$, the porous inclusion with $C = 0$ softens the elastic behavior, decreases both the yield stress and the plastic limit.

As observed in the linear case, the total number of iterations, reported in Table 5, is affected by the elastic contrast. If the increase, compared to $C = 1$, is quite low for a void inclusion ($C = 0$),

<table><caption>**Table 5** Total number of iterations for 20 loading steps for homogeneous and heterogeneous beams (square cross-section, cubic inclusion) with various contrasts. The first number corresponds to resolution 40, the second to resolution 80.</caption><thead><tr><th></th><th>Homogeneous</th><th>C = 0</th><th>C = 1</th><th>C = 20</th></tr></thead><tbody><tr><th>Bending</th><td>129/164</td><td>278/324</td><td>248/296</td><td>1011/1279</td></tr><tr><th>Torsion</th><td>261/211</td><td>600/667</td><td>509/583</td><td>1655/1782</td></tr></tbody></table>

the increase is significant for the stiff inclusion ($C = 20$). However, it must be mentioned that, in spite of the high non-linearity (perfect-plasticity), even in the worst case of the torsion test with ($C = 20$), the computation times for simulations performed on a single machine with 14 cores, are less than 2.5 minutes to obtain converged results with resolution 40, that is 77 440 voxels. In the better case, the flexion test with $C = 1$, the computation time to obtain converged results with resolution 40 is 0.5 minute.

## 4 Conclusion

The method described in the present paper extends the application domain of any existing FFT-based code, from homogeneous applied strain loading to strain gradient loadings, with a
minor and low invasive modification. Actually, the modification consists in adding the applied strain gradient field when adding the homogeneous applied strain field to the periodic strain fluctuation. This is achieved, first when initializing the strain field and then, at the end of each iteration. The impact on the memory footprint is null, as no additional field needs to be allocated. Among the 18 components of the strain gradient, it has been shown that only nine components can be really prescribed here with periodic boundary conditions. However, this subset has proven to be large enough to address the problems of heterogeneous beams and plates homogenization. As a practical result, the massively parallel code AMITEX_FFTP (Gélébart 2022) now offers a new functionality allowing the users to prescribe torsion and bending loadings to beam or plate heterogeneous unit-cells.

The approach proposed for the homogenization of beams and plates relies on this method to apply bending and torsion loadings to heterogeneous unit-cells enlarged with void voxels to account for traction-free boundary conditions. An FFT-based method dedicated to the homogenization of plates had previously been proposed in (Nguyen et al. 2008), and their simulations have been reproduced. These comparisons, together with analytical and Finite Element comparisons, validate the approach and demonstrate that a non-specific FFT-based solver, with a minor modification, can be almost as precise and efficient as a dedicated one.

The extension of the method proposed for plates in (Nguyen et al. 2008) to the homogenization of beams with any arbitrary cross-section does not seem straightforward. To the best of the author's knowledge, it was considered in the present study for the first time with an FFT-based solver. The approach is partially validated from simulations of homogeneous beams, with disk and square cross-sections, subject to bending and torsion. Simulations of heterogeneous beams with void or stiff inclusions (with a high elastic contrast of 1000) demonstrate the robustness of the solver. In addition, the use of composite voxels proves once again its ability to improve the numerical solution when voxels are crossed by an interface. Here, the interface is the beam boundary, and composite voxels consist of two phases: the void and the beam material. Applied to a beam with a disk cross-section, composite voxels with the Voigt homogenization rule improve the estimation of macroscopic moments (torsion, bending) and reduce spurious oscillations observed in torsion in the neighborhood of the beam boundary. Surprisingly, in spite of smoother fields, the use of composite voxels deteriorates the convergence of the algorithm. This point should be clarified. From the application point-of-view, it is noticeable that for the considered heterogeneities, located at the beam center, the additional stiffness induced by the presence of inclusions is quite moderate and independent of the beam cross-section.

As a final numerical experiment, the simulation of homogeneous and heterogeneous beams with an elasto-plastic behavior (perfect plasticity) demonstrates that the method can be used in a wide range of applications including non-linear analysis.

As future prospects, the approach proposed for the homogenization of beams and plates could be applied to more complex, yet of practical interest, microstructures such as stranded cables and corrugated beams or plates, as studied in (Cartraud and Messager 2006). Finally, another interesting application comes with the use of FFT-based solvers to simulate volume elements for which DIC or DVC (Digital Image or Volume Correlation) measurements have been performed (Chen et al. 2019a). Up to now, FFT-based solvers can apply the average strain measured experimentally. The present method will allow extending the loading to strain gradients components measured experimentally. It could be useful for example if the in-situ tensile test device induces spurious bending or torsion loadings.

# A Consequences of periodic fluctuations on average strain and strain gradient fluctuations

## A.1 Periodicity in one direction

Let us consider the periodicity condition of the fluctuation displacement $\tilde{u}$ on opposites faces with normal parallel to $e_1$ (i.e. $\tilde{u}$ periodic on $\partial_1 \Omega$) with the property $\tilde{u}(x_1 + h_1, x_2, x_3) = \tilde{u}(x_1, x_2, x_3)$, $\forall x_2, x_3$. Using $\bar{x}$, a linear average of $x$ along $e_1$, leads to $\tilde{u}_{1,1}(x_2, x_3) = \tilde{u}_{2,1}(x_2, x_3) = \tilde{u}_{3,1}(x_2, x_3) = 0$, $\forall x_2, x_3$. The consequences on the average strain fluctuations are that the conditions $\langle \tilde{u}_{1,1} \rangle =$
$\langle \tilde{u}_{2,1} \rangle = \langle \tilde{u}_{3,1} \rangle = 0$ without any constraint on the other components imply that $\langle \tilde{\varepsilon}_{11} \rangle = 0$ without any constraint on other components. Above, since $\tilde{\varepsilon}$ is the symmetric part of the gradient, shear components do not systematically vanish. It can be observed additionally that $\tilde{u}_{i,j}(x_1 + h_1, x_2, x_3) = \tilde{u}_{i,j}(x_1, x_2, x_3)$, $\forall x_2, x_3$ and $i = 1, 2, 3, j = 2, 3$. The vanishing components $ijk$ of $\langle \tilde{u}_{i,jk} \rangle$ arising from this periodicity are reported below (the second column is simply an inversion of the derivation order:

<table><tr><td>121</td><td>112</td></tr><tr><td>131</td><td>113</td></tr><tr><td>221</td><td>212</td></tr><tr><td>231</td><td>213</td></tr><tr><td>321</td><td>312</td></tr><tr><td>331</td><td>313</td></tr></table>

Above, since $\tilde{G}$ is symmetric with respect to the first two indices, the vanishing components of $\tilde{G}$, the average strain gradient fluctuations, are reported in bold character. Similar results can be obtained in the other two directions.

## A.2 Application to the full periodicity problem

Applying the periodicity conditions below in the three directions, the consequences on the average strain fluctuations are $\langle \tilde{\varepsilon}_{ij} \rangle = 0, \forall i, j$. The consequences on the average strain gradient fluctuations are:

$$ \begin{aligned} \tilde{G}_{221} &= \tilde{G}_{231} = \tilde{G}_{331} = 0 \\ \tilde{G}_{112} &= \tilde{G}_{132} = \tilde{G}_{332} = 0 \\ \tilde{G}_{113} &= \tilde{G}_{123} = \tilde{G}_{223} = 0 \end{aligned} \qquad (A.1) $$

with no constraint on other components. Hence, when using $\varepsilon = E + \varepsilon^* + \tilde{\varepsilon}$ in Problem (8):

* all components of E are prescribed for the average strain,

* only the nine components of the gradient of $\varepsilon^*$ corresponding to $\tilde{G}_{ijk} = 0$, are prescribed for the average strain gradient (components $G_{ijj}$ cannot be prescribed).

The other components of the strain gradient of $\varepsilon^*$ (i.e. $G_{ijj}^*$) can be set to arbitrary values.

## A.3 Application to the periodic beam problem

The problem involves a single periodicity condition along $e_1$. The conclusions of Appendix A.1 apply. When using $\varepsilon = E + \chi_3 x_3 + \chi_2 x_2 + \tilde{\varepsilon}$ in Problem (20):

* only the component $E_{11}$ is prescribed for the average strain.

* only the component $\chi_2$ of $\chi_2$ and $\chi_3$ are prescribed for the average strain gradient,

* components 221, 231, 331 of the average strain gradient are also implicitly prescribed to 0.

The other components can be set to any arbitrary values, they will be 'measured' quantities during post-treatment.

According to these arguments, components 12 of $\chi_3$ and 23 of $\chi_2$ cannot be prescribed. However, in the case of symmetric unit-cells, the axial component $\tilde{u}_1$ of the fluctuation displacement, is equal to 0 on the two opposite faces. In addition to $\langle \tilde{u}_{2,13} \rangle = \langle \tilde{u}_{3,12} \rangle = 0$ reported in Appendix A.1, we have $\langle \tilde{u}_{1,23} \rangle = \langle \tilde{u}_{1,32} \rangle = 0$, so that the components of the strain gradient fluctuation $\tilde{G}_{213}$ and $\tilde{G}_{312}$ vanish and the corresponding components of the average strain gradient can be prescribed. In the general case (no specific symmetry), we have $G_{213} = -\alpha + \frac{1}{2}\langle \tilde{u}_{1,23} \rangle$ and $G_{312} = \alpha + \frac{1}{2}\langle \tilde{u}_{1,32} \rangle$. If the quantity $\langle \tilde{u}_{1,23} \rangle = \langle \tilde{u}_{1,32} \rangle$ is not prescribed to 0 by periodic boundary conditions, it cannot be adjusted (by the numerical resolution) to vanish the applied load, as it is the case for components $\langle \tilde{u}_{i,jj} \rangle$ that is able to vanish the applied strain gradient $G_{ijj}$. Actually if $\langle \tilde{u}_{1,23} \rangle$ increases to reduce the magnitude of $G_{213}$, on the other hand, it will increase $G_{312}$.

## A.4 Application to the periodic plate

The problem involves two periodicity conditions along $e_1$ and $e_2$. In $\varepsilon = E + \chi x_3 + (\nabla \tilde{u})^{\text{sym}}$ in Problem (11):
* only the components 11, 22 and 12 of *E* are prescribed for the average strain,

* only the components 11, 22 and 12 of $\chi$ are prescribed for the average strain gradient,

* components 221, 231, 331, 112, 132, 332 of the average strain gradient are implicitly prescribed to 0.

The other components can be set to any arbitrary values. They will be 'measured' quantities evaluated as a post-treatment of the simulation results.

## B Laminate plate-additional case

Table B.1 shows that the method converges towards analytical results.

<table><caption>**Table B.1** Bending stiffness $D_{22}$ for the laminate case: comparison between analytical results and the proposed FFT-based method.</caption><thead><tr><th>Resolution</th><th>Analytical [GPa/m]</th><th>Present study</th><th>Rel. error</th></tr></thead><tbody><tr><td>16 × 16</td><td>19.949927685076197</td><td>19.897348729113546</td><td>-0.26%</td></tr><tr><td>32 × 32</td><td>19.949927685076197</td><td>19.936783319118270</td><td>-0.07%</td></tr><tr><td>64 × 64</td><td>19.949927685076197</td><td>19.946641214046725</td><td>-0.016%</td></tr><tr><td>128 × 128</td><td>19.949927685076197</td><td>19.949106561516157</td><td>-0.004%</td></tr></tbody></table>

## C Plate with periodic fibers: comparison with (Nguyen et al. 2008)

All the simulation cases reported in (Nguyen et al. 2008) have been reproduced with our approach and results are gathered in Table C.2 for various contrasts (contrast = $E_{inclusion}/E_{matrix}$ in [GPa]) and loading cases. The absolute values are respectively the tension and bending stiffness $A_{22} = N_{22}/E_{22}$ and $D_{22} = M_{22}/\chi_{22}$ in [GPa/m]. The column “Iter.” reports the number of iterations at convergence of our code AMITEX_FFTP, based on a field equilibrium criterion, while the column “Iter.*” indicates the number of iterations to reach the convergence criterion used in (Nguyen et al. 2008), with a relative difference of strain energy between two iterations less than $10^{-4}$.

<table><caption>**Table C.2** Comparison analysis. “Iter.” means “Iterations”. “Rel. err.” means “Relative error”. Starred columns refer to (Nguyen et al. 2008).</caption><thead><tr><th></th><th>Resolution</th><th>FEM*</th><th>Present study</th><th>Rel. err.</th><th>FFT*</th><th>Rel. err.</th><th>Iter.</th><th>Iter.*</th></tr></thead><tbody><tr><td rowspan="5">Contrast = 0</td><td>16 × 16</td><td>0.5840</td><td>0.5822</td><td>-0.3%</td><td>0.5826</td><td>-0.24%</td><td>28/8</td><td>8</td></tr><tr><td>32 × 32</td><td>0.5840</td><td>0.5834</td><td>-0.1%</td><td>0.5837</td><td>-0.05%</td><td>45/11</td><td>9</td></tr><tr><td>64 × 64</td><td>0.5840</td><td>0.5838</td><td>-0.07%</td><td>0.5839</td><td>-0.01%</td><td>44/14</td><td>9</td></tr><tr><td>128 × 128</td><td>0.5840</td><td>0.5838</td><td>-0.07%</td><td>0.5840</td><td>~ 0%</td><td>48/15</td><td>9</td></tr><tr><td>Bending</td><td>16 × 16</td><td>8.224 × 10<sup>-4</sup></td><td>8.204 × 10<sup>-4</sup></td><td>-0.24%</td><td>8.242 × 10<sup>-4</sup></td><td>0.22%</td><td>19/8</td><td>8</td></tr><tr><td rowspan="5">Contrast = 0.001</td><td>32 × 32</td><td>8.224 × 10<sup>-4</sup></td><td>8.219 × 10<sup>-4</sup></td><td>-0.06%</td><td>8.228 × 10<sup>-4</sup></td><td>0.05%</td><td>19/8</td><td>8</td></tr><tr><td>64 × 64</td><td>8.224 × 10<sup>-4</sup></td><td>8.223 × 10<sup>-4</sup></td><td>+0.015%</td><td>8.225 × 10<sup>-4</sup></td><td>0.01%</td><td>19/8</td><td>8</td></tr><tr><td>128 × 128</td><td>8.224 × 10<sup>-4</sup></td><td>8.224 × 10<sup>-4</sup></td><td>-0.004%</td><td>8.224 × 10<sup>-4</sup></td><td>~ 0%</td><td>21/8</td><td>8</td></tr><tr><td>Bending</td><td>16 × 16</td><td>58.495</td><td>58.329</td><td>-0.28%</td><td>58.363</td><td>-0.23%</td><td>28/8</td><td>8</td></tr><tr><td>32 × 32</td><td>58.495</td><td>58.455</td><td>-0.07%</td><td>58.472</td><td>-0.04%</td><td>39/11</td><td>9</td></tr><tr><td rowspan="5">Contrast = 10</td><td>64 × 64</td><td>58.495</td><td>58.484</td><td>-0.02%</td><td>58.502</td><td>+0.01%</td><td>45/14</td><td>9</td></tr><tr><td>128 × 128</td><td>58.495</td><td>58.487</td><td>-0.01%</td><td>58.506</td><td>+0.02%</td><td>48/15</td><td>9</td></tr><tr><td>Bending</td><td>16 × 16</td><td>0.08225</td><td>0.08206</td><td>-0.23%</td><td>0.08244</td><td>+0.23%</td><td>50/8</td><td>8</td></tr><tr><td>32 × 32</td><td>0.08225</td><td>0.08221</td><td>-0.05%</td><td>0.08230</td><td>+0.06%</td><td>35/8</td><td>8</td></tr><tr><td>64 × 64</td><td>0.08225</td><td>0.08224</td><td>-0.005%</td><td>0.08227</td><td>+0.02%</td><td>20/8</td><td>8</td></tr><tr><td rowspan="5">Contrast = 10</td><td>128 × 128</td><td>0.08225</td><td>0.08225</td><td>-0.006%</td><td>0.08226</td><td>+0.01%</td><td>20/8</td><td>8</td></tr><tr><td>Bending</td><td>16 × 16</td><td>0.16288</td><td>0.16267</td><td>-0.13%</td><td>0.16209</td><td>-0.49%</td><td>49/14</td><td>15</td></tr><tr><td>32 × 32</td><td>0.16288</td><td>0.16284</td><td>-0.02%</td><td>0.16257</td><td>-0.19%</td><td>39/17</td><td>16</td></tr><tr><td>64 × 64</td><td>0.16288</td><td>0.16287</td><td>-0.01%</td><td>0.16279</td><td>-0.06%</td><td>60/17</td><td>16</td></tr><tr><td>128 × 128</td><td>0.16288</td><td>0.16290</td><td>+0.01%</td><td>0.16287</td><td>-0.01%</td><td>76/17</td><td>16</td></tr><tr><td rowspan="5">Contrast = 10</td><td>Bending</td><td>16 × 16</td><td>1.0652 × 10<sup>-4</sup></td><td>1.0597 × 10<sup>-4</sup></td><td>-0.5%</td><td>1.0701 × 10<sup>-4</sup></td><td>+0.46%</td><td>40/17</td><td>18</td></tr><tr><td>32 × 32</td><td>1.0652 × 10<sup>-4</sup></td><td>1.0638 × 10<sup>-4</sup></td><td>-0.13%</td><td>1.0659 × 10<sup>-4</sup></td><td>+0.07%</td><td>42/17</td><td>18</td></tr><tr><td>64 × 64</td><td>1.0652 × 10<sup>-4</sup></td><td>1.0649 × 10<sup>-4</sup></td><td>-0.03%</td><td>1.0654 × 10<sup>-4</sup></td><td>+0.02%</td><td>51/17</td><td>18</td></tr><tr><td>128 × 128</td><td>1.0652 × 10<sup>-4</sup></td><td>1.0652 × 10<sup>-4</sup></td><td>-0.002%</td><td>1.0654 × 10<sup>-4</sup></td><td>+0.02%</td><td>45/17</td><td>18</td></tr></tbody></table>
<table><caption>**Table C.3** Glass Epoxy composite.</caption><thead><tr><th></th><th>Resolution</th><th>Present study</th><th>Iterations</th></tr></thead><tbody><tr><td rowspan="4">Tension</td><td>16 × 16</td><td>0.621 221 × 10<sup>9</sup></td><td>68/26</td></tr><tr><td>32 × 32</td><td>0.622 367 × 10<sup>9</sup></td><td>80/35</td></tr><tr><td>64 × 64</td><td>0.622 609 × 10<sup>9</sup></td><td>113/28</td></tr><tr><td>128 × 128</td><td>0.622 700 × 10<sup>9</sup></td><td>122/32</td></tr><tr><td rowspan="4">Bending</td><td>16 × 16</td><td>3.889 901 074 × 10<sup>5</sup></td><td>62/20</td></tr><tr><td>32 × 32</td><td>3.906 817 410 × 10<sup>5</sup></td><td>65/20</td></tr><tr><td>64 × 64</td><td>3.910 938 081 × 10<sup>5</sup></td><td>71/20</td></tr><tr><td>128 × 128</td><td>3.912 030 558 × 10<sup>5</sup></td><td>71/20</td></tr></tbody></table>

## D Characterization of composite voxels: beam with cylinder cross-section

The beam volume fraction (cylinder with disk cross-section) within composite voxels is approximated from the intersection area between two disks: the disk base of the cylinder and a disk of radius $R_{cv}$, radius of the encompassing circle, centered on the voxel. The expressions are as follows:

$$ 
\begin{aligned}
S &= S_1 + S_2 \\
S_1 &= R_0^2 \arccos\left(\frac{d_i}{R_0}\right) - d_i \sqrt{R_0^2 - d_i^2} \\
S_2 &= R_{cv}^2 \arccos\left(\frac{d_1 - d_i}{R_{cv}}\right) - (d_1 - d_i) \sqrt{R_{cv}^2 - (d_1 - d_i)^2} \\
d_i &= \frac{R_0^2 - R_{cv}^2 + d_1^2}{2d_1} \\
f &= \frac{S}{\pi R_{cv}^2}.
\end{aligned} 
\quad (D.1) $$

$d_1$ being the distance between the center of the beam section and the center of the voxel, and $R_0$, the radius of the disk cross-section. In practice, in order to avoid too low (or high) volume fractions, the detection of composite voxels uses a distance of 0.9345$R_{cv}$ which corresponds to a threshold of 1 % of volume fraction (i.e. below this value, a voxel is not considered as composite).

## References

Bertin, N. and L. Capolungo (2018). A FFT-based formulation for discrete dislocation dynamics in heterogeneous media. *Journal of Computational Physics* 355:366–384. [DOI].

Brenner, R. and J. Bravo-Castillero (2010). Response of multiferroic composites inferred from a fast-Fourier-transform-based numerical scheme. *Smart Materials and Structures* 19(11):115004. [DOI].

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

Bruhns, O. (2003). *Advanced Mechanics of Solids*. Springer. [DOI].

Caillerie, D. and J. C. Nedelec (1984). Thin elastic and periodic plates. *Mathematical Methods in the Applied Sciences* 6(1):159–191. [DOI], [OA].

Cartraud, P. and T. Messager (2006). Computational homogenization of periodic beam-like structures. *International Journal of Solids and Structures* 43(3-4):686–696. [DOI], [OA], [HAL].

Charière, R., A. Marano, and L. Gélébart (2020). Use of composite voxels in FFT based elastic simulations of hollow glass microspheres/polypropylene composites. *International Journal of Solids and Structures* 182-183:1–14. [DOI], [HAL].

Chateau, C., L. Gélébart, M. Bornert, and J. Crépin (2015). Micromechanical modeling of the elastic behavior of unidirectional CVI SiC/SiC composites. *International Journal of Solids and Structures* 58(0):322–334. [DOI], [OA].

Chen, Y., L. Gélébart, C. Chateau, M. Bornert, C. Sauder, and A. King (2019a). 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. *International Journal of Solids and Structures* 161:111–126. [DOI], [OA].

Chen, Y., D. Vasiukov, L. Gélébart, and C. H. Park (2019b). A FFT solver for variational phase-field modeling of brittle fracture. *Computer Methods in Applied Mechanics and Engineering* 349:167–190. [DOI], [HAL].

Dietrich, F., D. Merkert, and B. Simeon (2019). Derivation of Higher-Order Terms in FFT-Based Numerical Homogenization. *Numerical Mathematics and Advanced Applications: ENUMATH 2017* (University of Bergen, Norway, Sept. 25, 2017–Sept. 29, 2017). Vol. 126. Lecture Notes in Computational Science and Engineering. Springer, pp 289–297. [DOI], [ARXIV].

Geers, M., E. Coenen, and V. Kouznetsova (2007). Multi-scale computational homogenization of structured thin sheets. *Modelling and Simulation in Materials Science and Engineering* 15(4):S393–S404. [DOI].

Gélébart, L., C. Chateau, and M. Bornert (2009). Conditions aux limites mixtes normales. *19e Congrès Français de Mécanique* (Marseille, France, Aug. 24, 2009–Aug. 28, 2009). [HAL].

Gélébart, L. (2022). *AMITEX_FFTP*. URL: https://amitexfftp.github.io/AMITEX/index.html.

Gélébart, L. and R. Mondon-Cancel (2013). Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials. *Computational Materials Science* 77:430–439. [DOI].

Gélébart, L. and F. Ouaki (2015). Filtering material properties to improve FFT-based methods for numerical homogenization. *Journal of Computational Physics* 294:90–95. [DOI], [ARXIV].

Helfen, C. and S. Diebels (2014). Computational homogenisation of composite plates: Consideration of the thickness change with a modified projection strategy. *Computers & Mathematics with Applications* 67(5):1116–1129. [DOI], [OA].

Helfer, T., B. Michel, J.-M. Proix, M. Salvo, J. Sercombe, and M. Casella (2015). Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform. *Computers & Mathematics with Applications* 70(5):994–1023. [DOI], [OA].

Huet, C. (1990). Application of variational concepts to size effects in elastic heterogeneous bodies. *Journal of the Mechanics and Physics of Solids* 38(6):813–841. [DOI], [HAL].

Kabel, M., T. Böhlke, and M. Schneider (2014). Efficient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deformations. *Computational Mechanics* 54(6):1497–1514. [DOI], [HAL].

Kabel, M., A. Fink, and M. Schneider (2017). The composite voxel technique for inelastic problems. *Computer Methods in Applied Mechanics and Engineering* 322:396–418. [DOI].

Kabel, M., D. Merkert, and M. Schneider (2015). Use of composite voxels in FFT-based homogenization. *Computer Methods in Applied Mechanics and Engineering* 294:168–188. [DOI].

Kanit, T., S. Forest, I. Galliet, V. Mounoury, and D. Jeulin (2003). Determination of the size of the representative volume element for random composites: statistical and numerical approach. *International Journal of Solids and Structures* 40(13-14):3647–3679. [DOI].

Kochmann, J., S. Wulfinghoff, B. Svendsen, and S. Reese (2015). Phase-field modeling of martensitic phase transformations in polycrystals coupled with crystal plasticity: A spectral-based approach. *PAMM* 15(1):317–318. [DOI], [OA].

Lebensohn, R. and A. Needleman (2016). Numerical implementation of non-local polycrystal plasticity using fast Fourier transforms. *Journal of the Mechanics and Physics of Solids* 97:333–351. [DOI], [OA].

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

Nguyen, T.-K., K. Sab, and G. Bonnet (2008). Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates. *International Journal of Solids and Structures* 45(25-26):6518–6534. [DOI], [OA].

Schneider, M., D. Merkert, and M. Kabel (2017). FFT-based homogenization for microstructures discretized by linear hexahedral elements. *International Journal for Numerical Methods in Engineering* 109(10):1461–1489. [DOI].

Tran, T.-H., V. Monchiet, and G. Bonnet (2012). A micromechanics-based approach for the
derivation of constitutive elastic coefficients of strain-gradient media. *International Journal of Solids and Structures* 49(5):783-792. [DOI], [OA].

Willot, F. (2015). Fourier-based schemes for computing the mechanical response of composites with accurate local fields. *Comptes Rendus Mécanique* 343(3):232-245. [DOI], [OA].

Zeman, J., J. Vondřejc, J. Novák, and I. Marek (2010). Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients. *Journal of Computational Physics* 229(21):8065-8071. [DOI], [arXiv].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

**Supplementary Material** The code used to generate the results exposed in this paper is freely available at the address: https://amitexfftp.github.io/AMITEX/index.html.

**Acknowledgements** The author sincerely thanks the reviewers for their fruitful remarks.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.

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