Julien Réthoré

Lionel Gélébart

linear and non-linear behaviors

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

à 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.

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 [START_REF] Gélébart | Conditions aux limites mixtes normales. 19e Congrès Français de Mécanique[END_REF]). If the choice of periodic BC is obvious for periodic microstructures, it also provides the best estimate for random microstructures [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF][START_REF] Chateau | Micromechanical modeling of the elastic behavior of unidirectional CVI SiC/SiC composites[END_REF] compared to KUBC or SUBC, that Lionel Gélébart Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates provide upper and lower bounds, at least for linear material [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF]. 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 [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF].

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 [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF][START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF], algorithms replacing the initial fixed-point algorithm to improve convergence [START_REF] Zeman | Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients[END_REF][START_REF] Gélébart | Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials[END_REF][START_REF] Kabel | Efficient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deformations[END_REF]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 [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Gélébart | Filtering material properties to improve FFT-based methods for numerical homogenization[END_REF][START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF]. 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 [START_REF] Lebensohn | Numerical implementation of non-local polycrystal plasticity using fast Fourier transforms[END_REF] or damage phase field (Chen et al. 2019b), the coupling with metallurgical phase fields such as martensitic phase transformation [START_REF] Kochmann | Phase-field modeling of martensitic phase transformations in polycrystals coupled with crystal plasticity: A spectral-based approach[END_REF] or other physics such as magneto-electricity [START_REF] Brenner | Response of multiferroic composites inferred from a fast-Fourier-transform-based numerical scheme[END_REF], the coupling with Dislocation Dynamics (DD) codes [START_REF] Bertin | A FFT-based formulation for discrete dislocation dynamics in heterogeneous media[END_REF]. 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 [START_REF] Tran | A micromechanics-based approach for the Lionel Gélébart Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates derivation of constitutive elastic coefficients of strain-gradient media[END_REF][START_REF] Dietrich | Derivation of Higher-Order Terms in FFT-Based Numerical Homogenization[END_REF]. 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 [START_REF] Geers | Multi-scale computational homogenization of structured thin sheets[END_REF] for which strain gradients are associated to flexion and torsion loadings. For plates, if the classical numerical implementation relies on the finite element method [START_REF] Geers | Multi-scale computational homogenization of structured thin sheets[END_REF][START_REF] Helfen | Computational homogenisation of composite plates: Consideration of the thickness change with a modified projection strategy[END_REF], [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF] 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 [START_REF] Cartraud | Computational homogenization of periodic beam-like structures[END_REF], 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]. 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

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF], 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.

FFT-based method extension to strain gradient loadings 2.1 Periodicity condition

The description of the periodicity condition prescribed on a unit-cell Ω can be done equivalently from two different points of view. In the first case, fields are defined on Ω and boundary conditions have to be applied on its boundary 𝜕Ω. In the second case, fields are defined on an infinite space and Ω-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

div(𝜎 (𝒙)) = 0 𝜎 • 𝒏 antiperiodic on 𝜕Ω ⇔ div(𝜎 (𝒙)) = 0 𝜎 Ω-periodic (1)
where ł𝜎 • 𝒏 antiperiodic on 𝜕Ωž means that the traction vector 𝜎 • 𝒏 on opposite points of opposite faces of Ω are of opposite sign. The periodicity condition of a displacement field ũ in

