figurec

# Adaptation and validation of FFT methods for homogenization of lattice based
materials

S. Lucarini1, L. Cobian1,2, A. Voitus1 and J. Segurado1,2
1 Fundación IMDEA Materiales, C/ Eric Kandel 2, 28906, Getafe, Madrid, Spain
2 Universidad Politécnica de Madrid, Department of Materials Science,
E.T.S.I. Caminos, C/ Profesor Aranguren 3, 28040, Madrid, Spain

###### Abstract

An FFT framework which preserves a good numerical performance in the case of
domains with large regions of empty space is proposed and analyzed for its
application to lattice based materials. Two spectral solvers specially suited
to resolve problems containing phases with zero stiffness are considered (1) a
Galerkin approach combined with the MINRES linear solver and a discrete
differentiation rule and (2) a modification of a displacement FFT solver which
penalizes the indetermination of strains in the empty regions, leading to a
fully determined equation. The solvers are combined with several approaches to
smooth out the lattice surface, based on modifying the actual stiffness of the
voxels not fully embedded in the lattice or empty space. The accuracy of the
resulting approaches is assessed for an octet-lattice by comparison with FEM
solutions for different relative densities and discretization levels. It is
shown that the adapted Galerkin approach combined with a Voigt surface
smoothening was the best FFT framework considering accuracy, numerical
efficiency and $h$-convergence. With respect to numerical efficiency it was
observed that FFT becomes competitive compared to FEM for cells with relative
densities above $\approx$7%. Finally, to show the real potential of the
approaches presented, the FFT frameworks are used to simulate the behavior of
a printed lattice by using direct 3D tomographic data as input. The approaches
proposed include explicitly in the simulation the actual surface roughness and
internal porosity resulting from the fabrication process. The simulations
allowed to quantify the reduction of the lattice stiffness as well as to
resolve the stress localization of $\approx$ 50% near large pores.

##  1 Introduction

The latest improvements in additive manufacturing techniques (AM) have made
possible the fabrication of micro- and nano-architected metamaterials for
mechanical applications with tailored stiffness, strength, toughness and
energy absorption capacity [1]. In most cases, the topology of these
architected materials at the lower length scale is based on the periodic
repetition of a unit cell made up of bars or shells forming a lattice. These
lattice materials can reach high strength-to-weight ratios and other specific
properties as well as achieving non-standard elastic responses such as
negative compressibility or zero Poisson’s ratio [1, 2, 3]. In addition, unit
cells of lattice metamaterials can be designed using unstable structural
elements to achieve energy dissipation with fully reversible deformation or
programmable/switchable properties [4].

In order to design an optimal lattice for a specific property is fundamental
to perform accurate and computationally efficient simulations of the effective
response and the local fields developed within the microstructure. Although
discrete mechanical models are a good approach for capturing the overall
response [5], full-field homogenization approaches allow obtaining a more
complete, and often more accurate, result. Under this approach, a
micromechanical problem is solved on a representative volume element (RVE)
which contains a unit cell or a collection of unit cells in which the lattice
geometry is explicitly represented. In particular, the Finite Element Method
(FEM) is the most common approach to carry out these simulations since the
designed _ideal_ microstructure can be accurately reproduced with adaptative
meshes [6]. However, the real microstructures arising from an AM process can
differ from the _ideal_ ones designed in-silico, and usually present
geometrical deviations from the target geometry, porosity and surface
roughness [7]. Although these features can be neglected for the lattice design
phase, a reliable simulation of the material deformation in the non linear
regime which considers its irreversible response or fracture should take these
differences into account. FEM models of realistic geometries including surface
roughness or porosity are very difficult to generate and mesh. Moreover, the
real geometrical data is usually obtained from tomographic images [8, 7], and
transforming this data into a FEM model is a very complex process that results
in a very large number of elements.

In this context, the use of spectral methods, based on the Fast Fourier
Transform (FFT) algorithm, can be an interesting alternative. FFT methods do
not require meshing because the calculations are carried out on a regular grid
and the phase belonging to each point of the grid can be obtained directly
from digital images or tomographic data. In addition, FFT approaches are very
efficient, with a computational cost which grows as $n\log(n)$, improving the
FEM computational efficiency in the homogenization of bulk heterogeneous
materials by orders of magnitude [9, 10]. There are also other tangential
benefits coming from the use of spectral approaches for studying lattice
materials like the possibility of studying brittle fracture or ductile damage
using phase-field fracture or gradient damage approaches, methods that show a
very efficient performance in FFT [11, 12, 13].

The basic ideas of FFT homogenization were proposed by Moulinec and Suquet
[14]. Since then, different approaches have been proposed to improve the
convergence rate of the method. Some approaches are derived from the original
method, based on the use of Green’s functions for a reference medium and the
solution of the Lippmann-Schwinger equation [14, 15, 16, 17, 18, 19, 20, 21].
Alternative approaches have also been proposed, based on either using a
Galerkin approximation of the equilibrium using trigonometric polynomials to
discretize test and trial functions, [22, 23, 24, 25] or on solving the strong
form of equilibrium using displacements as unknown [26]. Despite the clear
potential benefits of FFT solvers for lattice materials, two main limitations
arise that have prevented their extensive use in this field. First, FFT
approaches present a convergence rate and accuracy strongly dependent on the
contrast between the phases represented in the domain. In the present case,
the contrast is infinite because a large amount of the RVE voxels are empty
space, with zero stiffness. Second, although the voxelized representation is
very convenient for using digital microstructure data, its use for
representing smooth geometries might result in poor results near the
boundaries. In addition, since FFT and inverse FFT transformations should
include all the points of the RVE including the empty space, the efficiency of
the method with respect to FEM will depend on the relative volume fraction,
and it will not be competitive for very small relative densities.

Regarding the first limitation, several studies can be found which aim at
overcoming the problems with the high phase stiffness contrast in FFT solvers
and try to extend their use to study materials with voids. Michel et al. [27]
developed an Augmented Lagrangian formulation to solve a non linear problem
including non-compatible fields. Although this formulation allows introducing
zero stiffness phases, it might require a very large number of iterations to
fulfil both stress equilibrium and strain compatibility, as noticed in [28].
Brisard and Dormieux [29] developed a variational formulation based on the
energy principle of Hashin and Shtrikman applied to a porous media. Their
approach allows to accurately predict the overall response of porous materials
but it involves the pre-computation of a consistent Green operator which is
computationally very expensive. More recently, a method for solving the
conductivity problem in the presence of voids has been developed by To and
Bonnet [30]. This approach is focused on solving the equilibrium only in the
bulk phases including a flux term at the inter-phase between bulk and void
phases. This method is suitable for scalar fields, but cannot be directly
extended to the vector and tensor fields that arise in the mechanical problem
since the flux term in the internal boundaries does not restrict tensor
components parallel to the interphase. Another recent method proposed by
Schneider [31] consists in searching solutions in a subspace of solutions on
which the homogenization problem is nondegenerate for the resolution of a
material with pores. In parallel to these techniques specifically suited for
porous materials, an efficient and simple alternative to improve the
convergence rate under very large phase contrast is the use of methods to
reduce the numerical oscillations that may occur due to Gibb’s phenomenon or
aliasing effects. A first possibility is filtering the high frequencies, as
proposed for example in [32, 33, 34]. A second possibility is replacing the
continuum differential operators in the formulation of the partial
differential equations in the real space by some finite difference
differentiation rule. This idea was first introduced in [35] who incorporated
the finite difference definitions of the derivatives in the real space using
the FFT algorithm and the definition of Fourier derivation using modified
frequencies. The finite difference stencil used in the real space defines the
particular form of the modified frequencies to be used. Among the different
discrete differential approaches, the so-called rotated scheme [36] shows a
significant reduction of noise and improves the convergence [10, 37]. Other
two finite differentiation schemes which are not based on the use of modified
frequencies were also proposed in [38] and [39] showing a clear improvement of
the spurious oscillations. A third approach to reduce the noise in the
solution was proposed by Eloh et. al, [40] who instead of using the DFT as the
discrete counterpart of the continuous Fourier transform, considered the
continuous Fourier transform of a piecewise constant operator in the real
space to derive _consistent periodized discrete Green operators._

The second disadvantage of using FFT homogenization for lattice materials
arises when, instead of using direct images or tomographic data to construct
the model, the objective is representing smoothly the lattice material
boundaries for general geometries defined analytically or through CAD models.
In these cases, the voxelized representation can be inaccurate for small
number of voxels, not achieving the actual relative density of the lattice.
This can lead to local inaccurate results due to the combination of a high
phase contrast and a non-smooth interface. This issue has been treated by
smoothening the sharp inter-phase using _composite voxels_ [33, 41]. These
approaches are based on defining the material response for the voxels crossed
by an internal interphase as the homogenization of the two phases contained.

Due to these limitations, only a few previous attempts of using FFT to study
lattice materials can be found on the literature. In all the cases the
motivation was to describe the actual lattice geometry, including
imperfections, using data obtained from tomographic images. In [42], the
augmented Lagrangian approach [27] was used only as a preprocessing step,
using tomographic images of a single strut as data, in order to determine an
equivalent diameter of the struts. Then, the simulation of full lattices was
done in FEM using ideal geometries with equivalent radius. In [8] and [43],
FFT was used to determine directly the response of lattice materials. In both
cases, the studies were limited to the linear elastic regime and the FFT
approaches used did not consider an infinite compliant phase for the empty
space, but used a material with small but finite stiffness to achieve
convergence instead. None of these works included an assessment of the
accuracy of the FFT approach in order to determine the effect of the
artificial stiffness used to represent the empty areas.

The present work presents a systematic and critical assessment of the accuracy
and efficiency of FFT approaches for predicting the mechanical response of
lattice based materials, both in the linear and the nonlinear regime, in order
to establish an optimal framework for the homogenization of this type of
materials. After a preliminary study, two linear FFT approaches are selected
as potential candidates, including a novel algorithm for RVEs including pores.
In parallel, different geometrical approaches to represent smoothly the
lattice geometries are combined with each approach. The accuracy and
efficiency of the different combinations of FFT solvers and geometrical
approaches are compared against FEM to model both the linear and non linear
responses of an ideal octet lattice. Finally, both frameworks are applied to
study a real octet lattice cell, including fabrication defects, obtained with
3D-tomography.

##  2 FFT homogenization for RVEs with empty areas

As previously discussed, a clear limitation for the use of FFT techniques for
homogenization is the deficient convergence rate in the case of
microstructures with a large contrast in the phase properties. This problem
becomes critical in the case of materials with voids where one of the phases
is infinitely compliant. In this case, the problem of the low convergence rate
(or no convergence in many cases) is superposed to the singularity of the
problem: the solution is not unique since any compatible strain field in the
empty phase is admissible. The problem has been circumvented on many occasions
by setting a very compliant elastic behavior for the void phase, but this
artificial stiffness might have an effect on the cell response, especially for
large volume fractions of empty voxels.

The proper adaptation of FFT homogenization algorithms to account for actual
zero stiffness has been studied almost from the first developments of this
numerical technique. The augmented Lagrangian approach [27] is a modification
of the original basic scheme to improve convergence for RVEs containing phases
with very low or zero stiffness. The algorithm is based on the combination of
two strain fields, one of them forced to be compatible, and two stress fields,
one of them forced to be in equilibrium. Then, the solution is obtained by the
iterative minimization of a Lagrangian. Although the method is in theory able
to resolve cases with infinite phase contrast, the convergence rate becomes
really small when controlling a residual which enforces both stress
equilibrium and strain compatibility in addition to the macroscopic
constraints [28]. Other acceleration methods also based on using additional
fields have been developed to account for infinite phase contrast as [44].
Nevertheless, these methods present a similar convergence rate when small
tolerances are imposed for both equilibrium and compatibility [28]. Another
potential method specifically developed for considering voids is the
variational approach by Brisard et.al. [29]. Although the idea presented is
very smart, the consistent differential operator derived is very complex to
compute and the authors themselves do not use it finally, but an approximate
version based on discrete derivatives.