ũ periodic on 𝜕Ω ⇔ ũ Ω-periodic ⇔        ε compatible ε Ω-periodic ⟨ ε⟩ = 0 (2)
is also equivalently formulated on Ω in terms of strain field ε, adding compatibility and null average conditions in that case. Below, ł ũ periodic on 𝜕Ωž means that the displacements ũ evaluated on opposite points of opposite faces are equal.

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 Ω described by its stiffness tensor field 𝑐, submitted to an average displacement (first order) gradient ∇𝑼 with periodic boundary conditions

           div(𝜎 (𝒙) = 0 𝜎 (𝒙) = 𝑐 (𝒙) : 𝜀 (𝒙) 𝜀 (𝒙) = (∇𝒖) sym (𝒙) 𝒖 (𝒙) = ∇𝑼 • 𝒙 + ũ (𝒙) with ũ periodic on 𝜕Ω 𝜎 • 𝒏 antiperiodic on 𝜕Ω (3)
The displacement field is written as the sum of the prescribed linear term and the unknown fluctuation field ũ. Using Equation (1) and Equation (2), Problem (3) written on Ω, is equivalent to

       div(𝜎) = 0 𝜎 = 𝑐 : 𝜀 𝜀 = 𝐸 + ε with        ε compatible ε Ω-periodic σ Ω-periodic (4)
Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

where the 𝒙 dependence is removed for the sake of concision, the average strain 𝐸 being equal to the symmetrized average gradient 𝐸 = ∇𝑼 sym . Also, ε being defined as a fluctuation term, its average vanishes, ⟨ ε⟩ = 0. Introducing a homogeneous reference stiffness tensor 𝑐 0 allows to equivalently define

       div(𝑐 0 : 𝜀 + 𝜏) = 0 𝜏 = 𝑐 : 𝜀 -𝑐 0 : ε 𝜀 = 𝐸 + ε with        ε compatible ε Ω-periodic 𝜎 = 𝑐 0 : ε + 𝜏 Ω-periodic (5)
Temporarily assuming the so-called polarization 𝜏 as a known Ω-periodic field yields the auxiliary problem

ε = -Γ 0 * 𝜏 (↔ Fourier) ⇔            div(𝑐 0 : ε + 𝜏) = 0 ε compatible ε Ω-periodic 𝜎 = 𝑐 0 : ε + 𝜏 Ω-periodic (6)
which can be solved by the application of the Green operator Γ 0 , straightforward in Fourier space as the stiffness 𝑐 0 is homogeneous (see the Appendix of [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] for a detailed description).

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

       𝜏 = 𝑐 : 𝜀 -𝑐 0 : ε, ε = -Γ 0 * 𝜏 (Fourier), 𝜀 = 𝐸 + ε. (7a) (7b) (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 [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. Initializing the strain field 𝜀 (with 𝜀 = 𝐸) 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.

A straightforward extension to strain gradient loadings

The purpose of this section is to add a compatible non-uniform applied strain 𝜀 * to the homogeneous applied strain 𝐸, modifying the initial Problem (4) in Problem (8):

       div(𝜎) = 0 𝜎 = 𝑐 : 𝜀 𝜀 = 𝐸 + 𝜀 * + ε with        ε compatible ε Ω-periodic σ Ω-periodic (8)
Following the different steps introduced in Section 2.2, Problem ( 8) is equivalent to

       𝜏 = 𝑐 : 𝜀 -𝑐 0 : ε ε = -Γ 0 * 𝜏 (Fourier) 𝜀 = 𝐸 + 𝜀 * + ε (9a) (9b) (9c)
It depicts a fixed-point algorithm (succession of steps (𝑎) → (𝑏) → (𝑐) → (𝑎) → (𝑏) . . .). The only simple difference with Equation ( 7) is step (𝑐) 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 𝐺 as 𝜀 * (𝒙) = 𝐺 • 𝒙 or equivalently 𝜀 * 𝑖 𝑗 = 𝐺 𝑖 𝑗𝑘 𝑥 𝑘 , assuming 𝒙 = 0 is the center of the unit-cell. Note that 𝐺 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.

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

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 𝐺 𝑖 𝑗 𝑗 applied strain gradient, a periodic and compatible fluctuation displacement component ũ𝑖 can be found that modifies its effect. For example, a fluctuation displacement ũ1 = G122 𝑥 2 2 can modify the effect of the applied strain gradient 𝐺 122 . Actually, ũ1 satisfies the periodicity condition (i.e ũ1 (ℎ 2 /2) = ũ1 (-ℎ 2 /2), with ℎ 2 the cell size in direction 2) and it is associated to a fluctuation strain field ε12 = G122 𝑥 2 , so that the average strain gradient seen by the unit-cell is not 𝐺 122 but 𝐺 122 + G122 . Hence, these strain gradient components 𝑖 𝑗 𝑗 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 𝐺 𝑖 𝑗 𝑗 . The nine remaining components are associated to six bending loadings (components 𝐺 𝑖𝑖 𝑗 with 𝑖 ≠ 𝑗) and three torsion loading (components 𝐺 𝑖 𝑗𝑘 with 𝑖 ≠ 𝑗, 𝑖 ≠ 𝑘, and 𝑗 ≠ 𝑘).

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, ε stores the strain Algorithm 1 Standard FFT-based algorithm with strain gradient: linear behaviour.

1 Initialization:

𝜀 ← 𝐸 + 𝜀 * , ε ← 0 2 for 𝑖 do 3 𝜎 ← 𝑐 : 𝜀 4 σ ← FFT(𝜎) 5 ε ← -Γ 0 : ( σ -𝑐 0 : ε) 6 ε (0) ← 𝐸 7 ε ← iFFT( ε) 8
𝜀 ← ε + 𝜀 * 9 end for fluctuation added to the average strain 𝐸 using ε (0) = 𝐸. In addition to this fixed-point algorithm, a convergence acceleration procedure is implemented in the code AMITEX_FFTP [START_REF] Gélébart | AMITEX_FFTP[END_REF]. Its description can be found in (Chen et al. 2019a). It consists of storing 𝑁 couples of solutions ε𝑖 and residual ε𝑖 -ε𝑖-1 associated to the last four iterations, and, every 𝑛 iterations, suggests a solution deduced from a combination of the stored solutions (in AMITEX_FFTP, 𝑁 = 4 and the default value for 𝑛 is 3).

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 𝐸 and average strain gradient 𝐺, is discretized into time steps. The non-linear behavior of materials is classically introduced through internal variables 𝜶 and the stress and internal variables between time steps 𝑖 and 𝑖 + 1, associated to an applied strain increment ε𝑖+1 -ε𝑖 , is formally given by

[𝜎 𝑖+1 , 𝜶 𝑖+1 ] = F ( [𝜎 𝑖 , 𝜶 𝑖 , 𝜀 𝑖 , 𝜀 𝑖+1 ]) (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 [START_REF] Helfer | Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform[END_REF]) 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.

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

A slight modification of Algorithm 1, dedicated to linear behaviours, is used to determine the fields 𝜎 𝑖+1 , 𝜶 𝑖+1 , 𝜀 𝑖+1 from 𝜎 𝑖 , 𝜶 𝑖 , 𝜀 𝑖 and the average applied load 𝐸 𝑖+1 , 𝜀 * 𝑖+1 in Algorithm 2. The two minor adaptations to non-linearity are: (1) the linear behaviour 𝜎 = 𝑐 : 𝜀 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 𝑖 and 𝑖 -1: 𝜀 𝑖+1 ← 𝜀 𝑖 + (𝜀 𝑖 -𝜀 𝑖-1 ) (𝑡 𝑖+1 -𝑡 𝑖 )/(𝑡 𝑖 -𝑡 𝑖-1 ).

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

1 Initialization:

𝜀 𝑖+1 ← 𝜀 𝑖 + (𝜀 𝑖 -𝜀 𝑖-1 )(𝑡 𝑖+1 -𝑡 𝑖 )/(𝑡 𝑖 -𝑡 𝑖-1 ), ε ← 0 2 for 𝑖 do 3 [𝜎 𝑖+1 , 𝜶 𝑖+1 ] ← F ([𝜎 𝑖 , 𝜶 𝑖 , 𝜀 𝑖 , 𝜀 𝑖+1 ]) 4 σ𝑖+1 ← FFT(𝜎 𝑖+1 ) 5 ε ← -Γ 0 : ( σ𝑖+1 -𝑐 0 : ε) 6 ε (0) ← 𝐸 𝑖+1 7 ε ← iFFT( ε) 8 𝜀 𝑖+1 ← ε + 𝜀 * 𝑖+1 9 end for

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. [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF] 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]). 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.

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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]. 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.

Method

General setting We consider a thin heterogeneous plate described by a unit cell Ω periodically repeated in the plane (𝒆 1 , 𝒆 2 ), with direction 𝒆 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF] citing [START_REF] Caillerie | Thin elastic and periodic plates[END_REF]:

       div(𝜎) = 0 𝜎 = 𝑐 : 𝜀 𝜀 = 𝐸 + 𝜒𝑥 3 + (∇ ũ) sym with        ũ periodic on 𝜕 12 Ω 𝜎 • 𝒏 antiperiodic on 𝜕 12 Ω 𝜎 • 𝒏 = 0 on 𝜕 3 Ω (11)
where 𝜕 12 Ω denotes the four faces of normal directions 𝒆 1 and 𝒆 2 , submitted to periodic boundary conditions and 𝜕 3 Ω the two faces submitted to a traction-free boundary condition. The prescribed

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 𝐸 33 and 𝜒 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 𝐸 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. 𝜎 • 𝒏 = 0 on 𝜕 3 Ω). For example, the transverse strain 𝐸 33 (but also 𝜒 33 ) arising when applying a uniaxial strain 𝐸 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF] for heterogeneous plates.

FFT-based resolution on an enlarged unit-cell In order to solve Problem (11), [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF] follow exactly the reasoning described in Section 2. However, by keeping the traction-free condition on 𝜕 3 Ω, 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 Ω but on an enlarged unit-cell Ω * adding to Ω one layer with null elastic properties on each side of 𝜕 3 Ω, see Figure 1. In practice the added layers are one voxel thickness in the discretized problem and using one layer on each side allows to keep the same middle surface for Ω and Ω * . Due to the added void layers, the in-plane components of 𝐸 * and 𝜒 * applied to Ω * are the same as applied to the sub-domain Ω but the out-of-plane components of 𝐸 * and 𝜒 * , can be chosen arbitrarily. For the out-of-plane components of 𝜒 * , i.e. 𝐺 * 𝑖33 , 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:

𝐸 * = 𝑬 11 𝑬 12 0 𝑬 12 𝑬 22 0 0 0 0 and 𝜒 * = 𝝌 11 𝝌 12 0 𝝌 12 𝝌 22 0 0 0 0 = 𝑮 113 𝑮 123 0 𝑮 123 𝑮 223 0 0 0 0 . ( 13
)
In order to clarify this aspect, let us consider the in-plane component 22 and the out-of-plane component 33 successively. Let ε * 22 (𝑥 1 , 𝑥 3 ) be the linear average strain (component 22) in direction

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

2 so that ∀𝒙 ∈ Ω * , ε * 22 (𝑥 1 , 𝑥 3 ) = 𝐸 * 22 + 𝐺 *
223 𝑥 3 (the linear average contribution of strain fluctuation ε22 vanishes). As, the linear average in direction 2 is performed in a single domain (whether Ω or in the void) then ε22 (𝑥 1 , 𝑥 3 ), the linear average strain in domain Ω, follows the same relation

ε22 (𝑥 1 , 𝑥 3 ) = ε * 22 (𝑥 1 , 𝑥 3 ) = 𝐸 22 + 𝐺 223 𝑥 3 with 𝐸 22 = 𝐸 * 22 and 𝐺 223 = 𝐺 * 223 .
In other words, the macroscopic components 𝐸 22 and 𝐺 223 applied to Ω are the same as the macroscopic components 𝐸 * 22 and 𝐺 * 223 applied to Ω * . The same idea can be followed with components 11 and 12. Now, ε * 33 (𝑥 1 , 𝑥 2 ) is the linear average strain (component 33) in direction 3 so that ∀𝒙 ∈ Ω * , ε * 33 (𝑥 1 , 𝑥 2 ) = 𝐸 * 33 . In that case, ε33 (𝑥 1 , 𝑥 2 ), the linear average strain in domain Ω, is different from

ε * 33 (𝑥 1 , 𝑥 2 ) as ε * 33 (𝑥 1 , 𝑥 2 ) = 𝑓 ε33 (𝑥 1 , 𝑥 2 ) + (1 -𝑓 ) εvoid 33 (𝑥 1 , 𝑥 2 )
, where εvoid 33 (𝑥 1 , 𝑥 2 ) represents the linear average strain in the void layers and 𝑓 the volume fraction of Ω in Ω * . As a consequence, ε33 (𝑥 1 , 𝑥 2 ) ≠ 𝐸 * 33 and then 𝐸 33 ≠ 𝐸 * 33 . In other words, the additional void layers are free to deform so that the average 𝐸 * 33 can be chosen arbitrarily. The same idea can be followed with components 13 and 23 .

Finally, as the out-of-plane components of 𝜒 * cannot be prescribed on Ω * (see Section 2.4 and Appendix A), they cannot be prescribed on Ω, a subset of Ω * . To make it short, in-plane components of 𝐸 * and 𝜒 * , prescribed on Ω * , are also prescribed on Ω. On the other hand, the out-of-plane components of 𝐸 and 𝜒 cannot be prescribed on Ω 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 𝜎 𝛼𝛽 (𝛼), 𝛽 = 1 or 2) as a function of 𝑥 3 , the two macroscopic 'stress' tensors are the linear force tensor 𝑁 , expressed in [N/m], and the linear moment tensor 𝑀, in [Nm/m], defined by:

𝑁 𝛼𝛽 = ∫ ℎ -ℎ 𝜎 𝛼𝛽 d𝑥 3 and 𝑀 𝛼𝛽 = ∫ ℎ -ℎ 𝑥 3 𝜎 𝛼𝛽 d𝑥 3 . (14)
Now, considering a heterogeneous unit-cell (and 𝜎 𝛼𝛽 as a function of 𝑥 1 , 𝑥 2 and 𝑥 3 ), the following averaging procedure is used, with 𝑆 the area of 𝑆 12 the cross-sections of Ω in plane (1, 2):