As mentioned in the introduction, an alternative (or additional) way of
improving the convergence rate under large phase property contrast is using
methods that reduce the numerical oscillations. In particular, the standard
continuous differentiation scheme can be replaced by finite difference schemes
through different approaches as using staggered grids [38] or modifying the
Fourier derivative definition with modified frequencies which corresponds to
different finite difference stencils [35, 36]. These alternative
differentiation schemes can be combined with different FFT solvers, improving
the convergence of the original ones.

In this paper the Galerkin FFT approach [22, 23] combined with the use of
mixed control [25] and the rotated finite difference scheme [36] has been
chosen as a first candidate for the lattice material homogenization. This
combination has been made after a preliminary study and is based on its very
fast convergence rate but also on its relative simplicity and the ability of
the method to be efficiently extended for non linear cases. This approach does
not break the underdetermination of the solution, but allows to converge in
the presence of regions with zero stiffness to an equilibrated stress and
compatible strain with relatively low noise. Note also that any other Krylov
based approach with the appropriate reference medium and combined with the
same discretization scheme would eventually lead to similar results.

As second candidate, a new method based on a modification of the DBFFT
approach [26] is proposed. In this method, the standard equilibrium is
augmented with additional conditions for the void regions and interfaces to
break the underdetermination leading to a non-singular discrete system of
equations. In the next subsections, both methods will be presented including
their extension for non linear response.

###  2.1 Galerkin-FFT with discrete differences and mixed loading control

The FFT-Galerkin method was initially developed by Vondrějc et al. [22] to
homogenize the elastic behavior of fully dense heterogeneous materials. This
approach presents a very fast convergence rate for limited phase contrast, but
is not able by itself to converge in the presence of infinitely compliant
phases. In order to extend this scheme for RVEs with empty regions, the
original method is adapted by changing the iterative linear solver,
introducing an alternative differentiation scheme [36] and using mixed
macroscopic control [26]. Both the original method and these modifications
will be presented below.

Following [22], starting from the weak form of the equilibrium in small
strains for a given heterogeneous periodic domain $\Omega$ the next equation
can be derived

$\mathcal{F}^{-1}\left\\{\widehat{\mathbb{G}}\left(\boldsymbol{\xi}\right):\mathcal{F}\left\\{\boldsymbol{\sigma}\left(\mathbf{x}\right)\right\\}\right\\}=\mathbf{0}\text{
,}$ (1)

where $\mathbf{x}$ represents the spatial position, $\boldsymbol{\xi}$ the
spatial frequency vector, $\boldsymbol{\sigma}$ is the Cauchy stress
—determined by the local constitutive equations— and $\widehat{\mathbb{G}}$ is
the Fourier transform of a linear operator which projects any arbitrary tensor
field into its compatible part. The Fourier transform and the inverse Fourier
transform are represented in eq. (1) by $\mathcal{F}$ and $\mathcal{F}^{-1}$
respectively.

In equation (1) the domain is discretized in a regular grid in which each
voxel center $\mathbf{x}=(x_{1},x_{2},x_{3})$ is given by

$x_{p}=\frac{L_{p}}{2N_{p}}+\frac{L_{p}}{N_{p}}n_{p},\text{with}\
n_{p}=0,\dots,N_{p}-1;\ p=1,2,3\text{ ,}$

where $L_{p}$ and $N_{p}$ stand for the length of the cell edge and number of
voxels in direction $p$. The discrete form of the equation (1) is a linear
system of algebraic equations in which the unknown is the value of the strain
$\boldsymbol{\epsilon}$ at the center of each voxel. The frequency vector
$\boldsymbol{\xi}$ is given by

$\xi_{p}=iq_{p}N_{p}/L_{p}\quad\text{with}\quad
q_{p}=2\pi\frac{n_{p}-N_{p}/2}{N_{p}}\in[-\pi,\pi]\text{ ,}$ (2)

where $i=\sqrt{-1}$ is the imaginary unit. The Fourier transforms correspond
to the direct and inverse Discrete Fourier Transform that are carried out
using the FFT algorithm. The ability of the FFT algorithm to reduce the
computational cost of the DFT transforms from $\mathcal{O}(n^{2})$ to
$\mathcal{O}(n\log n)$ is the main reason behind the high performance of
spectral solvers.

The macroscopic state is provided as a combination of macroscopic strain
components $kl$,
$\overline{\boldsymbol{\varepsilon}}=\overline{\varepsilon}_{kl}(\mathbf{e}_{k}\otimes\mathbf{e}_{l})^{sym}$
or macroscopic stress components $KL$
$\overline{\boldsymbol{\sigma}}=\overline{\sigma}_{KL}(\mathbf{e}_{K}\otimes\mathbf{e}_{L})^{sym}$
with $kl\cap KL=\emptyset$. For these general loading conditions [25] the
expression of the projector operator in the Fourier space follows