𝑁 𝛼𝛽 = ∫ ℎ -ℎ 1 𝑆 ∫ 𝑛 𝑆 12 𝜎 𝛼𝛽 d𝑆 d𝑥 3 = 1 𝑆 ∫ Ω 𝜎 𝛼𝛽 d𝑉 𝑀 𝛼 𝛽 = ∫ ℎ -ℎ 1 𝑆 ∫ 𝑆 12 𝑥 3 𝜎 𝛼 𝛽 d𝑆 d𝑥 3 = 1 𝑆 ∫ Ω 𝑥 3 𝜎 𝛼 𝛽 d𝑉 (15)
From the in-plane periodicity of ũ and the equilibrium condition ( ∫ 𝑆 12 𝜎 𝑖3 d𝑆 = 0 for all sections 𝑆 12 (𝑥 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

1 𝑆 ∫ Ω 𝜎 : 𝜀 d𝑉 = 1 𝑆 ∫ Ω 𝜎 : (𝐸 + 𝜒𝑥 3 + ε) d𝑉 = 𝑁 : 𝐸 𝑝 + 𝑀 : 𝜒 𝑝 . ( 16
)
Finally, the homogenized elastic behavior relates the kinematic tensors (𝐸 𝑃 , 𝜒 𝑃 ), restrictions of the tensors (𝐸, 𝜒) to their in-plane components, to the macroscopic 'stress' tensors (𝑁 , 𝑀) as

𝑁 𝑀 = 𝐴 𝐵 𝐵 𝐷 𝐸 𝑃 𝜒 𝑃 . ( 17
)
Superscript ł𝑃ž is associated to in-plane components. Tensor 𝐴 describes the membrane behavior, 𝐷 the bending and torsion behaviors and 𝐵 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:

1 𝑆 ∫ Ω 𝜎 : 𝜀 d𝑉 = 𝐸 𝑃 • 𝐴 • 𝐸 𝑃 + 2𝐸 𝑃 • 𝐵 • 𝜒 𝑃 + 𝜒 𝑃 • 𝐷 • 𝜒 𝑃 . ( 18
)
Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

In practice, the identification of the homogenized behavior results from six simulations corresponding to six independent loadings (𝐸 𝑃 , 𝜒 𝑃 ) applied to the unit-cell. The mechanical approach post-treats the stress fields to evaluate the corresponding resultant stress tensors (𝑁 , 𝑀) 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF] is strictly equivalent to the mechanical approach followed in the present paper.

Laminate plate

Following the validation test used by [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF], we consider the case of a two phases laminate plate for which a closed form solution exists. The unit-cell Ω represented in Figure 1, is symmetric so that the coupling tensor 𝐵 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF], the spatial resolution used for the simulations are 2 𝑛 × 2 𝑛 (×1) voxels (the code is three-dimensional) for the unit cell Ω and two additional layers enlarge the unit-cell to simulate traction-free boundary condition in Ω tension: 𝐸 * = 0 0 0 0 𝑬 22 0 0 0 0 , 𝜒 * = 0 and bending:

𝐸 * = 0, 𝜒 * = 0 0 0 0 𝝌 22 0 0 0 0 . ( 19
)
The choice made for the reference material's Lamé coefficients

(𝜆 0 , 𝜇 0 ) is 𝜆 0 = (𝜆 min + 𝜆 max )/2
with 𝜆 min = 0 (for the void layers) and 𝜇 0 = (𝜇 min + 𝜇 max )/2 with 𝜇 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]). The comparison with analytical result is exact up to the 8 th digit. The comparison with [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]) 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 × 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 𝑀, depends on the mesh size. Table 2 compares the analytical results with the numerical results in [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]) 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 × 32. For resolution 16 × 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 2 Bending stiffness 𝐷 22 for laminate case: comparison between analytical results and two FFT-based methods.

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.

Plate with periodic fibers

This second validation test proposed in [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]) 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 𝐸 * 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 𝑛 × 2 𝑛 (×1) voxels for the unit cell Ω and two additional layers enlarge the unit-cell to simulate traction-free boundary condition in Ω * , see Section 2.

All the simulations reported in [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]) have been reproduced and all the results are gathered in tables of Appendix C, in which the membrane stiffness 𝐴 22 , bending stiffness 𝐷 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF], 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF].

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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]. Second, the results are also quite comparable to [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]) 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]). Finally, using the strain energy based criterion, the number of iterations at convergence are rather similar to [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF].

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.

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

To conclude, the minor modification applied to our general FFT-based solver, allows reproducing the numerical results obtained for a matrix/inclusion plate, almost as precisely and as efficiently as a solver dedicated to the homogenization of plates.

Homogenization of beams

If the extension of the approach proposed for plates in [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF] 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 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.

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 [START_REF] Gélébart | Filtering material properties to improve FFT-based methods for numerical homogenization[END_REF][START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF][START_REF] Kabel | The composite voxel technique for inelastic problems[END_REF][START_REF] Charière | Use of composite voxels in FFT based elastic simulations of hollow glass microspheres/polypropylene composites[END_REF]) 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 (𝜆 0 , 𝜇 0 ) is 𝑋 0 = (𝑋 min +𝑋 max )/2, with 𝑋 equal to 𝜆 or 𝜇. Note that the min and max are evaluated over the beam material, excluding the additional void volume in Ω * : for contrasts 0 and 1000, the differences are negligible, and for contrast 1, the reference material corresponds to the behavior of the homogeneous beam.

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

Method

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

       div(𝜎) = 0 𝜎 = 𝑐 : 𝜀 𝜀 = 𝐸 + 𝜒 3 𝑥 3 + 𝜒 2 𝑥 2 + (∇ ũ) sym with        ũ periodic on 𝜕 1 Ω 𝜎 • 𝒏 antiperiodic on 𝜕 1 Ω 𝜎 • 𝒏 = 0 on 𝜕 23 Ω (20)
where 𝜕 1 Ω denotes the two faces of outward normal 𝒆 1 , subject to periodic boundary conditions, and 𝜕 23 Ω the external contour of the beam subject to a traction-free boundary condition. 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 Ω * with a void volume surrounding Ω in order to fulfill the traction free boundary condition on 𝜕 23 Ω. Due to the added void layers, the out-of-plane components prescribed to Ω * are the same as prescribed to Ω (in bold) but the in plane components can be chosen arbitrarily. They are set to 0:

𝐸 * =
𝑬 11 0 0 0 0 0 0 0 0