$\widehat{G}_{ijkl}(\boldsymbol{\xi})=\left\\{\begin{array}[]{lr}I^{s}_{ijKL}&\text{if
}\boldsymbol{\xi}=\mathbf{0}\text{ for components $KL$ }\\\ 0_{ijkl}&\text{if
}\boldsymbol{\xi}=\mathbf{0}\text{ for components $kl$}\\\ 0_{ijkl}&\text{for
Nyquist frequencies}\\\
\left[I^{s}_{ijkl}\xi_{j}\xi_{l}\right]^{-1}\xi_{j}\xi_{l}&\text{for
}\boldsymbol{\xi}\neq\mathbf{0}\end{array}\right.\text{ ,}$ (3)

where $\mathbb{I}^{s}$ is the fourth order symmetric identity tensor and
$\widehat{\mathbb{G}}$ accounts for major and minor symmetries.

In the case of linear elastic phases eq. (1) yields

$\mathcal{F}^{-1}\left\\{\widehat{\mathbb{G}}\left(\boldsymbol{\xi}\right):\mathcal{F}\left\\{\mathbb{C}\left(\mathbf{x}\right):\boldsymbol{\varepsilon}\left(\mathbf{x}\right)\right\\}\right\\}=-\mathcal{F}^{-1}\left\\{\widehat{\mathbb{G}}\left(\boldsymbol{\xi}\right):\mathcal{F}\left\\{\mathbb{C}\left(\mathbf{x}\right):\overline{\boldsymbol{\varepsilon}}\left(\mathbf{x}\right)-\overline{\boldsymbol{\sigma}}\right\\}\right\\}\text{
,}$ (4)

where $\mathbb{C}$ is the local fourth order stiffness tensor. Equation (4)
represents a linear system of equations in which the left-hand side is a
symmetric semidefinite positive linear operator acting on a discrete strain
field

$\mathcal{A}(\bullet)=\mathcal{F}^{-1}\left\\{\widehat{\mathbb{G}}\left(\boldsymbol{\xi}\right):\mathcal{F}\left\\{\mathbb{C}\left(\mathbf{x}\right):\bullet\right\\}\right\\}$
(5)

and the right-hand side is the independent term $\mathbf{b}$

$\mathbf{b}=-\mathcal{F}^{-1}\left\\{\widehat{\mathbb{G}}\left(\boldsymbol{\xi}\right):\mathcal{F}\left\\{\mathbb{C}\left(\mathbf{x}\right):\overline{\boldsymbol{\varepsilon}}\left(\mathbf{x}\right)-\overline{\boldsymbol{\sigma}}\right\\}\right\\}.$

The system $\mathcal{A}(\boldsymbol{\varepsilon})=\mathbf{b}$ can be solved
efficiently using the Conjugate Gradient method (CG) method for domains with a
relatively low stiffness contrast. The residual of this system of equations is
defined as the $L_{2}$ norm of the difference between the linear operator
applied on the candidate solution and the right-hand side over the norm of the
right-hand side,

$r_{lin}=\frac{\left\|\mathcal{A}(\boldsymbol{\varepsilon})-\mathbf{b}\right\|_{L_{2}}}{\left\|\mathbf{b}\right\|_{L_{2}}}$
(6)

reaching the solution when $r_{lin}$ is below a tolerance. Note that eq. (4)
is undetermined independently of the phase properties because any incompatible
field $\boldsymbol{\varepsilon}^{INC}$ added to the solution still fulfils eq.
(4),

$\mathcal{A}(\boldsymbol{\varepsilon}+\boldsymbol{\varepsilon}^{INC})=\mathcal{A}(\boldsymbol{\varepsilon})=\mathbf{b}.$

Nevertheless, the solution of eq. (4) is unique in the subspace of compatible
strain fields, so the CG is able to handle the underdetermination and recovers
the unique compatible solution.

In the case of lattices, there are empty regions in the RVE, which do not
transfer stresses and that should be properly accounted for with zero
stiffness. The solution of the problem in the full RVE is singular even in the
subspace of compatible strain fields, because any compatible strain field
which is zero outside of the empty region can be added to the solution
providing the same equilibrated stress. This singularity is transferred to the
numerical method increasing the underdetermination of the system also to
compatible fields. As a result, it is observed that the CG method is not able
to reach convergence and the Minimal Residual Method (MINRES), an alternative
Krylov subspace solver able to handle efficiently singular systems, is used as
linear solver to ovecome this limitation.

Finally, in order to improve both smoothness of the solution and convergence
rate, an alternative discretization scheme is used, the rotated forward finite
difference rule [36]. This discretization is introduced through modified
frequencies in the Fourier derivation leading to an alternative projection
operator (eq. 3). The modified frequencies correspond to

$\xi^{\prime}_{p}=i2\frac{N_{p}}{L_{p}}\tan{\left(q_{p}/2\right)}\prod_{p=1}^{d}\frac{1}{2}\left(1+e^{iq_{p}}\right)\quad\text{with}\quad
q_{p}=2\pi\frac{n_{p}-N_{p}/2}{N_{p}}\in[-\pi,\pi]\text{ ,}$ (7)

where $d=1,2$ or 3 is the space dimension. It is important to remark that
standard Galerkin approach does not converge using standard Fourier
discretization and a conjugate gradient solver. On the contrary, the
combination of the alternative differentiation scheme and the use of MINRES
allows the Galerkin approach to reach a solution in a relatively small number
of iterations, as it will be shown in the numerical results.

####  2.1.1 Non linear extension

In order to take into account material non linearities on the Galerkin FFT
formulation the macroscopic strain/stress history is applied as function of
the time (or pseudo-time for rate independent materials) in several time
increments. The non linear equilibrium at each increment is solved iteratively
using the Newton-Raphson method, as proposed in [23]. If the solution at time
$t$ and the macroscopic stress and strain components applied at $t+\Delta t$,
$\overline{\boldsymbol{\varepsilon}}^{t+\Delta
t},\overline{\boldsymbol{\sigma}}^{t+\Delta t}$ are known, the non linear
equation at time $t+\Delta t$ is linearized at each iteration $i$ around the
strain field at previous iteration $\boldsymbol{\varepsilon}^{i-1}$. Let
$\delta\boldsymbol{\varepsilon}\left(\mathbf{x}\right)$ be the strain field
correction to be obtained at iteration $i$, then the linearized stress
corresponds to

$\boldsymbol{\sigma}^{i}\left(\mathbf{x}\right)=\boldsymbol{\sigma}^{i-1}\left(\mathbf{x}\right)+\left.\frac{\partial\boldsymbol{\sigma}\left(\mathbf{x}\right)}{\partial\boldsymbol{\varepsilon}}\right|_{\epsilon=\epsilon^{i-1}}:\delta\boldsymbol{\varepsilon}\left(\mathbf{x}\right)=\boldsymbol{\sigma}^{i-1}\left(\mathbf{x}\right)+\mathbb{C}^{i-1}\left(\mathbf{x}\right):\delta\boldsymbol{\varepsilon}\left(\mathbf{x}\right)\text{
,}$ (8)

where $\mathbb{C}^{i-1}$ is the material consistent tangent evaluated using
the solution of the previous iteration $i-1$. The equilibrium equation
linearized at $\boldsymbol{\varepsilon}^{i-1}$ reads

$\mathcal{F}^{-1}\left\\{\widehat{\mathbb{G}}\left(\boldsymbol{\xi}\right):\mathcal{F}\left\\{\mathbb{C}^{i-1}\left(\mathbf{x}\right):\delta\boldsymbol{\varepsilon}\left(\mathbf{x}\right)\right\\}\right\\}=-\mathcal{F}^{-1}\left\\{\widehat{\mathbb{G}}\left(\boldsymbol{\xi}\right):\mathcal{F}\left\\{\boldsymbol{\sigma}\left(\boldsymbol{\varepsilon}^{i-1}\left(\mathbf{x}\right)\right)-\overline{\boldsymbol{\sigma}}^{t+\Delta
t}\right\\}\right\\}\text{ ,}$ (9)

where the solution at the previous time step enters in the first iteration as
$\boldsymbol{\varepsilon}^{0}\left(\mathbf{x}\right)=\boldsymbol{\varepsilon}^{t}\left(\mathbf{x}\right)+\overline{\boldsymbol{\varepsilon}}^{t+\Delta
t}-\overline{\boldsymbol{\varepsilon}}^{t}$.

The left-hand side of equation (9) corresponds to a linear operator
$\mathcal{A}_{i}$, acting on the correction strain field
$\delta\boldsymbol{\varepsilon}$, that is equivalent to the one defined in eq.
(5) using the tangent stiffness at iteration $i-1$ instead of the elastic one.
Similarly, the right-hand side forms an independent vector $\mathbf{b}_{i}$.
The solution of the non linear problem at each time increment is obtained
solving the linear equation (9) for each Newton iteration, and adding the
successive solution corrections until the convergence is reached.

Special care has to be paid to the definition of the residuals in the non
linear case. As in the original approach [23], two residuals are proposed for
the non linear solver, but they are redefined in order to avoid oversolving.
The linear residual controls the accuracy of the solution in the linear system
resulting from each Newton iteration. This linear system
$\mathcal{A}_{i}(\delta\boldsymbol{\varepsilon})=\mathbf{b}_{i}$ is solved
iteratively up to a given tolerance and the non linear algorithm becomes an
inexact damped Newton method [45]. The standard non-linear approach for
Galerkin FFT establish the convergence criterion on the relative residual,
defined as

$\frac{\left\|\mathcal{A}_{i}\left(\mathbf{x}\right)-\mathbf{b}_{i}\right\|_{L_{2}}}{\left\|\mathbf{b}_{i}\right\|_{L_{2}}}\text{
,}$

where the norm used is the $L_{2}$ norm. However in this case, since the
right-hand-side $\mathbf{b}_{i}$ changes at each Newton iteration and should
converge towards $\mathbf{0}$, the relative residual becomes too restrictive
near the solution and results in additional meaningless linear iterations. To
avoid this problem, the norm in the first iteration
$\left\|\mathbf{b}_{0}\right\|$ is used here to normalize the residual for the
rest of the Newton iterations. The convergence criterion for the linear solver
can be then rewritten as

$r_{lin}=\frac{\left\|\mathcal{A}_{i}(\delta\boldsymbol{\varepsilon})-\mathbf{b}_{i}\right\|_{L_{2}}}{\left\|\mathbf{b}_{0}\right\|_{L_{2}}}.$
(10)

Note that if the increment is linear, this expression corresponds to eq. (6).
Moreover, since $\mathbf{b}_{i}$ is zero for an equilibrated stress field, the
ratio $\|\mathbf{b}_{i}\|/\left\|\mathbf{b}_{0}\right\|$ is a relative measure
of the internal equilibrium and as a result the number of iterations required
for the linear solver decreases with the number of Newton iterations. The
second residual is the Newton residual for the non linear equation. It is
defined as the ratio between the infinity norm of the last deformation
gradient correction and the infinity norm of the change in the total strain
field within each time increment

$r_{newton}=\frac{\left\|\delta\boldsymbol{\varepsilon}\right\|_{\infty}}{\left\|\boldsymbol{\varepsilon}^{i}-\boldsymbol{\varepsilon}^{t}\right\|_{\infty}}.$
(11)

The solution is accepted only when both of the residuals (eqs. (10) and (11))
are below their respective tolerances. Note that although the choice of the
Newton forcing term has not been optimized (as proposed in [45]), it is
observed that when a consistent tangent is used the number of Newton
iterations per strain increment was quite small, always less than five.

###  2.2 Modified displacement based FFT for infinite contrast (MoDBFFT)

A modification of the displacement based FFT approach presented in [26],
called from now on MoDBFFT for brevity, is proposed here for simulating
lattice materials. The objective is to derive a method that presents accurate
results and a good convergence rate for infinite phase contrast maintaining
standard continuum discretization and derivation in the Fourier space, without
the need of using alternative discretization approaches and modified
frequencies.

The starting point of the method is the strong formulation of the conservation
of linear momentum on a periodic domain which can be divided into two sub-
domains (see Figure 1(a)) representing the two phases: $\Omega_{m}$, the
matrix, and $\Omega_{v}$, the void. The interphase between the two phases is
$\Gamma$.

(a) Two domain scheme

(b) MoDBFFT scheme

Figure 1: Schematic representation of the domains for the standard FFT and the
present approach

The real boundary value problem is defined only in the domain $\Omega_{m}$ and
corresponds to finding the displacement field $\mathbf{u}\in\Omega_{m}$ such
that

$\left\\{\begin{array}[]{c}\nabla\cdot\boldsymbol{\sigma}\left(\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)\right)=\mathbf{0}\quad\text{in}\quad\Omega_{m}\\\
\boldsymbol{\sigma}(\mathbf{x})\cdot\mathbf{n}_{\Gamma}=\mathbf{0}\ \text{in}\
\mathbf{x}\in\Gamma\\\ \boldsymbol{\sigma}\ \text{and}\
\boldsymbol{\varepsilon}\ \text{periodic}\ \text{in}\
\Omega\end{array}\right..$ (12)

Due to the periodicity in $\Omega$, the weak formulation of this problem with
free Neumann boundary conditions on $\Gamma$ is simply

$\int_{\Omega_{m}}\boldsymbol{\sigma}(\mathbf{u}):\nabla^{s}\delta\mathbf{u}\
\mathrm{d}\Omega_{m}=0\text{ ,}$ (13)

where $\mathbf{u}$ and $\delta\mathbf{u}$ are the trial and virtual
displacement fields respectively.

FFT methods need to resolve the fields in the full cell, including the points
in $\Omega_{v}$ where no material exists and fields are not defined. On the
contrary, the balance equation (13) is not defined in the entire unit cell
$\Omega$ so any compatible strain field is acceptable for the stress
equilibrium in the region $\Omega_{v}$ and the system is underdetermined.

In the new approach, the standard formulation of the two-phase domain is
modified (see Figure 1(b)) in order to extend the weak formulation to every
point of the domain $\Omega$ , including the free Neumann conditions. To this
aim, first, an artificial elastic energy density is defined in $\Omega_{v}$
which aims to prevent the indeterminacy of the displacement in that region.
This energy density depends on the square of the displacement gradient is
weighted by a numerical parameter $\alpha$, with stress dimensions and which
represents an artificial stiffness in the interior of $\Omega_{v}$. The first
variation of the total energy in $\Omega_{v}$ leads to a new term in the weak
formulation which is defined only in that region and corresponds to

$\int_{\Omega_{v}}\alpha\nabla^{s}\mathbf{u}:\nabla^{s}\delta\mathbf{u}\
\mathrm{d}\Omega_{v}\text{ ,}$ (14)

where $\alpha$ is assumed to be small. Note that when $\alpha=0$, the original
underdetermined formulation is recovered.

To extend the weak formulation of the original problem (eq. 13) in
$\Omega_{v}$ to the full domain $\Omega$, the indicator function of the voided
region $\chi_{v}(\mathbf{x})$ which defines the microstructure, is introduced
(eq. 15)

$\displaystyle\chi_{v}=\begin{cases}1&\forall\mathbf{x}\in\Omega_{v}\\\
0&\forall\mathbf{x}\in\Omega_{m}\\\ \end{cases}.$ (15)

The eqs. (13) and (14), defined in $\Omega_{m}$ and $\Omega_{v}$ respectively,
are then premultiplied by their corresponding indicator functions to be
defined in a unique domain $\Omega$, leading to

$\displaystyle\int_{\Omega}(1-\chi_{v})\boldsymbol{\sigma}(\mathbf{u}):\nabla^{s}\delta\mathbf{u}+\chi_{v}\
\alpha\ \nabla^{s}\mathbf{u}:\nabla^{s}\delta\mathbf{u}\ \mathrm{d}\Omega=0.$
(16)

Using the divergence theorem and removing boundary terms because of the
periodicity, the corresponding strong formulation of eq. (16) follows

$\nabla\cdot\left(1-\chi_{v}\left(\mathbf{x}\right)\right)\boldsymbol{\sigma}\left(\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)\right)+\nabla\cdot\chi_{v}\left(\mathbf{x}\right)\alpha\left(\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)\right)=0$

and if the chain rule is applied to the previous equation and terms are
regrouped, the result is

$\left(1-\chi_{v}\right)\nabla\cdot\boldsymbol{\sigma}+\chi_{v}\nabla\cdot\alpha\left(\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)\right)+\nabla\chi_{v}\cdot(\alpha\left(\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)\right)-\boldsymbol{\sigma})=0.$
(17)

Second, the stress free condition at the interphase $\Gamma$ of the boundary
value problem, defined in the second equation of eq. (12) should be imposed in
equation (17). To this aim, the boundary condition on $\Gamma$ is diffused
over a thin volume which is defined by the surface $\Gamma$ with an
infinitesimal thickness [30]. This translation can be done using the surface
delta function $\delta(\mathbf{x})_{\Gamma}$, defined as [46]

$\delta(\mathbf{x})_{\Gamma}=\mathbf{n}_{\Gamma}\cdot\nabla\chi_{v}(\mathbf{x})\text{
,}$ (18)

where the direction of $\mathbf{n}_{\Gamma}$ is represented in Fig. 1(b). The
gradient of the indicator function vanishes everywhere except near the surface
$\Gamma$, where it points in the normal direction [46]. Therefore, multiplying
the previous equation by the normal vector leads to

$\mathbf{n}_{\Gamma}\delta(\mathbf{x})_{\Gamma}=\nabla\chi_{v}(\mathbf{x}).$
(19)

Using the surface delta (eqs. 18 and 19), the stress free condition can be
expressed as a volume integral as

$\int_{\Gamma}\boldsymbol{\sigma}(\mathbf{x})\cdot\mathbf{n}_{\Gamma}\mathrm{d}\Gamma=\int_{\Omega}\boldsymbol{\sigma}(\mathbf{x})\cdot\mathbf{n}_{\Gamma}\delta(\mathbf{x})_{\Gamma}\mathrm{d}\Omega=\int_{\Omega}\boldsymbol{\sigma}(\mathbf{x})\cdot\nabla\chi_{v}\mathrm{d}\Omega.$
(20)

In order to apply the stress free condition, eq. (20) has to be incorporated
to eq. (17). Since for a non-vanishing $\alpha$ and continuous displacement
field in the surface $\Gamma$ it is fulfilled

$\boldsymbol{\sigma}\cdot\mathbf{n}=\alpha\left(\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)\right)\cdot\mathbf{n}\
\ \text{in}\ \Gamma$

then, the volume counterparts of any of these terms are also identical and any
of them can be set to zero in eq. (17).

$\boldsymbol{\sigma}(\mathbf{x})\cdot\nabla\chi_{v}=\alpha\left(\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)\right)\cdot\nabla\chi_{v}\
\text{in}\ \Omega.$

Choosing for simplicity
$\alpha\left(\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)\right)\cdot\nabla\chi_{v}=0$
leads to the final expression of the strong formulation