, 𝜒 * 3 = 𝑮 113 -𝛼 0 -𝛼 0 0 0 0 0 and 𝜒 * 2 = 𝑮 112 0 𝛼 0 0 0 𝛼 0 0 . ( 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 𝐸 * 𝑖 𝑗 = 1 and 𝐺 * 223 = 𝐺 * 332 = 1, that the strain field is homogeneous with an average on Ω corresponding to 𝐸 11 = 1, 𝐸 22 = 𝐸 33 = -𝑣𝐸 11 (and null shear strain components). The arbitrary choice of the components 𝐺 * 𝑖 𝑗 𝑗 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 Ω * is actually prescribed on Ω, 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 𝜎 • 𝒆 1 over the cross-section provide the macroscopic force 𝑵 (expressed in [N]) and moments 𝑴 (in [N m]):

𝑵 = ∫ 𝑆 23 𝜎 • 𝒆 1 d𝑆 = ∫ 𝑆 23 𝜎 11 𝜎 12 𝜎 13 d𝑆, 𝑴 = ∫ 𝑆 23 0 𝑥 2 𝑥 3 ∧ (𝜎 • 𝒆 1 ) d𝑆 = ∫ 𝑆 23 𝜎 13 𝑥 2 -𝜎 12 𝑥 3 𝜎 11 𝑥 3 -𝜎 11 𝑥 2 d𝑆. ( 23
)
Component 𝑁 1 of the force vector is the axial force and the other components are shear forces.

Component 𝑀 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 𝒆 1 of the unit-cell Ω so that they can be evaluated from volume integrations over

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates the complete unit-cell Ω:

𝑵 = 1 ℓ 1 ∫ 𝐿 1 ∫ 𝑆 23 𝜎 • 𝒆 1 d𝑆 d𝑥 1 = 1 ℓ 1 ∫ Ω 𝜎 • 𝒆 1 d𝑉 , 𝑴 = 1 ℓ 1 ∫ 𝐿 1 ∫ 𝑆 23 0 𝑥 2 𝑥 3 ∧ 𝜎 • 𝒆 1 d𝑆 d𝑥 1 = 1 ℓ 1 ∫ Ω 0 𝑥 2 𝑥 3 ∧ 𝜎 • 𝒆 1 d𝑉 . (24)
Using the axial periodicity of ũ and the equilibrium condition

∫ 𝑆 12 (𝑥 3 ) 𝜎 𝑖3 d𝑆 = ∫ 𝑆 13 (𝑥 2 ) 𝜎 𝑖2 d𝑆 = 0 (25)
for all longitudinal sections 𝑆 12 (𝑥 3 ) and 𝑆 13 (𝑥 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

1 ℓ 1 ∫ Ω 𝜎 : 𝜀 d𝑉 = 1 ℓ 1 ∫ Ω 𝜎 : (𝐸 + 𝜒 2 𝑥 2 + 𝜒 3 𝑥 3 + ε) d𝑉 = 𝑁 1 𝐸 11 + 𝑀 2 𝐺 113 -𝑀 3 𝐺 112 + 2𝑀 1 𝛼 . ( 26
)
The different terms are respectively associated to the uniaxial, the two bendings and the torsion loadings. The minus sign of third term comes from the minus sign of 𝑀 3 in Equation ( 23).

Finally, the homogenized elastic behavior relates the kinematic quantities 𝑼 = (𝐸 11 , 𝛼, 𝐺 112 , 𝐺 113 ) to the forces and moments 𝑭 = (𝑁 1 , 𝑀 1 , 𝑀 2 , 𝑀 3 ) with the linear relationship 𝑭 = 𝑲 • 𝑼 .

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 𝑲 .

Beam bending

The bending loading is applied to the unit-cell Ω * with 𝐺 * 112 = 𝐺 112 . Components 𝐺 * 222 and 𝐺 * 332 arbitrarily set to 0 in Equation ( 22) are now set to -𝑣𝐺 * 112 . Even if components 222 and 332 are not prescribed to the unit-cell Ω, 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 (1𝑗) 𝑗=1,2,3 and average stress components (𝑖 𝑗) 𝑖=1,2𝑗=1,2 equal to zero. All the propositions provide the same bending moment 𝑀 3 (up to a relative difference of ∼ 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, 𝐺 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 𝐺 112 and 𝐺 * 222 = 𝐺 * 332 = -𝑣𝐺 112 (see introduction of Section 3.2.2). Hence, the initialization of the strain field imposes 𝜀 12 = 𝜀 13 = 𝜀 23 = 0 for shear strains, 𝜀 11 = 𝐺 112 𝑥 2 for the axial strain and 𝜀 22 = 𝜀 33 = -𝑣𝜀 11 for the transverse strains. This strain field initialized on the whole unit-cell Ω * also applies on the sub-volume Ω 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:

𝑀 3 = 𝐸 𝜋𝐷 4 64 𝐺 112 (disk cross-section) 𝑀 3 = 𝐸 𝐷 4 12 𝐺 112 (square cross-section) (27) Lionel Gélébart
Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

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

Spatial resolution

Bending moment relative error square disk disk+composite voxels is very low even for the lowest resolution 6 × 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 (𝜎 11 𝑥 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 × 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). Heterogeneous beam bending Figure 5 (right) represents the evolution of the bending moments (or bending stiffness as 𝐺 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) × 10 7 for the

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates 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. 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,

Square

Disk Disk with Composite Voxel Table 3 Bending loading case: number of iterations at convergence for the various beams, contrasts, resolutions and use of composite voxels.

Resolution 𝐶 = 0 𝐶 = 1 𝐶 = 1000 𝐶 = 0 𝐶 = 1 𝐶 = 1000 𝐶 = 0 𝐶 = 1 𝐶 = 1000
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

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates 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 𝜎 1 1 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).

-500 -250 0 250 500

𝜎 11 [MPa]
-2,000 -1,000 0 1,000 2,000

𝜎 11 [MPa] Figure 6 Axial stress field 𝜎 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.

Beam torsion

The torsion loading applies to the unit-cell Ω * with 𝛼 = 𝐺 132 = -𝐺 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 (1𝑗) 𝑗=1,2,3 and average stress components (𝑖 𝑗) 𝑖=1,2𝑗=1,2 equal to zero. All the propositions provide the same torsion moment 𝑀 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

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates practice, 𝐺 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

𝑀 1 = 𝜇 𝜋𝐷 4 32 𝛼 (disk cross-section) 𝑀 1 = 𝜇𝐷 4 𝛼 × 0.1406 (square cross-section). ( 28
)
For the square cross-section, the expression is deduced from numerical simulation or from Fourier series evaluated numerically [START_REF] Bruhns | Advanced Mechanics of Solids[END_REF].

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

Spatial resolution

Torsion moment relative error square disk disk+composite voxels 5 × 10 -4 for the lowest resolution 40 and decreases up to 2 × 10 -4 for the highest resolution 160. The shear stress field 𝜎 13 , observed in Figure 8 (down), is not a simple linear function of 𝑥 2 but also depends on 𝑥 3 , emphasizing that the solution, obtained after 9 iterations, is quite different from the initial stress field (linear function of 𝑥 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 𝑥 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 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

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates composite voxels: decreasing from 12 × 10 -4 up to 4 × 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. 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 × 10 8 N m to 1.41 × 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.

Heterogeneous beam torsion

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 Square Disk Disk with Composite Voxel Table 4 Torsion loading case: number of iterations at convergence for the different beams, contrasts, resolutions and use of composite voxels.

Resolution 𝐶 = 0 𝐶 = 1 𝐶 = 1000 𝐶 = 0 𝐶 = 1 𝐶 = 1000 𝐶 = 0 𝐶 = 1 𝐶 =
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

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates 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.

Non-linear beam application

The aim of this final subsection is to demonstrate the ability of the method to account for nonlinear behaviors. For that purpose, a beam with a square cross-section (𝐷 × 𝐷 with 𝐷 = 100 mm) and a cube inclusion (𝑑 × 𝑑 × 𝑑 with 𝑑 = 75 mm) is used. Note that, in order to emphasize the effect of the inclusion on the beam behavior, the inclusion size 𝑑 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 𝑅 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 𝐶 = 0 for a void inclusion), 10 MPa (no elastic contrast with the matrix, 𝐶 = 1), 200 MPa (elastic contrast 𝐶 = 20 as observed

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates in glass/epoxy composites for example).

For bending tests, the applied load corresponds to a null average strain 𝐸 and the average strain gradient component 𝐺 * 112 = 𝐺 112 is set to 0.001 mm with 𝐺 132 = -𝐺 123 = 0.001 mm -1 at the end of the loading. Finally, in both tests, the load is discretized in 20 steps.

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 𝑥 2 ∈ [-𝑦 0 , 𝑦 0 ]) where the material behaves elastically. The second comes from the rest of the beam where the material flows plastically with 𝜎 11 = 𝑅 0 . The coordinate 𝑦 0 that separates the two regions corresponds to the coordinate for which the stress in the elastic region reaches 𝑅 0 = 𝐸𝜀 0 = 𝐸𝐺 112 𝑦 0 . If the loading is too low, the whole beam remains elastic and 𝑦 0 = 𝐷/2, so that the expression can be summarized by 𝑦 0 = min(𝐷/2, 𝑅 0 /(𝐸𝐺 112 )). Then, the evaluation of the bending moment from the two contributions reads

𝑀 = 𝑀 1 + 𝑀 2 = 𝐸 𝐷 (2𝑦 0 ) 3 12 𝐺 112 + 𝐷𝑅 0 𝐷 2 4 -𝑦 2 0 with 𝑦 0 = min 𝐷 2 , 𝑅 0 𝐸𝐺 112 . (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

𝑀 = 𝑀 1 + 𝑀 2 = 𝜇𝐺 123 𝜋𝑟 4 0 + 2𝜋𝑅 0 √ 3 𝐷 3 24 - 𝑟 3 0 3 with 𝑟 0 = min 𝐷 2 , 𝑅 0 / √ 3 2𝜇𝐺 112 . ( 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. 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.

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

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 𝐶 = 0 (void inclusion), 𝐶 = 1 (no elastic contrast) and 𝐶 = 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 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 𝐶 = 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 𝐶 = 1, the stiff inclusion with 𝐶 = 20 stiffens the elastic behavior, increases the yield stress but does not change the plastic limit, • compared to the heterogeneous beam with 𝐶 = 1, the porous inclusion with 𝐶 = 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 𝐶 = 1, is quite low for a void inclusion (𝐶 = 0), the increase is significant for the stiff inclusion (𝐶 = 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 (𝐶 = 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 𝐶 = 1, the computation time to obtain converged results with resolution 40 is 0.5 minute.

Homogeneous 𝐶 = 0 𝐶 = 1 𝐶 = 20

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

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates 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 [START_REF] Gélébart | AMITEX_FFTP[END_REF]) 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF], 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 [START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF] 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 [START_REF] Cartraud | Computational homogenization of periodic beam-like structures[END_REF]. 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 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 ũ on opposites faces with normal parallel to 𝒆 1 (i.e. ũ periodic on 𝜕 1 Ω) with the property ũ (𝑥 1 + ℎ 1 , 𝑥 2 , 𝑥 3 ) = ũ (𝑥 1 , 𝑥 2 , 𝑥 3 ), ∀𝑥 2 , 𝑥 3 . Using x, a linear average of 𝑥 along 𝒆 1 , leads to ũ1,1 (𝑥 2 , 𝑥 3 ) = ũ2,1 (𝑥 2 , 𝑥 3 ) = ũ3,1 (𝑥 2 , 𝑥 3 ) = 0, ∀𝑥 2 , 𝑥 3 . The consequences on the average strain fluctuations are that the conditions ⟨ ũ1,1 ⟩ =

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

• only the components 11, 22 and 12 of 𝐸 are prescribed for the average strain,

• only the components 11, 22 and 12 of 𝜒 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. D Characterization of composite voxels: beam with cylinder crosssection

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 𝑅 cv , radius of the encompassing circle, centered on the voxel. The expressions are as follows:

𝑆 = 𝑆 1 + 𝑆 2 𝑆 1 = 𝑅 2 0 arccos 𝑑 i 𝑅 0 -𝑑 i 𝑅 2 0 -𝑑 2 i 𝑆 2 = 𝑅 2 cv arccos 𝑑 1 -𝑑 i 𝑅 cv -(𝑑 1 -𝑑 i ) 𝑅 2 cv -(𝑑 1 -𝑑 i ) 2 𝑑 i = 𝑅 2 0 -𝑅 2 cv + 𝑑 2 1 2𝑑 1 𝑓 = 𝑆 𝜋𝑅 2 cv . (D.1)
𝑑 1 being the distance between the center of the beam section and the center of the voxel, and 𝑅 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𝑅 cv which corresponds to a threshold of 1 % of volume fraction (i.e. below this value, a voxel is not considered as composite).

Figure 1

 1 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.

Figure 3

 3 Figure3Relative 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.

Figure 4

 4 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.

Figure 5

 5 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.

Figure 7

 7 Figure7Relative 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.

Figure 8

 8 Figure8Homogeneous beam under torsion loading: shear stress field 𝜎 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.

Figure 9

 9 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.

Figure 10

 10 Figure 10 Shear stress 𝜎 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.

  Figure 11 (a) Bending moment for a square cross-section. (b) Torsion moment for a disk cross-section.

Figure 12

 12 Figure12Moments for a square cross-section and two resolutions 𝑟 . Black curves correspond to homogeneous beams (no inclusion) obtained analytically for bending and numerically for torsion. Legend applies to both figures.

  𝑬 11 𝑬 12 𝐸 13 𝑬 12 𝑬 22 𝐸 23 𝐸 13 𝐸 23 𝐸 33 and 𝜒 = 𝝌 11 𝝌 12 𝜒 13 𝝌 12 𝝌 22 𝜒 23 𝜒 13 𝜒 23 𝜒 33 = 𝑮 113 𝑮 123 𝐺 133 𝑮 123 𝑮 223 𝐺 233 𝐺 133 𝐺 133 𝐺 333 .

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates components of both the average strain 𝐸 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: 𝐸 =

  , see Section 2. The non-vanishing components (𝐸 22 , 𝜒 22 ) of 𝐸 𝑃 and 𝜒 𝑃 define the two applied loadings: (1, 0) for the tensile loading and (0, 1) for the bending loading. They directly allow to evaluate 𝐴 22 = 𝑁 22 /𝐸 22 and 𝐷 22 = 𝑀 22 /𝜒 22 . The corresponding loadings applied to the enlarged unit-cell Ω

* * are:

Table 1

 1 Axial stiffness 𝐴 22 for laminate case: comparison between analytical results and two FFT-based methods.

	Resolution Analytical	Present work (Nguyen et al. 2008)
	8 × 8 4 × 4	3.07692307 3.07692303 3.07692307 3.07692303	3.077

  Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates

	Lionel Gélébart
	Resolution Analytical [GPa/m]	Present study	Rel. error	(Nguyen et al. 2008)	Rel. error
	16 × 16 32 × 32 64 × 64 128 × 128 256 × 256 512 × 512 1024 × 1024 0.0038003663 0.0038003663 0.0038003663 0.0038003663 0.0038003663 0.0038003663 0.0038003663	0.0037903502 0.0037978623 0.0037997403 0.0038002098 0.0038003272 0.0038003566 0.0038003638	26 × 10 -4 7 × 10 -4 2 × 10 -4 4 × 10 -5 1 × 10 -5 2 × 10 -6 7 × 10 -7	0.00382 0.003805 0.003802	52 × 10 -4 13 × 10 -4 5 × 10 -4

  𝑬 11 𝐸 12 𝐸 13 𝐸 12 𝐸 22 𝐸 23 𝐸 13 𝐸 23 𝐸 33 , 𝜒 3 = 𝑮 113 -𝛼 𝐺 133 -𝛼 𝐺 223 𝐺 233 𝐺 133 𝐺 233 𝐺 333 The component 𝐸 11 prescribes a tensile loading, the components 𝐺 112 and 𝐺 113 prescribe two bending loadings respectively, and 𝛼 = 𝐺 132 = -𝐺 123 prescribes the torsion loading. According to the arguments proposed in Appendix A.1 based on the vanishing components of 𝐸 and 𝐺 imposed by the periodicity condition in direction 1, only the component 11 of the 𝐸, 𝜒 2 and 𝜒 3 tensors can be prescribed. However, the components 𝐺 132 = -𝐺 123 = 𝛼 were successfully prescribed for the torsion loading in the simulations below. Arguments regarding the symmetry of the problem are given in Appendix A.1.

						The
	macroscopic applied strain is limited to components in boldface in
			𝑮 112 𝐺 122	𝛼
	𝐸 =	and 𝜒 2 =	𝐺 122 𝐺 222 𝐺 232	.	(21)
			𝛼	𝐺 232 𝐺 332

  -1 at the end of the loading. If the gradients 𝐺 222 and 𝐺 332 , cannot be prescribed (see Appendix A), the choice of 𝐺 * 222 and 𝐺 * 332 influences the convergence of the iterative algorithm. The choice 𝐺 * 222 = 𝐺 * 332 = -𝑣𝐺 112 made in the previous section with linear behaviors is kept. Other strain gradient components 𝐺 * 𝑖 𝑗𝑘 are set to zero. For torsion tests, the non-zero components of the strain gradient are 𝐺 * 132 = 𝐺 132 and 𝐺 * 123 = 𝐺 123

Table 5

 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.

  Table B.1 Bending stiffness 𝐷 22 for the laminate case: comparison between analytical results and the proposed FFT-based method.All the simulation cases reported in[START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF] have been reproduced with our approach and results are gathered in Table C.2 for various contrasts (contrast = 𝐸 inclusion /𝐸 matrix in [GPa]) and loading cases. The absolute values are respectively the tension and bending stiffness 𝐴 22 = 𝑁 22 /𝐸 22 and 𝐷 22 = 𝑀 22 /𝜒 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[START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF], with a relative difference of strain energy between two iterations less than 10 -4 . Table C.2 Comparison analysis. łIter. ž means łIterationsž. łRel. err. ž means łRelative errorž. Starred columns refer to (Nguyen et al. 2008). Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates Table C.3 Glass Epoxy composite.

	Lionel Gélébart
							Resolution	Present study	Iterations
						Tension Bending	16 × 16 32 × 32 64 × 64 128 × 128 16 × 16 32 × 32 64 × 64 128 × 128	0.621 221 × 10 9 0.622 367 × 10 9 0.622 609 × 10 9 0.622 700 × 10 9 3.889 901 074 × 10 5 3.906 817 410 × 10 5 3.910 938 081 × 10 5 3.912 030 558 × 10 5	68/26 80/35 113/28 122/32 62/20 65/20 71/20 71/20
					Resolution Analytical [GPa/m]	Present study	Rel. error
					16 × 16 32 × 32 64 × 64 128 × 128	19.949927685076197 19.897348729113546 -0.26% 19.949927685076197 19.936783319118270 -0.07% 19.949927685076197 19.946641214046725 -0.016% 19.949927685076197 19.949106561516157 -0.004%
			Resolution FEM*	Present study Rel. err.		FFT* Rel. err.	Iter.	Iter. *
	Contrast = 0	Tension Bending	16 × 16 32 × 32 64 × 64 128 × 128 16 × 16 32 × 32 64 × 64 128 × 128	0.5840 0.5840 0.5840 0.5840 8.224 × 10 -4 8.224 × 10 -4 8.224 × 10 -4 8.224 × 10 -4	0.5822 -0.3% 0.5834 -0.1% 0.5838 -0.07% 0.5838 -0.07% 8.204 × 10 -4 -0.24% 8.219 × 10 -4 -0.06% 8.223 × 10 -4 +0.015% 8.224 × 10 -4 -0.004%	0.5826 -0.24% 0.5837 -0.05% 0.5839 -0.01% 0.5840 ∼ 0% 0.22% 8.242 × 10 -4 0.05% 8.228 × 10 -4 0.01% 8.225 × 10 -4 8.224 × 10 -4 ∼ 0 %	28/8 45/11 44/14 48/15 19/8 19/8 19/8 21/8	8 9 9 9 8 8 8 8
	Contrast = 0.001	Tension Bending	16 × 16 32 × 32 64 × 64 128 × 128 16 × 16 32 × 32 64 × 64 128 × 128	58.495 58.495 58.495 58.495 0.08225 0.08225 0.08225 0.08225	58.329 -0.28% 58.455 -0.07% 58.484 -0.02% 58.487 -0.01% 0.08206 -0.23% 0.08221 -0.05% 0.08224 -0.005% 0.08225 -0.006%	58.363 -0.23% 58.472 -0.04% 58.502 +0.01% 58.506 +0.02% 0.08244 +0.23% 0.08230 +0.06% 0.08227 +0.02% 0.08226 +0.01%	28/8 39/11 45/14 48/15 50/8 35/8 20/8 20/8	8 9 9 9 8 8 8 8
	Contrast = 10	Tension Bending	16 × 16 32 × 32 64 × 64 128 × 128 16 × 16 32 × 32 64 × 64 128 × 128	0.16288 0.16288 0.16288 0.16288 1.0652 × 10 -4 1.0652 × 10 -4 1.0652 × 10 -4 1.0652 × 10 -4	0.16267 -0.13% 0.16284 -0.02% 0.16287 -0.01% 0.16290 +0.01% 1.0597 × 10 -4 -0.5% 1.0638 × 10 -4 -0.13% 1.0649 × 10 -4 -0.03% 1.0652 × 10 -4 -0.002% 1.0654 × 10 -4 0.16209 -0.49% 0.16257 -0.19% 0.16279 -0.06% 0.16287 -0.01% 1.0701 × 10 -4 +0.46% 1.0659 × 10 -4 +0.07% 1.0654 × 10 -4 +0.02% +0.02%	49/14 39/17 60/17 76/17 40/17 42/17 51/17 45/17	15 16 16 16 18 18 18 18

C Plate with periodic fibers: comparison with

[START_REF] Nguyen | Green's operator for a periodic medium with traction-free boundary conditions and computation of the effective properties of thin plates[END_REF]

Journal of Theoretical, Computational and Applied Mechanics July 2022 jtcam.episciences.org

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

Ethics approval and consent to participate Not applicable.

Consent for publication Not applicable.

Lionel Gélébart

Extension of FFT-based methods to strain gradient loadings ś Homogenization of beams and plates ⟨ ũ2,1 ⟩ = ⟨ ũ3,1 ⟩ = 0 without any constraint on the other components imply that ⟨ ε11 ⟩ = 0 without any constraint on other components. Above, since ε is the symmetric part of the gradient, shear components do not systematically vanish. It can be observed additionally that ũ𝑖,𝑗 (𝑥 1 + ℎ 1 , 𝑥 2 , 𝑥 3 ) = ũ𝑖,𝑗 (𝑥 1 , 𝑥 2 , 𝑥 3 ), ∀𝑥 2 , 𝑥 3 and 𝑖 = 1, 2, 3, 𝑗 = 2, 3. The vanishing components 𝑖 𝑗𝑘 of ⟨ ũ𝑖,𝑗𝑘 ⟩ arising from this periodicity are reported below (the second column is simply an inversion of the derivation order: Above, since G is symmetric with respect to the first two indices, the vanishing components of 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 ⟨ ε𝑖 𝑗 ⟩ = 0, ∀𝑖, 𝑗. The consequences on the average strain gradient fluctuations are:

with no constraint on other components. Hence, when using 𝜀 = 𝐸 + 𝜀 * + ε in Problem ( 8): • all components of 𝐸 are prescribed for the average strain, • only the nine components of the gradient of 𝜀 * corresponding to G𝑖 𝑗𝑘 = 0, are prescribed for the average strain gradient (components 𝐺 𝑖 𝑗 𝑗 cannot be prescribed).

The other components of the strain gradient of 𝜀 * (i.e. 𝐺 * 𝑖 𝑗 𝑗 ) can be set to arbitrary values.

A.3 Application to the periodic beam problem

The problem involves a single periodicity condition along 𝒆 1 . The conclusions of Appendix A.1 apply. When using 𝜀 = 𝐸 + 𝜒 3 𝑥 3 + 𝜒 2 𝑥 2 + ε in Problem (20): • only the component 𝐸 11 is prescribed for the average strain.

• only the component 11 of 𝜒 2 and 𝜒 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 𝜒 3 and 23 of 𝜒 2 cannot be prescribed. However, in the case of symmetric unit-cells, the axial component ũ1 of the fluctuation displacement, is equal to 0 on the two opposite faces. In addition to ⟨ ũ2,13 ⟩ = ⟨ ũ3,12 ⟩ = 0 reported in Appendix A.1, we have ⟨ ũ1,23 ⟩ = ⟨ ũ1,32 ⟩ = 0, so that the components of the strain gradient fluctuation G213 and G312 vanish and the corresponding components of the average strain gradient can be prescribed. In the general case (no specific symmetry), we have 𝐺 213 = -𝛼 + 1 2 ⟨ ũ1,23 ⟩ and 𝐺 312 = 𝛼 + 1 2 ⟨ ũ1,32 ⟩. If the quantity ⟨ ũ1,23 ⟩ = ⟨ ũ1,32 ⟩ 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 ⟨ ũ𝑖,𝑗 𝑗 ⟩ that is able to vanish the applied strain gradient 𝐺 𝑖 𝑗 𝑗 . Actually if ⟨ ũ1,23 ⟩ increases to reduce the magnitude of 𝐺 213 , on the other hand, it will increase 𝐺 312 .

A.4 Application to the periodic plate

The problem involves two periodicity conditions along 𝒆 1 and 𝒆 2 . In 𝜀 = 𝐸 + 𝜒𝑥 3 + (∇ 𝒖) sym in Problem (11):

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.