$\nabla\cdot\left(1-\chi_{v}\left(\mathbf{x}\right)\right)\boldsymbol{\sigma}\left(\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)\right)+\chi_{v}\left(\mathbf{x}\right)\nabla\cdot\alpha\nabla^{s}\mathbf{u}\left(\mathbf{x}\right)=\mathbf{0}$
(21)

together with periodic boundary conditions in $\nabla^{s}\mathbf{u}$.

In order to impose periodicity conditions in the strain field, the
displacement field is split into two contributions as

$\mathbf{u}(\mathbf{x})=\widetilde{\mathbf{u}}(\mathbf{x})+\overline{\boldsymbol{\varepsilon}}\cdot\mathbf{x}\text{
,}$ (22)

where $\widetilde{\mathbf{u}}(\mathbf{x})$, the fluctuation of the
displacement field which is periodic and has zero average, becomes the new
unknown of the problem. If this new formulation is particularized for a linear
elastic matrix using as input a general mixed macroscopic state as a
macroscopic strain $\overline{\boldsymbol{\varepsilon}}$ or macroscopic stress
$\overline{\boldsymbol{\sigma}}$, the final equation to solve the fluctuations
in the displacement reads

$\nabla\cdot\left[\left(1-\chi_{v}\left(\mathbf{x}\right)\right)\mathbb{C}\left(\mathbf{x}\right):\left(\nabla^{s}\widetilde{\mathbf{u}}\left(\mathbf{x}\right)+\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\right)\right]+\chi_{v}\left(\mathbf{x}\right)\nabla\cdot\left(\alpha\nabla^{s}\widetilde{\mathbf{u}}\left(\mathbf{x}\right)\right)=\\\
-\nabla\cdot\left(\mathbb{C}\left(\mathbf{x}\right):\overline{\boldsymbol{\varepsilon}}-\overline{\boldsymbol{\sigma}}\right)\text{
,}$ (23)

where the field $\widetilde{\mathbf{u}}$ is the unknown, that has to be solved
together with the $IJ$ components of the overall strain tensor,
$\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}$, which
are conjugate of the applied macroscopic stress. For this last unknown, an
extra equation is used linking the $IJ$ components of
$\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}$ with
the corresponding components of the macroscopic stress,

$\frac{1}{V_{\Omega}}\int_{\Omega}\left[\mathbb{C}\left(\mathbf{x}\right):\left(\nabla^{s}\widetilde{\mathbf{u}}\left(\mathbf{x}\right)+\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\right)\right]_{IJ}\mathrm{d}\Omega=-\frac{1}{V_{\Omega}}\int_{\Omega}\left[\mathbb{C}\left(\mathbf{x}\right):\overline{\boldsymbol{\varepsilon}}\right]_{IJ}\mathrm{d}\Omega+\overline{\boldsymbol{\sigma}}_{IJ}.$
(24)

As usual in FFT methods, the differential operators can be defined by their
Fourier space counterparts using spatial frequencies as

$\displaystyle\nabla^{s}(\bullet)=\mathcal{F}^{-1}\left\\{\mathcal{F}\left\\{\bullet\right\\}\otimes^{s}\boldsymbol{\xi}\right\\}=\mathcal{F}^{-1}\left\\{\frac{1}{2}\left(\mathcal{F}\left\\{\bullet\right\\}\otimes\boldsymbol{\xi}+\boldsymbol{\xi}\otimes\mathcal{F}\left\\{\bullet\right\\}\right)\right\\}$
(25)
$\displaystyle\nabla\cdot(\bullet)=\mathcal{F}^{-1}\left\\{\mathcal{F}\left\\{\bullet\right\\}\cdot\boldsymbol{\xi}\right\\}$
(26)

To solve (23) in the Fourier space, the spatial domain is discretized in a
standard regular grid (section 2.1). The result is a linear system of
equations that now is fully determined. Nevertheless, to solve the system
iteratively, the use of a preconditioner is unavoidable for a competitive
convergence rate. The linear operator
$\mathbb{M}\left(\mathbf{x}\right)(\bullet)$, proposed for the DBFFT approach
[26], is also used here as preconditioner

$\mathbb{M}(\bullet)=\left[\nabla\cdot\overline{\mathbb{C}}\
\nabla\right]^{-1}(\bullet)\text{ ,}$ (27)

where $\overline{\mathbb{C}}$ is the volume averaged stiffness tensor,

$\overline{\mathbb{C}}=\frac{1}{V_{\Omega}}\int_{\Omega}\mathbb{C}\left(\mathbf{x}\right)\mathrm{d}\Omega$
(28)

and $V_{\Omega}$ represents the volume of the entire domain $\Omega$.

The resulting equilibrium equation written in Fourier space yields

$\mathcal{F}\left\\{\left(1-\chi_{v}\right)\mathbb{C}:\left(\mathcal{F}^{-1}\left\\{\widehat{\widetilde{\mathbf{u}}}\otimes^{s}\boldsymbol{\xi}\right\\}+\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\right)\right\\}\cdot\boldsymbol{\xi}+\mathcal{F}\left\\{\chi_{v}\mathcal{F}^{-1}\left\\{\alpha\widehat{\widetilde{\mathbf{u}}}\otimes^{s}\boldsymbol{\xi}\cdot\boldsymbol{\xi}\right\\}\right\\}=-\mathcal{F}\left\\{\mathbb{C}:\overline{\boldsymbol{\varepsilon}}-\overline{\boldsymbol{\sigma}}\right\\}\cdot\boldsymbol{\xi}\text{
,}$ (29)

and the Fourier transform of the equation to impose the macroscopic stress
components $IJ$ (eq.(24)) corresponds to

$\mathcal{F}\left\\{\left(1-\chi_{v}\right)\mathbb{C}:\left(\mathcal{F}^{-1}\left\\{\widehat{\widetilde{\mathbf{u}}}\otimes^{s}\boldsymbol{\xi}\right\\}+\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\right)\right\\}\left(\mathbf{0}\right)=\mathcal{F}\left\\{\mathbb{C}:\overline{\boldsymbol{\varepsilon}}-\overline{\boldsymbol{\sigma}}\right\\}\left(\mathbf{0}\right)\text{
.}$ (30)

The equations (29,30) are linear and can be solved iteratively using a Krylov
solver. To improve the convergence, a preconditioner $\widehat{\mathbb{M}}$ is
used for eq. (29), which is the Fourier transform of the preconditioner
defined in eq. (27), and is given by

$\widehat{\mathbb{M}}(\ast)=\left[\boldsymbol{\xi}\cdot\overline{\mathbb{C}}\cdot\boldsymbol{\xi}\right]^{-1}\cdot\ast\text{
,}$ (31)

where $\ast$ represents a complex valued vector defined in the Fourier space
for all non-zero frequencies. For preconditioning eq. (30), the inverse of the
volume averaged stiffness tensor (eq. 28) is directly used. Contrary to the
Galerkin approach using standard Fourier differentiation, the linear system of
equations defined in eq. (29) is non-singular, and therefore the Conjugate
Gradient method is able to converge efficiently and provide the solution of
the system. This is a potential benefit of this approach with respect to the
modified Galerkin which relies on the use of the more memory demanding MINRES
solver.

For the implementation, the problem unknowns are joined forming a vector
composed of the fluctuating displacement field and the components of the
macroscopic strain where the stress is imposed (eq. 30),
$\left\\{\begin{array}[]{c}\widehat{\widetilde{\mathbf{u}}}\\\
\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\end{array}\right\\}$.
The left-hand side of eqs. (29) and (30) can be expressed as a linear operator
that acts over the composed vector following

$\widehat{\mathcal{A}}^{\alpha}\left(\left\\{\begin{array}[]{c}\widehat{\widetilde{\mathbf{u}}}\\\
\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\end{array}\right\\}\right)=\left\\{\begin{matrix}\mathcal{F}\left\\{\left(1-\chi_{v}\right)\mathbb{C}:\left(\mathcal{F}^{-1}\left\\{\widehat{\widetilde{\mathbf{u}}}\otimes^{s}\boldsymbol{\xi}\right\\}+\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\right)\right\\}\cdot\boldsymbol{\xi}+\mathcal{F}\left\\{\chi_{v}\mathcal{F}^{-1}\left\\{\alpha\widehat{\widetilde{\mathbf{u}}}\otimes^{s}\boldsymbol{\xi}\cdot\boldsymbol{\xi}\right\\}\right\\}\\\
\mathcal{F}\left\\{\left(1-\chi_{v}\right)\mathbb{C}:\left(\mathcal{F}^{-1}\left\\{\widehat{\widetilde{\mathbf{u}}}\otimes^{s}\boldsymbol{\xi}\right\\}+\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\right)\right\\}\left(\mathbf{0}\right)\end{matrix}\right\\}$
(32)

and the right-hand side can be written as a vector reading as

$\widehat{\mathbf{b}}^{\alpha}=\left\\{\begin{matrix}-\mathcal{F}\left\\{\mathbb{C}:\overline{\boldsymbol{\varepsilon}}-\overline{\boldsymbol{\sigma}}\right\\}\cdot\boldsymbol{\xi}\\\
\mathcal{F}\left\\{\mathbb{C}:\overline{\boldsymbol{\varepsilon}}-\overline{\boldsymbol{\sigma}}\right\\}\left(\mathbf{0}\right)\end{matrix}\right\\}$
(33)

It should be remarked that the linear operator (eq. 32) has a significantly
higher computational cost (around 1.5x) compared to the linear operator of the
Galerkin method (eq. 3) since it requires performing the additional Fourier
transforms of a vector field in the extra term. The equilibrium is reached
when a linear residual defined as

$r_{lin}=\frac{\left\|\widehat{\mathcal{A}}^{\alpha}\left(\left\\{\begin{array}[]{c}\widehat{\widetilde{\mathbf{u}}}\\\
\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\end{array}\right\\}\right)-\widehat{\mathbf{b}}^{\alpha}\right\|_{L_{2}}}{\left\|\widehat{\mathbf{b}}^{\alpha}\right\|_{L_{2}}}$
(34)

is lower than a given tolerance.

####  2.2.1 Non linear extension

In the case of non linear materials, a linearization of equation (21) is done
similarly to the Galerkin method (Section 2.1.1). The stresses and strains are
linearized following eq. (8). The non linear problem is divided into time
increments and an iterative Newton method is used at each time increment to
solve the problem. The linearization at each Newton iteration leads to a
system of equations

$\mathcal{F}\left\\{\left(1-\chi_{v}\right)\mathbb{C}^{i-1}:\left(\mathcal{F}^{-1}\left\\{\widehat{\delta\widetilde{\mathbf{u}}}\otimes\boldsymbol{\xi}\right\\}+\delta\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\right)\right\\}\cdot\boldsymbol{\xi}+\mathcal{F}\left\\{\chi_{v}\mathcal{F}^{-1}\left\\{\alpha\widehat{\delta\widetilde{\mathbf{u}}}\otimes\boldsymbol{\xi}\cdot\boldsymbol{\xi}\right\\}\right\\}\\\
=-\mathcal{F}\left\\{\left(1-\chi_{v}\right)\boldsymbol{\sigma}\left(\mathcal{F}^{-1}\left\\{\widehat{\widetilde{\mathbf{u}}}^{i-1}\otimes\boldsymbol{\xi}\right\\}+\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}^{i-1}+\overline{\boldsymbol{\varepsilon}}^{t+\Delta
t}\right)-\overline{\boldsymbol{\sigma}}^{t+\Delta
t}\right\\}\cdot\boldsymbol{\xi}\\\
-\mathcal{F}\left\\{\chi_{v}\mathcal{F}^{-1}\left\\{\alpha\widehat{\widetilde{\mathbf{u}}}^{i-1}\otimes\boldsymbol{\xi}\cdot\boldsymbol{\xi}\right\\}\right\\}$
(35)

in which the displacement correction and the average strain correction,
$\widehat{\delta\widetilde{\mathbf{u}}}$ and
$\delta\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}$,
are the unknowns and the solution for iteration $i$ is updated as

$\widetilde{\mathbf{u}}^{i}=\widetilde{\mathbf{u}}^{i-1}+\delta\widetilde{\mathbf{u}};\quad\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}^{i}=\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}^{i-1}+\delta\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}.$

In equation (35) the macroscopic prescribed strain and stress fields enter in
the definition of the first iteration as
$\widetilde{\mathbf{u}}^{0}=\widetilde{\mathbf{u}}^{t}$ and
$\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}^{0}=\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}^{t}$.
The linear equation can be translated into a linear operator applied to the
unknown
$\widehat{\mathcal{A}}^{\alpha}_{i}\left(\left\\{\begin{array}[]{c}\delta\widehat{\widetilde{\mathbf{u}}}\\\
\delta\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\end{array}\right\\}\right)$
and an independent right-hand side vector $\widehat{\mathbf{b}}^{\alpha}_{i}$,
both defined in the $i$-th iteration. At each Newton iteration, the
preconditioner given by eq. (27) is recomputed using the average tangent
stiffness.

Analogous to the Galerkin scheme, two residuals are used to solve the non
linear problem. The linear solver residual is defined as

$r_{lin}=\frac{\left\|\widehat{\mathcal{A}}_{i}^{\alpha}\left(\left\\{\begin{array}[]{c}\delta\widehat{\widetilde{\mathbf{u}}}\\\
\delta\overline{\boldsymbol{\varepsilon}}_{\overline{\boldsymbol{\sigma}}}\end{array}\right\\}\right)-\widehat{\mathbf{b}}_{i}^{\alpha}\right\|_{L_{2}}}{\left\|\widehat{\mathbf{b}}_{0}^{\alpha}\right\|_{L_{2}}}$
(36)

which normalizes the absolute error in the linear problem by the norm of the
right-hand side vector of the first Newton iteration. The second residual is
the Newton residual, which is formulated in strains and is identical to the
one used in the Galerkin FFT in eq. (11).

##  3 Geometrical adaptation

FFT methods rely on a regular discretization of a hexahedral domain in which
the cell of the lattice material is embedded. This voxelized representation
allows a simple generation of the microstructure and the direct use of
image/tomographic data but presents two disadvantages when considering an
ideal cell, especially for coarse discretizations. First, the final density
represented can differ slightly from the designed one and, second, the
voxelized representation of the struts might impact the overall behavior of
the lattice.

In order to alleviate these problems, the properties assigned to the voxels
near the strut surface can be adapted to better capture the smoothness of the
surfaces. In this section, different approaches for determining the properties
of the voxels near the interfaces are presented. The geometry definition of
the cross section of a circular strut is represented schematically in Fig. 2.

###  3.1 Plain Voxelized representation

The plain voxelized (PV) geometry approach is the most simple representation:
it assigns lattice material or empty space to a voxel based on whether the
center of that voxel is inside or outside of the lattice struts. This method
generates sharp boundaries for the struts and the resulting relative densities
can slightly differ from the designed one in the case of coarse
discretizations.

###  3.2 Phase-field smoothening

The phase-field smoothening (PFS) method consists in creating a smooth
property transition, from the lattice material properties to zero, across the
lattice interfaces. The property profile is dictated by the minimization of a
phase-field functional and is controlled by a smoothening length scale $\ell$.
The result of the phase-field minimization is a phase map $\phi(\mathbf{x})$
that assigns to each voxel of the cell a phase value between 0 (empty space)
and 1 (lattice material). Similar to damage models [47], the stress resulting
from applying the constitutive equation of the pristine material is multiplied
by the value of the phase map at that point.

Let $\omega(\mathbf{x})$ be a function which represents the exact geometry of
the lattice, and which in this case corresponds to the indicator function of
the lattice material, $\omega(\mathbf{x})=1-\chi_{v}(\mathbf{x})$ (eq. 15)
where the value of $1$ is attributed to the points belonging to the material
and $0$ to the empty space. Then, the value of $\phi(\mathbf{x})$ is obtained
as the minimizer of the functional defined by eq. (37).

$E[\phi]=\int_{\Omega}\frac{1}{2}\ell^{2}\|\nabla\phi\|^{2}+\frac{\epsilon}{2}(\phi-\omega)^{2}\mathrm{d}\Omega\text{
,}$ (37)

where $\ell$ is the characteristic length, a parameter which defines the width
of the smoothening region and $\epsilon$ is a weight which penalizes the
difference between the exact geometry $\omega$ and its smoothened counterpart
$\phi$. The functional assumes periodicity of all the fields. The first term
of the functional penalizes the smoothened area and it is modulated by the
length of the diffusion and the second term penalizes the difference between
the initial field and the smoothened one. The result of this minimization
corresponds to the solution of the partial differential equation described in
eq. (38)

$\frac{\ell}{\epsilon}^{2}\nabla^{2}\phi-\phi=-\omega$ (38)

under periodic boundary conditions in $\phi$. If the problem is discretized in
a regular grid and $\omega(\mathbf{x})$ is replaced by its discrete
counterpart, defined by the value of the indicator function at the center of
each voxel of the grid, the equation can be explicitly solved on the Fourier
space as

$\widehat{\phi}\left(\boldsymbol{\xi}\right)=\frac{1}{1-\frac{\ell^{2}}{\epsilon}\boldsymbol{\xi}\cdot\boldsymbol{\xi}}\widehat{\omega}\left(\boldsymbol{\xi}\right)\text{
,}$ (39)

where the frequency vector $\boldsymbol{\xi}$ is given in eq. (2) and all the
fields involved are periodic.

The numerical implementation of this approach is done using two different
discretization levels. In a pre-processing step a very fine grid is used with
voxel size $d_{fine}$ to discretize the real geometry and solve the phase-
field problem to obtain $\phi(\mathbf{x})$. During the simulation of the
mechanical problem a coarser grid is used, $d_{coarse}$. The ratio between
$d_{coarse}/d_{fine}$ is between 2 and 20. The fine discretization is used for
having an accurate representation of the indicator function
$\omega(\mathbf{x})$. This fine grid is used for solving eq. (38) using as
characteristic length of diffusion $\ell$ the half of the length of the voxel
in the coarse discretization ($\ell=0.5d_{coarse}$) which encompasses several
voxels of the fine grid. Finally, the resulting field $\phi(\mathbf{x})$ in
the fine grid is averaged for each voxel of the coarse discretization to
define the phase map to be used during the mechanical simulations. Note that
the values of the phase field below 5% are taken as 0, to prevent the spurious
presence of material detached from the truss. In all the phase-field
smoothening cases th e characteristic weight $\epsilon$ is taken equal to 1.

###  3.3 Voigt approaches

The third approach to smooth out the lattice surfaces is based on the Voigt
rule, following [33]. Under this approach, the stiffness of the voxels
partially occupied by the lattice material is obtained using the Voigt
homogenization approach. This rule establishes that the effective elastic
stiffness of that composite voxel is the volume average of the stiffness of
the materials present in the voxel. Therefore, since the stiffness of the
empty phase is zero, the effective stiffness corresponds to the product of the
volume fraction of lattice material in the voxel, $\phi$, multiplied by the
stiffness of the lattice material ($\mathbb{C}_{m}$)

$\mathbb{C}=\phi\ \mathbb{C}_{m}.$

The volume fraction $\phi\in[0,1]$ is then equivalent to a phase map, as the
one generated using phase-field smoothening. In [33], to compute the volume
fraction of each phase contained in the voxels partially occupied by different
phases, it was proposed the use of a subgrid to count the number of points in
the subgrid belonging to each phase. We have followed this approach first,
using the same ratio for coarse and fine grids used in PFS, and have named
this approach as Voigt fine grid (VFG). In parallel, since the surface of the
lattice is known either by its mathematical expression or by an .stl file, we
propose an alternative way to compute the phase map $\phi$ which does not
require the use of a second grid and is just based on the distance of the
voxel center to the lattice surface. This approach is named Voigt analytic
(VA) smoothening. The method assigns $\phi=1$ or $\phi=0$ to the voxels whose
centers are respectively inside or outside of the geometry considering an
offset of $\ell/2$ with respect to the boundary of the struts $\Gamma$. For
those voxels centers whose distance to the boundary is smaller than $\ell/2$ ,
the value of $\phi$ is obtained as a linear function of the signed distance of
that center to the boundary, $D$, following eq. (40),

$\phi\left(\mathbf{x}\right)=1-\frac{D+\ell/2}{\ell}\quad\text{for}\quad
D=\left\\{\begin{array}[]{c}d(\mathbf{x},\Gamma)\ \text{if}\
\mathbf{x}\in\Omega_{m}\\\ -d(\mathbf{x},\Gamma)\ \text{if}\
\mathbf{x}\in\Omega_{v}\end{array}\right.\text{ ,}$ (40)

where $d$ denotes the distance between a point and a surface.The
characteristic length considered, $\ell$, is the length of one voxel. Note
that this definition of $\phi$ corresponds exactly to the volume fraction of
lattice material in the case of a planar interface. Since the strut curvature
is normally small with reference to the voxel length, the values obtained
using eq. (40) are almost identical to the VFG in the case of fine
discretizations. This equivalence has been assessed quantitatively and the
phase map generated using the Voigt rule with a finer grid (the same used for
the phase-field smoothening) was almost identical to the one obtained by VA
(average difference below 0.2%). For clarity, both methods are only considered
for coarse grids while for finer discretization, where the results are almost
identical, only the results of VA are represented. The practical benefit of
the VA definition of $\phi$ with respect to [33] is that the smoothening is
obtained by an analytical expression using the exact geometry and not
requiring the use of a smaller grid.

###  3.4 Combined smoothening

This approach (CS) consists in applying the phase-field smoothening (Section
3.2) to the coefficients map resulting from the Voigt analytic method (Section
3.3). The combination of these two approaches will result in a very smoothened
phase map of weight coefficients that will be multiplied by the stresses
within the different algorithms.

(a) Plain voxelized representation (PV).

(b) Voigt approaches (VA and VFG).

(c) Phase-field smoothening (PFS).

(d) Combined smoothening (CS).

Figure 2: Local densities (phase maps $\phi$) near the cross section of a
strut for the different geometrical approaches. In the grey scale $\phi=1$
corresponds to black and 0 to white.

##  4 Validation for elastic materials

The numerical methods proposed to homogenize the mechanical behavior of
lattice based materials and their combination with the different geometrical
representations are studied in this section for elastic materials. In order to
evaluate the accuracy and efficiency of the two FFT methods, several numerical
tests have been carried out, and both the macroscopic result and microscopic
fields have been compared with FEM simulations. For the microscopic solution,
the relative $L_{2}$ norm of the difference between the local fields of the
solution in a given method, $f(\mathbf{x})$ compared to a reference solution
($f_{ref}(\mathbf{x})$) is used as metric of the error

$\text{Local\ diff.}\
\left[\%\right]=\frac{\left\|f-f_{ref}\right\|_{L_{2}}}{\left\|f_{ref}\right\|_{L_{2}}}\text{
.}$ (41)

###  4.1 Lattice geometry and simulation parameters

The octet-truss lattice has been selected for the numerical studies. This
structure is one of the most interesting lattice based materials since it
presents both bending-dominated and stretching-dominated responses, depending
on the strut thicknesses and the loading conditions. Relative densities
ranging from $0.5$% to $30$% are considered.

All the FFT simulations are performed using the FFTMAD code [10], to which the
new algorithms have been added. FEM simulations are done using the commercial
code ABAQUS to serve as reference solutions in order to assess the accuracy
and efficiency of FFT approaches. The FEM model consists of a geometrically
conforming mesh of quadratic tetrahedral elements (C3D10 in ABAQUS). The
element size is controlled by the number of elements occupying the lattice
strut diameter and is taken to be equal to the corresponding FFT voxel size.
An iterative solver (CG) has been selected to carry out a fair comparison with
FFT approaches which are based on the same type of iterative solver. Periodic
boundary conditions are used and introduced using multipoint linear
constraints [48].

An example of a FEM mesh of the unit cell of an octet lattice with 10%
relative density is represented in Figure 3 together with two FFT voxelized
models of the same geometry with different discretization levels

Figure 3: Octet lattice with 10% of relative density, FEM model with 15
elements per diameter and FFT discretizations with $54^{3}$ and $269^{3}$
voxels

All simulations have been done in a single node workstation Dual 10 core
Intel(R) Xeon(R) CPU E5-2630 v4 @ 2.20GHz with 64GB RAM memory. Both ABAQUS
and FFTMAD use parallelization by threading (20 threads) both for the
evaluation of the constitutive equations and for linear algebra operations.

Regarding the elastic properties of the material, the parameters correspond to
a typical polyamide PA12 with isotropic linear elastic behavior with
$E=1.7$GPa and $\nu=0.4$. The numerical tolerances for the relative errors in
the linear iterative solvers in the Galerkin FFT, MoDBFFT and FEM have been
set to $10^{-6}$.

###  4.2 Analysis of the numerical performance of FFT approaches

In this section, the convergence rate of the different adaptations of FFT
approaches for infinite phase contrast will be studied. To this aim, the
evolution of the residual of the linear iterative solvers of eq. (4) for the
adapted Galerkin FFT (with modified frequencies and use of MINRES) and of eq.
(29) for the MoDBFFT (using standard frequencies and CG) will be compared.

Although the MoDBFFT method results in a fully determined system of linear
equations with a unique solution, the well-posedness of the resulting
coefficient matrix depends on the numerical parameter $\alpha$ . The system
becomes ill-posed when $\alpha$ decreases and gets very low values with
respect to the stiffness of the material domain. On the other hand, larger
values of $\alpha$ induce an artificial stiffness in the BVP that affects the
computed effective properties. In this case, although the null traction at the
void interface is still considered explicitly by eqs. (19-21), the overall
results are slightly affected. This effect is only relevant for large values
of $\alpha$ and is due to the components of the stress tensor which are not
contained in the surface traction. Therefore, those components which are not
directly cancelled by eqs. (19-21) are indirectly affecting the overall
response.

In this study, a linear elastic test under uniaxial tension has been simulated
for different alphas in a $\overline{\rho}=0.1$ relative density octet-truss
lattice discretized using $215^{3}$ voxels (20 voxels/diameter) and PV
representation. Figure 4 represents the relative residual values for the
different linear equilibrium equations and different $\alpha$ values, being
$E$ is the Young’s modulus of the lattice material.

Figure 4: Residual evolution and effective Young’s modulus on an octet truss
lattice for the different $\alpha$s considered

Lower values of $\alpha$ in the MoDBFFT lead to worse convergence compared to
the adapted Galerkin FFT method. As a trade-off, large $\alpha$ values affect
the macroscopic properties calculated inducing an artificial stiffness, and in
Figure 4 it can be observed that the induced differences can go up to 10% in
terms of the effective stiffness. In this work, $\alpha=10^{-4}E$ is selected
as a compromise of convergence rate and accuracy of computed effective
properties. It is interesting to note that the use of discrete frequencies is
the most important ingredient for the success of the adapted Galerkin scheme.
With discrete frequencies, the use of MINRES improves the performance of the
CG version, but the reduction in computation time is only around 5%. On the
contrary, when standard Fourier discretization is used the difference in the
performance between both solvers becomes substantial.

###  4.3 Analysis of the surface smoothening approaches

The regular discretization used in FFT can lead to actual densities slightly
different from the target one, especially for coarse discretizations. The use
of different geometrical representations (section 3) to smooth out the surface
also has an impact on the actual value of the relative density of the model.
To quantify this geometrical misrepresentation, the actual relative density
considered on the RVE has been calculated for each geometrical approach as the
volume integral of the phase map, corresponding to

$\overline{\rho}=\frac{1}{N_{vox}}\sum_{i\in N_{vox}}\phi_{i}.$

A target relative density of $\overline{\rho}=0.1$ has been analyzed for the
four geometrical representations, plain voxelized approach (PV), Voigt fine
grid (VFG), Voigt analytic (VA), phase-field smoothening (PFS), and combined
smoothening (CS). A range of discretizations from $54^{3}$ to $269^{3}$ voxels
are considered, which corresponds approximately to a range from 5 to 25 voxels
per diameter in this particular case. The resulting densities for different
levels of the RVE discretization are represented in Figure 5.

Figure 5: Resulting RVE relative density for the different surface smoothening
techniques with target relative density of $10$%.

From Figure 5 it can be observed that the misrepresentation of the density is
quite limited and the maximum error is below 0.7% in all the cases. If the
discretization is refined up to $162^{3}$ then the error is reduced below 0.1%
in the worst case. The maximum deviations occur for the plain voxelized
representation. On the contrary, the Voigt approaches give a fairly good
approximation of the density for all the mesh sizes considered, always below
0.1%. In the case of VFG and PFS the small deviations from the target density
are caused by the change on the ratio between the discretization used for the
simulation and the finer one to compute the phase map. Finally, it can be
observed that the relative densities obtained using phase-field smoothening
(before thresholding) do not modify the relative density of the geometry
function used as input. This mass conservation in the phase-field smoothening
is due to the periodic boundary conditions which makes that volume integration
of eq. (38) leads to $\overline{\phi}=\overline{\omega}$. Therefore PFS gives
the same density as VFG and CS the same density as VA.

###  4.4 Accuracy of the methods

The accuracy of the macroscopic and microscopic numerical results obtained
with the FFT approaches combined with the different smoothening techniques
will be assessed for different discretization levels, ranging from $54^{3}$ to
$269^{3}$ voxels (Fig. 3). To this end, the uniaxial tensile deformation of an
octet truss lattice with a relative density of $\overline{\rho}=0.1$ is
simulated for the combinations of FFT solvers and surface smoothening. The RVE
is deformed along one of its edges and stress free conditions are imposed in
the perpendicular directions. To assess the result of the simulations, the FFT
results are compared with FEM results with an equivalent discretization level
in terms of number of elements per truss diameter.

The macroscopic strain and stress tensors are extracted from the simulation
results to obtain the effective properties. The effective Young’s modulus and
Poisson’s ration have been represented in Figure 6 together with the
corresponding FEM results. First, it is observed that the Voigt fine grid
(VFG) approach shows almost identical behavior than the Voigt analytic (VA)
and therefore, for the shake of clarity, only VA will be considered for the
rest of the discussion.

Figure 6 shows that the convergence of FFT results with the discretization is
slower than FEM results, except the combination of the Galerkin approach with
the VA that provides a solution almost independent on the grid for model sizes
greater than 10 elements per diameter ($108^{3}$ voxels). The largest
oscillations of the effective response with the discretization are obtained
for the plain voxelized models (PV) and are a direct result of the variations
in the actual relative density of the cells. The results of the phase-field
smothening (PFS) and combined smoothening (CS) converge better than the plain
voxelized representations but worse than the VFG and VA approaches separately.
Therefore, it can be concluded that using phase-field smoothening has a non-
negligible effect on the resulting stiffness for the same average density. In
the case of plain voxel approaches, the MoDBFFT method provides the same
tendency as the Galerkin approach with a small offset, being the first one
slightly stiffer. In the case of Poisson’s ratio, which is less dependent on
the cell density, the convergence with the discretization is faster and
smoother. The maximum difference between FEM and FFT results is below 0.6% for
every discretization. In all the cases, the effective FFT response converges
to the same value that coincides with the FEM results.

Figure 6: Effective Young’s modulus and Poisson’s ratio for different
discretizations.

The accuracy of the microscopic fields has also been analyzed and compared
with FEM results. In Figure 7, the microscopic stress in the loading direction
has been superposed to the deformed geometry, magnified by a factor $\times
20$, for the two most representative FFT approaches (Galerkin with Voigt
analytic smoothening and MoDBFFT with a plain voxel approach). Qualitatively
it can be observed that the deformed shapes are almost identical and the
concentrations of stress fields are very similar both in location and
intensity.

Figure 7: Local stress fields in loading direction ($\sigma_{zz}$) on the
deformed configuration (x$20$) for FEM, Galerkin FFT (VA) and MoDBFFT (PV).

To quantify this difference, the $L_{2}$ norm of the difference in the stress
in the loading direction $\sigma_{zz}$ (eq. 41) is computed with respect to
the local fields of the FEM method and the result has been represented in Fig.
8. It can be observed that, as it happened with the effective response, FFT
results converge to FEM solutions. It is also remarkable that for
discretizations finer than $15$ voxels/diameter, the differences are always
below $10$% except for those methods where the phase-field smoothening
technique is used.

Figure 8: Local stress differences in the loading direction ($\sigma_{zz}$) as
function of the discretization level.

The solutions using the plain voxels geometric approach show a slightly better
microscopic response but need a larger number of voxels to accurately predict
the overall behavior due to the density variations. It is interesting to note
that, for a PV geometrical representation, both modified Galerkin (which uses
a rotated scheme) and MoDBFFT —which has a standard discretization— provide
very similar microfields. This result indicates that the terms included in the
MoDBFFT to approach the free surface condition have a similar effect in
smoothening the response to the use of discrete frequencies in the modified
Galerkin, as it can be observed in Fig. 8. Phase-field smoothening alleviates
the noise efficiently but, for the value of $\ell$ considered here, induces
non-negligible changes in both macroscopic and microscopic responses. As a
summary, Galerkin method combined with Voigt analytic smoothening shows the
best combination of accuracy in the microfields and effective response.

###  4.5 Effect of the relative density and numerical efficiency

The effect of the relative density on the accuracy and efficiency of the
different FFT methods is studied in this section. Relative densities
$\overline{\rho}$ ranging from $0.5$% to $30$% are considered for the octet
truss lattice. A discretization of $15$ voxels(elements)/diameter is used for
every volume fraction, leading to models with different total number of
voxels. The loading case applied is uniaxial stress, which is accounted for
using macroscopic mixed boundary conditions. For comparison purposes, finite
element simulations with the same conditions are also performed for every
cell, using in this case 15 elements per diameter.

The macroscopic specific stiffness, $E/\bar{\rho}$, and Poisson’ s ratio,
$\nu$, obtained using the different approaches are represented in Figure 9 as
functions of the cell relative density. It can be observed that, in most of
the cases, the specific Young’s modulus ($E/\bar{\rho}$) is very close to the
FEM results. The maximum relative difference with respect to FEM is $10$% for
the elastic modulus in the case of the phase-field smoothening method, showing
again that although near-to-surface local fields are smoothened with this
approach, the macroscopic response is slightly altered. On the other hand, the
difference of the response obtained using the Galerkin approach with discrete
frequencies and Voigt analytic smoothening (VA) with respect to FEM is always
below 2%. The prediction of the Poisson’s ratio was very accurate for all the
methods and densities considered, with maximum differences below 1.5%.

Figure 9: Effective Young’s modulus and Poisson’s ratio time for different
relative densities.

The microscopic stresses obtained are also analyzed and compared with the FEM
counterparts. The localization of stresses and their intensity are strongly
dependent on the relative density due to the change from a stretch to a
bending dominated behavior. In all the cases the response of all the FFT
approaches considered was very similar to the FEM results. This different
behavior is also reflected in the deformation modes, which were also
accurately captured for the FFT approaches for all the densities. As an
illustration, the diagonal stress component in the loading direction,
$\sigma_{zz}$, is represented in Figure 10 for the cell with $\bar{\rho}=30\%$
superposed to the deformation of the cell with a magnification of $\times 20$.
Qualitatively, it can be observed how both stresses and deformed shape of the
cell are very similar.

Figure 10: Local stress fields in loading direction ($\sigma_{zz}$) on the
deformed configuration (x$20$) for FEM, Galerkin FFT (VA) and MoDBFFT (PV).

From a quantitative viewpoint, the $L_{2}$ norm of the difference (eq. 41)
between the FFT microscopic stress component $\sigma_{zz}$ and the FEM value
was computed and represented in Fig. 11. The norm of the difference was around
15% in most of the cases. Again, the microscopic response of PV approaches are
more near to the FEM results than the other smoothening approaches.

Figure 11: Local stress differences in the loading direction ($\sigma_{zz}$)
for different relative densities.

As stated in the study of the cell with $\bar{\rho}=0.1$ for different
discretization levels, the Galerkin FFT combined with the Voigt smoothening
shows the best compromise between the accuracy of effective and the local
properties.

Finally, the numerical performance of the FFT approaches is analyzed for the
different relative densities. The time spent on the simulations was obtained
for different FFT approaches and FEM and was represented in Figure 12. It can
be observed first that all FFT approaches were more efficient than FEM method
for relative densities greater than $\approx 5$%. Note that this comparison is
made for a particular choice of 15 elements/voxels per diameter and, if this
number were increased, better performance of the FFT solver would be expected.
Second, curves in Fig. 12 show that if the number of elements/voxel per strut
diameter is kept constant, the simulation times decrease with the density for
both FEM and FFT due to the reduction in the total number of elements/voxels
in the lattice. Nevertheless, the time reduction grows faster in FFT than in
FEM, and the FFT simulations for a relative density of $\bar{\rho}=0.3$ were 4
to 8 times faster than FEM ones. The improvement of the performance ratio
FFT/FEM with the density can be easily explained by the number of voxels of
the full RVE that belong to the interior of the lattice. For low densities
most of the voxels of the RVE belong to the empty space, not contributing to
the cell response but having to be considered for FFT operations. Therefore,
the use of FFT for densities below $\bar{\rho}<7\%$ is not competitive with
respect to FEM. On the contrary, it is remarkable that even with this strong
disadvantage, FFT becomes clearly more efficient for relative densities
exceeding 10% making the approaches here proposed very competitive for foams
and porous materials. As a conclusion, among the different FFT approaches and
smoothening techniques proposed, the Galerkin FFT combined with the Voigt
analytical smoothening is the most interesting one since it combines very
accurate results with the best numerical performance. The displacement
approach developed, the MoDBFFT, can provide smooth results which are as
accurate as the modified Galerkin, but is not competitive in terms of
efficiency since for obtaining such accurate results a small parameter
$\alpha$ is required (i.e. $\alpha=10^{-4}\ E$) and, for this value, the
number of iterations is larger than the adapted Galerkin approach.

Figure 12: Simulation time for different relative densities, using models with
around 15 elements/voxels per diameter.

##  5 Validation for non linear material response

To assess the non linear extension of the methods proposed, simulations have
been made using a Von Mises J2 plasticity model as lattice material behavior.
The elastic constants are the same as in the previous section and perfectly
plastic behavior is considered (no strain hardening), being the yield stress
$\sigma_{y}=70$MPa. Uniaxial compression tests have been carried out with a
maximum strain of $10$%. The strain is applied using a ramp divided in $20$
regular strain increments. The tolerance used for the Newton-Raphson method is
$5\cdot 10^{-3}$. The approaches selected for this study are the Galerkin with
Voigt analytic smoothening (VA) and the MoDBFFT with a plain voxelized (PV)
representation, considered as the most representative methods from previous
results. Octet truss lattices with relative densities of 10%, 20% and 30% are
studied for a fixed discretization of 15 voxels(elements)/diameter. FEM
simulations with equivalent discretization are performed for comparison
purposes.

The resulting macroscopic stress-strain curves are represented in Figure 13.
The simulations predict a large elastic region followed by a plastic regime
with a very small hardening rate. The elastic-plastic transition is smooth and
the strength reached increases with the relative density. The results of all
the simulations are very similar, being the maximum difference between FEM and
FFT results smaller than $5$% for all the densities considered.

Figure 13: Stress-strain curves with non linear material

The microscopic accumulated equivalent plastic strain field $\varepsilon_{p}$,
defined as

$\varepsilon_{p}=\int_{t}\left(\dot{\boldsymbol{\varepsilon}}^{P}:\dot{\boldsymbol{\varepsilon}}^{P}\right)^{1/2}\mathrm{d}t$

with $\dot{\boldsymbol{\varepsilon}}^{P}$ the plastic strain rate tensor, has
been represented in Figure 14 for both FFT and FEM approaches and time
corresponding to a total compressive strain of -10%. The iso-plots are
represented in the deformed cell (x2) to observe the deformed cell shape.
During the simulations, it was observed that the plastification started in the
strut joints and after plastifying these joints behave as ball joints
resulting in an almost uniaxial stretching of the struts, as it can be
observed in Figure 14. The deformed configurations and stress distributions
predicted by both FFT approaches are very similar to the FEM results.

Figure 14: Equivalent plastic strain for FEM, Galerkin FFT with Voigt analytic
smoothening and the MoDBFFT with plain voxelized approach.

For a quantitative measure of the local difference between FFT and FEM
microscopic results, the $L_{2}$ norm of the differences in the stress
component in the loading direction is represented in Fig. 15. In all the cases
the agreement is good also from a microscopic viewpoint and the differences
were always below 15%. These differences vary with the volume fraction and
also with the geometrical representation, since plastic strain is very
localized and a smoothened surface representation might affect the intensity
of the localization in those regions.

Figure 15: Local stress differences in the loading direction ($\sigma_{zz}$)
for the different relative densities with non linear material.

##  6 Application to a real 3D tomography

The principal application of the FFT framework adaptation for lattices is the
ability to analyze directly the actual topology obtained by 3D tomography or a
similar approach. In this section, the potential of the FFT framework proposed
will be shown by simulating and comparing the responses of an ideal lattice
and its real counterpart considering the fabrication defects resulting from
the additive manufacturing process. This example illustrates that full-field
simulation of the cell microstructure obtained by tomography is an extremely
powerful technique to quantify the changes in the cell response due to
fabrication defects. Under this framework, the real porosity can be considered
explicitly inside the RVE, without performing any post-processing to obtain
averaged porosities and without the need of using homogenization models to
account for the effect of the average porosity [7].

The unit cell selected is a cubic-diagonal lattice with a designed relative
density of 14.2% and manufactured in PA12 by Selective Laser Melting powder
deposition. The cell edge length is 6.2mm and the trust nominal diameter is
0.9mm. The cell was fabricated by CIRP
`https://www.cirp.de/comp/comp\_EN.php5` following standard fabrication
parameters.

The cell microstructure was analyzed by 3D tomography using a GE (Phoenix)
Nanotom 160 kV with a Hamamatsu 7942-25SK detector and nanofocus X-ray source.
The resolution of each voxel was $\approx 4\mu$m and the tomographic data
included $1551\times 1557\times 1581$ voxels. An image of the full resolution
3D tomography is shown in Figure 16. The analysis of the 3D tomography data of
the actual cell microstructure shows a volume fraction of porosity of around
3.6%. The presence of that porosity inside the struts might affect the
properties of the cell and the actual behavior can differ from the one
expected for the design geometry.

Figure 16: Left: Full resolution 3D tomography of the cell, Right: FFT model
from tomographic data.

To quantify the effect of the porosity in the elastic response, the design and
actual geometries have been subjected to uniaxial test of 1% of deformation
using the two different FFT approaches. For the design geometry, a
discretization of $25$ voxels per diameter was used ($216^{3}$ voxel RVEs). In
the case of the real tomography, the original $1500^{3}$ pixel 3D tomography
image has been compressed to $256^{3}$ voxelized model by averaging the
densities obtained from the tomography (Fig. 16, right figure). In the case of
the Galerkin FFT, the smooth map of averaged densities has been directly
introduced as the phase-map $\phi$ into the simulation. In the case of the
MoDBFFT, thresholding of the densities has been performed, distinguishing as a
material point all the local densities above a value that enforces an average
relative density equal to the measured one, and as empty region all the rest
of the points.

The effective response obtained are given in Table 1. The FFT simulations
predict a decay in the effective Young’s modulus of around 10%, a relatively
large reduction considering the low porosity volume fraction measured (3.6%).
The prediction of the overall stiffness reduction is a very interesting
characteristic of the FFT framework because it cannot be accurately obtained
using a mean field approach since the location of the porosity within the
struts influences the macroscopic response of the cell.

| Galerkin FFT (density map) | MoDBFFT (threshold map)
---|---|---
| $E\ [MPa]$ | $\nu$ | $E\ [MPa]$ | $\nu$
Design geometry | 45.76 | 0.314 | 46.25 | 0.313
3D Tomography | 41.06 | 0.292 | 43.00 | 0.291
Table 1: Effective properties extracted from uniaxial tests on cubic-diagonal
lattice design and actual geometries.

In addition to the changes of the macroscopic response, thanks to the
resolution of the local fields, the FFT analysis can be used to estimate the
microscopic fields and hot spots of the structure and, using damage
indicators, the reduction of the lattice strength due to the presence of
defects. As an example of these microscopic fields, the stress component in
the loading direction $\sigma_{zz}$ obtained in FFT simulations has been
represented in Fig. 17. It must be noted that the differences in local fields
using Galerkin FFT (using a density map as phase map) and MoDBFFT (using a
pure voxelized approach) are below 10%. In the cell with a perfect
microstucture, stress concentrates in the trust joints and varies smoothly
through the geometry of the bars. On the contrary, on the simulations with the
actual microstructure, large stress concentrations localized near bigger pores
are observed superposed to the concentrations near the joints. The maximum
stresses found with the real microstructures are approximately 50% larger than
the ones obtained with the designed cell. Due to these stress concentrations,
if the maximum local stress were taken as a rough estimation of the fracture
initiation, the real structure would fail at stress level 50% lower than the
design geometry.

Figure 17: Local stress fields in loading direction ($\sigma_{zz}$) on the
deformed configuration (x$20$) for Galerkin FFT (top) and MoDBFFT (bottom),
design and real geometries.

##  7 Conclusions

In this paper, an optimal FFT framework for the homogenization of lattice
materials has been searched and validated. The challenge was finding an FFT
approach that preserves the accuracy, good numerical performance and ability
to use images/tomographies as direct input, in the case of domains with large
regions of empty space. To this aim, two different FFT approaches able to
solve problems containing phases with zero stiffness were combined with
several approaches to smooth out the lattice surface in order to improve its
geometrical representation and reduce the noise in the microscopic solution.

Regarding FFT solvers, after a first analysis, two algorithms have been
selected as the best options for infinite phase contrast. The first one is an
adaptation of the Galerkin FFT approach using MINRES as linear solver and
modified Fourier frequencies to consider a discrete differentiation scheme,
the rotated forward approach. The second one, the MoDBFFT, is a method based
on the displacement FFT approach in [26] which eliminates the indeterminacy of
strains in the empty regions leading to a fully determined system of
equilibrium equations which allows the use of standard Fourier discretization
and differentiation. The accuracy of the two FFT solvers considered has been
validated by comparison with FEM simulations of an octet cell for several
volume fractions and discretization levels. The homogenized response of both
FFT approaches was almost identical to FEM macroscopic response for linear
elastic and elastoplastic materials, and differences in microfields were below
20%.

Regarding the surface smoothening, several approaches have been considered
based on modifying the actual stiffness of the voxels not fully embedded in
the lattice or empty space. The impact of these geometrical representations in
the effective response and local fields has been analyzed. The Voigt analytic
smoothening technique, which interpolates the stiffness of the interfacial
voxels with the distance to the real lattice surface, was the best option
since it allowed to represent exactly the relative density of the cell
allowing to use coarser grids with very accurate macroscopic response.

In terms of numerical efficiency, both FFT solvers succeed in converging in a
relatively small number of iterations considering actual zero stiffness for
the empty regions. Nevertheless, the adaptation of the Galerkin framework
convergence rate was slightly superior in all the cases. When comparing the
efficiency with respect to FEM with the same discretization in the interior of
the cell, FFT became competitive for relative densities greater than 7%. For
relative densities of 30% (70% of porosity), the simulation of this last FFT
approach was 4 to 8 times faster than FEM. As a conclusion, the modified
Galerkin approach combined with Voigt analytic smoothening was the best FFT
framework considering accuracy, numerical efficiency, and best h-convergence.

Finally, to show the real potential of the approaches presented, both FFT
frameworks are used to simulate the behavior of an actual printed lattice by
using direct 3D tomographic data as input. The simulation volume element
explicitly included the actual surface roughness and internal porosity (around
3.6%) resulting from the fabrication process. The macroscopic elastic response
was around 10% more compliant than the ideal designed geometry, and local
stress concentrations of 50% were found near large pores. As a summary, it is
shown that this technology can help to optimize the lattice fabrication
parameters as well as accurately determine the actual lattice response taking
into account the real fabrication defects.

## Acknowledgements

The authors gratefully acknowledge the support provided by the Luxembourg
National Research Fund (FNR), Reference No. 12737941 and the European Union’s
Horizon 2020 research and innovation programme for the project “Multi-scale
Optimisation for Additive Manufacturing of fatigue resistant shock-absorbing
MetaMaterials (MOAMMM)”, grant agreement No. 862015, of the H2020- EU.1.2.1. -
FET Open Programme.

## References

  * [1] X. Yu, H. Zhou, H. Liang, Z. Jiang, and L. Wu.  Mechanical metamaterials associated with stiffness, rigidity and compressibility: A brief review.  Progress in Materials Science, 94:114 – 173, 2018.
  * [2] V.S. Deshpande, N.A. Fleck, and M.F. Ashby.  Effective properties of the octet-truss lattice material.  Journal of the Mechanics and Physics of Solids, 49(8):1747 – 1769, 2001.
  * [3] S.Xu, J. Shen, S.Zhou, X. Huang, and Y. M. Xie.  Design of lattice structures with controlled anisotropy.  Materials & Design, 93:443 – 447, 2016.
  * [4] C. Findeisen, J. Hohe, M. Kadic, and P. Gumbsch.  Characteristics of mechanical metamaterials based on buckling elements.  Journal of the Mechanics and Physics of Solids, 102:151 – 164, 2017\.
  * [5] C. Lestringant and D. M. Kochmann.  Modeling of flexible beam networks and morphing structures by geometrically exact discrete beams.  Journal of Applied Mechanics, 87(8), 05 2020.  081006\.
  * [6] H. Wadley J. Berger and R. McMeeking.  Mechanical metamaterials at the theoretical limit of isotropic elastic stiffness.  Nature, 543:533–537, 2017.
  * [7] Y. Amani, S. Dancette, P. Delroisse, A. Simar, and E. Maire.  Compression behavior of lattice structures produced by selective laser melting: X-ray tomography based experimental and finite element approaches.  Acta Materialia, 159:395 – 407, 2018.
  * [8] P. Lhuissier, C. de Formanoir, G. Martin, R. Dendievel, and S. Godet.  Geometrical control of lattice structures produced by EBM through chemical etching: Investigations at the scale of individual struts.  Materials & Design, 110:485 – 493, 2016.
  * [9] A. Prakash and R.A. Lebensohn.  Simulation of micromechanical behavior of polycrystals: finite elements versus fast Fourier transforms.  Modelling and Simulation in Materials Science and Engineering, 17(6):064010, 2009.
  * [10] S. Lucarini and J. Segurado.  On the accuracy of spectral solvers for micromechanics based fatigue modeling.  Computational Mechanics, 63:365 – 382, 2019.
  * [11] R. Ma and W. Sun.  FFT-based solver for higher-order and multi-phase-field fracture models applied to strongly anisotropic brittle materials.  Computer Methods in Applied Mechanics and Engineering, 362:112781, 2020.
  * [12] F. Ernesti, M. Schneider, and T. Böhlke.  Fast implicit solvers for phase-field fracture problems on heterogeneous microstructures.  Computer Methods in Applied Mechanics and Engineering, 363:112793, 2020.
  * [13] M. Magri, S. Lucarini, G. Lemoine, L. Adam, and J. Segurado.  An FFT framework for simulating non-local ductile failure in heterogeneous materials.  Computer Methods in Applied Mechanics and Engineering, 380:113759, 2021.
  * [14] H. Moulinec and P. Suquet.  A FFT-based numerical method for computing the mechanical properties of composites from images of their microstructures.  In R. Pyrz, editor, IUTAM Symposium on Microstructure-Property Interactions in Composite Materials, pages 235–246, Dordrecht, 1995. Springer Netherlands.
  * [15] H. Moulinec and P. Suquet.  A numerical method for computing the overall response of nonlinear composites with complex microstructure.  Computer Methods in Applied Mechanics and Engineering, 157(1):69 – 94, 1998.
  * [16] D. J. Eyre and G. W. Milton.  A fast numerical scheme for computing the response of composites using grid refinement.  European Physics Journal Applied, 6(1):41–47, 1999.
  * [17] Sven Kaßbohm, Wolfgang H. Müller, and Robert Feßler.  Fourier series for computing the response of periodic structures with arbitrary stiffness distribution.  Computational Materials Science, 32(3):387–391, 2005.
  * [18] J. Zeman, J. Vondřejc, J. Novák, and I. Marek.  Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients.  Journal of Computational Physics, 229(21):8065 – 8071, 2010.
  * [19] M. Kabel, T. Böhlke, and M. Schneider.  Efficient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deformations.  Computational Mechanics, 54:1497–1514, 12 2014.
  * [20] M. Schneider.  On the Barzilai-Borwein basic scheme in FFT-based computational homogenization.  International Journal for Numerical Methods in Engineering, 118:482–494, 2019.
  * [21] D. Wicht, M. Schneider, and T. Böhlke.  Anderson-accelerated polarization schemes for fast Fourier transform-based computational homogenization.  International Journal for Numerical Methods in Engineering, 122(9):2287–2311, 2021.
  * [22] J. Vondřejc, J. Zeman, and I. Marek.  An FFT-based Galerkin method for homogenization of periodic media.  Computers & Mathematics with Applications, 68(3):156 – 173, 2014\.
  * [23] J. Zeman, T.W.J. de Geus, J. Vondřejc, R. H. J. Peerlings, and M.G.D. Geers.  A finite element perspective on non-linear FFT-based micromechanical simulations.  International Journal for Numerical Methods in Engineering, 111:903–, 09 2017.
  * [24] T.W.J. de Geus, J. Vondřejc, J. Zeman, R.H.J. Peerlings, and M.G.D. Geers.  Finite strain FFT-based non-linear solvers made simple.  Computer Methods in Applied Mechanics and Engineering, 318:412 – 430, 2017.
  * [25] S. Lucarini and J. Segurado.  An algorithm for stress and mixed control in Galerkin based FFT homogenization.  International Journal for Numerical Methods in Engineering, 119:797–805, 2019.
  * [26] S. Lucarini and J. Segurado.  DBFFT: A displacement based FFT approach for non-linear homogenization of the mechanical behavior.  International Journal of Engineering Science, 144:103131, 2019.
  * [27] J. C. Michel, H. Moulinec, and P. Suquet.  A computational scheme for linear and non-linear composites with arbitrary phase contrast.  International Journal for Numerical Methods in Engineering, 52(1‐2):139–160, 2001.
  * [28] H. Moulinec and F. Silva.  Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials.  International Journal for Numerical Methods in Engineering, 97(13):960–985, 2014.
  * [29] S. Brisard and L. Dormieux.  FFT-based methods for the mechanics of composites: A general variational framework.  Computational Materials Science, 49(3):663 – 671, 2010.
  * [30] Q.D.To and G. Bonnet.  FFT based numerical homogenization method for porous conductive materials.  Computer Methods in Applied Mechanics and Engineering, 368:113160, 2020.
  * [31] M. Schneider.  Lippmann-Schwinger solvers for the computational homogenization of materials with pores.  International Journal for Numerical Methods in Engineering, 121(22):5017–5041, 2020.
  * [32] S. Kaßbohm, W. H. Müller, and R. Feßler.  Improved approximations of Fourier coefficients for computing periodic structures with arbitrary stiffness distribution.  Computational Materials Science, 37(1):90–93, 2006.  Proceedings of the 14th International Workshop on Computational Mechanics of Materials.
  * [33] L. Gélébart and F. Ouaki.  Filtering material properties to improve FFT-based methods for numerical homogenization.  Journal of Computational Physics, 294:90–95, 2015.
  * [34] P. Shanthraj, P. Eisenlohr, M. Diehl, and F. Roters.  Numerically robust spectral methods for crystal plasticity simulations of heterogeneous materials.  International Journal of Plasticity, 66:31–45, 2015.  Plasticity of Textured Polycrystals In Honor of Prof. Paul Van Houtte.
  * [35] W. H. Müller.  Fourier transforms and their application to the formation of textures and changes of morphology in solids.  In Yehia A. Bahei-El-Din and George J. Dvorak, editors, IUTAM Symposium on Transformation Problems in Composite and Active Materials, pages 61–72. Kluwer, 1998.
  * [36] F. Willot.  Fourier-based schemes for computing the mechanical response of composites with accurate local fields.  Comptes Rendus Mécanique, 343(3):232 – 245, 2015.
  * [37] K. S. Djaka, S. Berbenni, V. Taupin, and R. A. Lebensohn.  A FFT-based numerical implementation of mesoscale field dislocation mechanics: Application to two-phase laminates.  International Journal of Solids and Structures, 184:136–152, 2020\.
  * [38] M. Schneider, F. Ospald, and M. Kabel.  Computational homogenization of elasticity on a staggered grid.  International Journal for Numerical Methods in Engineering, 105(9):693–720, 2016.
  * [39] M. Schneider, D. Merkert, and M. Kabel.  FFT-based homogenization for microstructures discretized by linear hexahedral elements.  International Journal for Numerical Methods in Engineering, 109(10):1461–1489, 2017.
  * [40] K. S. Eloh, A. Jacques, and S. Berbenni.  Development of a new consistent discrete Green operator for FFT-based methods to solve heterogeneous problems with eigenstrains.  International Journal of Plasticity, 116:1–23, 2019.
  * [41] M. Kabel, D. Merkert, and M. Schneider.  Use of composite voxels in FFT-based homogenization.  Computer Methods in Applied Mechanics and Engineering, 294:168 – 188, 2015.
  * [42] M. Suard, G. Martin, P. Lhuissier, R. Dendievel, F. Vignat, J.-J. Blandin, and F. Villeneuve.  Mechanical equivalent diameter of single struts for the stiffness prediction of lattice structures produced by Electron Beam Melting.  Additive Manufacturing, 8:124 – 131, 2015.
  * [43] Z. Chen, Y. M. Xie, X. Wu, Z. Wang, Q. Li, and S. Zhou.  On hybrid cellular materials based on triply periodic minimal surfaces with extreme mechanical properties.  Materials & Design, 183:108109, 2019.
  * [44] V. Monchiet and G. Bonnet.  A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast.  International Journal for Numerical Methods in Engineering, 89(11):1419–1436, 2012.
  * [45] D. Wicht, M. Schneider, and T. Böhlke.  On Quasi-Newton methods in fast Fourier transform-based micromechanics.  International Journal for Numerical Methods in Engineering, 121(8):1665–1694, 2020.
  * [46] R. J. Lange.  Potential theory, path integrals and the Laplacian of the indicator.  Journal of High Energy Physics, 32(11):1029–8479, 2012.
  * [47] J. Lemaitre.  A continuous damage mechanics model for ductile fracture.  Journal of Engineering Materials and Technology-transactions of The Asme, 107(1):83–89, 1985.
  * [48] J. Segurado and J. Llorca.  A numerical approximation to the elastic properties of sphere-reinforced composites.  Journal of the Mechanics and Physics of Solids, 50(10):2107 – 2121, 2002.