FFT based approaches in micromechanics: fundamentals, methods and applications

In the one-dimensional (1D) case, the continuous Fourier transform (CFT) or simply the Fourier transform of an integrable function $f(x):\mathbb{R}\to \mathbb{C}$ is written as

$$\begin{equation}\boxed{\mathcal{F}(f)=\hat{f}(\xi ){:=}{\int }_{-\infty }^{\infty }f(x){\mathrm{e}}^{-2\pi \text{i}x\xi }\enspace \mathrm{d}x} ,\end{equation} \tag{ 1 }$$

where $\xi \in \mathbb{R}$ is the frequency in Fourier space, $i=\sqrt{-1}$, and $\hat{f}(\xi )$ and $\mathcal{F}(f)$ are both the Fourier transforms of f(x).

The inverse Fourier transform of $\hat{f}(\xi )$ is defined as

$$\begin{equation}\begin{array}{c}{\mathcal{F}}^{-1}(\hat{f})=f(x){:=}{\int }_{-\infty }^{\infty }\hat{f}(\xi ){\mathrm{e}}^{2\pi \text{i}x\xi }\enspace \mathrm{d}\xi \end{array},\end{equation} \tag{ 2 }$$

where the exponential term corresponds to the Euler formula,

$$\begin{equation*}{\text{e}}^{\text{i}2\pi \theta }=\mathrm{cos}(2\pi \theta )+i\enspace \mathrm{sin}(2\pi \theta ).\end{equation*}$$

Properties

We recall some properties of the CFT that are useful to this work. Assume that f(x) and g(x) are integrable functions in $\mathbb{R}$. Then, the following properties hold:

• Linearity: for any $a,b\in \mathbb{C}$,
  $$\begin{equation*}\mathcal{F}(af(x)+bg(x))=a\hat{f}(\xi )+b\hat{g}(\xi ).\end{equation*}$$
• Integration: substituting ξ = 0 in equation (1) we get,
  $$\begin{equation*}\hat{f}(0)={\int }_{-\infty }^{\infty }f(x)\mathrm{d}x,\end{equation*}$$
  which implies that $\mathcal{F}(f)$ evaluated at the zeroth frequency corresponds to an integration over the entire domain in real space.
• Differentiation: assuming that f(x) is an absolutely continuous differentiable function in $\mathbb{R}$ and using integration by parts we have
  $$\begin{equation}\hat{{f}^{\prime }}(x)=2\pi \text{i}\xi \hat{f}(x).\end{equation} \tag{ 3 }$$
  More generally, for the nth order derivative ${f}^{(n)}(x)=\frac{{\mathrm{d}}^{n}(f(x))}{{(\mathrm{d}x)}^{n}}$, we have
  $$\begin{equation}\hat{{f}^{(n)}}(x)={(2\pi \text{i}\xi )}^{n}\hat{f}(x).\end{equation} \tag{ 4 }$$
• Convolution:
  $$\begin{align}\hfill \begin{aligned}\hfill (f\ast g)(x)& ={\int }_{-\infty }^{\infty }f(x-s)g(s)\mathrm{d}s\hfill \\ \hfill \hat{f\ast g}(\xi )& =\hat{f}(\xi )\hat{g}(\xi ).\hfill \end{aligned}\end{align} \tag{ 5 }$$
• Shift theorem: for an ${x}_{0}\in \mathbb{R}$,
  $$\begin{equation}\mathcal{F}(f(x-{x}_{0}))={\text{e}}^{-\text{i}2\pi \xi {x}_{0}}\mathcal{F}(f(x)).\end{equation} \tag{ 6 }$$
• Transform of the Dirac delta function δ(x − x₀)
  $$\begin{equation}\mathcal{F}(\delta (x-{x}_{0}))={\text{e}}^{-\text{i}2\pi \xi {x}_{0}}\end{equation} \tag{ 7 }$$
  and if x₀ = 0 then, $\hat{f}=1$.

If an integrable function f(x) is periodic in the interval [0, L] (see figure 1) such that

$$\begin{equation*}f(x)=f(x+nL)\quad \text{for}\enspace n\in \mathbb{Z},\end{equation*}$$

then the range of frequencies of the Fourier transform is discrete.

Zoom In Zoom Out Reset image size

This can be demonstrated using the definition of the inverse transform. Consider

$$\begin{equation}f(x+nL)={\int }_{-\infty }^{\infty }\hat{f}(\xi ){\text{e}}^{\text{i}2\pi \xi (x+nL)}\enspace \mathrm{d}\xi ={\int }_{-\infty }^{\infty }\hat{f}(\xi ){\text{e}}^{\text{i}2\pi \xi x}\enspace {\text{e}}^{\text{i}2\pi \xi nL}\enspace \mathrm{d}\xi .\end{equation} \tag{ 8 }$$

Because f(x) − f(x + nL) = 0, we have

$$\begin{equation}{\int }_{-\infty }^{\infty }\hat{f}(\xi ){\text{e}}^{\text{i}2\pi \xi x}(1-{\text{e}}^{\text{i}2\pi \xi nL})\mathrm{d}\xi =0\end{equation} \tag{ 9 }$$

and, since x and $\hat{f}$ are arbitrary, then the expression between the parenthesis must be zero. The only values of ξ that respect this condition are

$$\begin{equation*}\xi =k\frac{1}{L},\quad \forall \enspace k\in \mathbb{Z}.\end{equation*}$$

In this case, the integral in the inverse Fourier transform becomes a sum and f can be expressed as

$$\begin{equation}\boxed{f(x)=\sum\limits _{k=-\infty }^{\infty }\frac{1}{L}{\hat{f}}_{\hspace{-2.5pt}k}\enspace {\text{e}}^{\text{i}2\pi xk/L}} ,\end{equation} \tag{ 10 }$$

which is the general expression of a Fourier series where $(1/L){\hat{f}}_{\hspace{-2.5pt}k}$ are the Fourier coefficients which correspond to the value of the inverse Fourier transform $\hat{f}$. For a discrete frequency ξ = k/L, the corresponding coefficient is obtained following the definition of the integral Fourier transform equation (1) in the infinite spatial domain as

$$\begin{align}\hfill {\hat{f}}_{\hspace{-2.5pt}k}& =\sum\limits _{n=-\infty }^{\infty }{\int }_{0}^{L}f(x+nL){\text{e}}^{-\text{i}2\pi k(x+nL)/L}\enspace \mathrm{d}x\hfill \\ \hfill & ={\int }_{0}^{L}f(x)\enspace {\text{e}}^{-\text{i}2\pi kx/L}\left(\sum\limits _{n=-\infty }^{\infty }{\text{e}}^{-\text{i}2\pi kn}\right)\mathrm{d}x={\int }_{0}^{L}f(x){\text{e}}^{-\text{i}2\pi kx/L}\enspace \mathrm{d}x.\hfill \end{align} \tag{ 11 }$$

Note that equations (10) and (11) correspond to the exact expressions of the Fourier and inverse Fourier transforms, respectively, in the case that f is a periodic integrable function.

Let the vector f_n , n = 0, 1, ..., N − 1 be a discrete periodic function defined in an interval of length L = 1 where f_n = f(x_n ), and x_n , n = 0, 1, ..., N − 1, is a set of N equally spaced points in that interval. The discrete Fourier transform (DFT) is the counterpart of the CFT for discrete periodic functions. Its definition is

$$\begin{equation}\boxed{{\mathcal{F}}^{d}{(f)}_{k}=\hat{f}({\xi }_{k})={\hat{f}}_{\hspace{-2.5pt}k}=\sum\limits _{n=0}^{N-1}{f}_{n}\enspace {\text{e}}^{-\text{i}\frac{2\pi }{N}kn}} \end{equation} \tag{ 12 }$$

and the inverse DFT corresponds to

$$\begin{equation}\boxed{{{\mathcal{F}}^{d}}^{-1}{(\hat{f})}_{n}=f({x}_{n})={f}_{\hspace{-1.5pt}n}=\frac{1}{N}\sum\limits _{k=0}^{k=N-1}{\hat{f}}_{\hspace{-2.5pt}k}\enspace {\text{e}}^{\text{i}\frac{2\pi }{N}kn}} ,\end{equation} \tag{ 13 }$$

where the notation ${\mathcal{F}}^{d}(f)$ is used here to emphasize the difference between the DFT (equations (12) and (13)) and the CFT (equations (1) and (2)).

Since Fourier transforms can be defined for periodic functions $f:\mathbb{R}\to \mathbb{C}$ in a closed interval $[a,b]\in \mathbb{R}$, one can derive the DFT as a special case of the CFT of a discrete approximation of f. In the following, this derivation is performed in order to understand the errors that are made when substituting the CFT by the DFT in the application of a differential operator.

Consider a periodic continuous function f(x) on an interval [0, L]. Let [0, L] be a domain discretized into N equal intervals that are defined, for example, as the value of the left side of each interval: ${x}_{n}=\frac{nL}{N},n=0,1,\dots ,N-1$. An approximation of f(x), say f^N (x), can then be built as a piecewise constant function (figure 2)

$$\begin{equation*}{f}^{N}(x)=f\left({x}_{n},\enspace \text{with}\enspace n=\text{int}\left(n\frac{x}{L}\right)\right).\end{equation*}$$

Zoom In Zoom Out Reset image size

Note that in the limit case of N → ∞, we get f^N (x) → f(x). The discrete function f^N (x) can be represented by the array f_n . An example of these functions can be found in figure 2 where nine intervals are used to discretize L.

The Fourier transform of the piecewise constant function f^N (x) is given by equation (11),

$$\begin{equation}\hat{{f}^{N}}({\xi }_{k})={\hat{{f}^{N}}}_{\hspace{-2.5pt}k}={\int }_{0}^{L}{f}^{N}(x){\text{e}}^{-\text{i}2\pi xk/L}\mathrm{d}x,\quad k\in \mathbb{Z}.\end{equation} \tag{ 14 }$$

The graph includes the original periodic function f(x), its discrete representation f_n and the piecewise constant function used to approximate the Fourier integrals, f^N .

Neglecting terms |k| ⩾ N, taking the value of the integral in each interval as its value in the position of x_n multiplied by the interval length dx = L/N leads to

$$\begin{align}\mathcal{F}{({f}^{N})}_{k}& ={\hat{{f}^{N}}}_{\hspace{-2.5pt}k}={\int }_{0}^{L}{f}^{N}(x){\text{e}}^{-\text{i}2\pi xk/L}\mathrm{d}x\approx \sum\limits _{n=0}^{N-1}{f}_{\hspace{-1.5pt}n}\enspace {\text{e}}^{-\text{i}2\pi nk/N}\frac{L}{N}\\ & =\frac{L}{N}{\mathcal{F}}^{d}{(f)}_{k},\quad k\in [0,N-1].\end{align} \tag{ 15 }$$

If the normalization factor L/N is taken as one, the continuous expressions are directly replaced by their discrete counterparts.

A more refined analysis can be done (Eloh et al 2019), if the integral in equation (14) is solved exactly,

$$\begin{equation*}\hat{{f}^{N}}({\xi }_{k})={\hat{{f}^{N}}}_{\hspace{-2.5pt}k}=\sum\limits _{n=0}^{N-1}{f}_{\hspace{-1.5pt}n}{\int }_{{x}_{n}}^{{x}_{n}+L/N}{\text{e}}^{-\text{i}2\pi xk/L}\enspace \mathrm{d}x,\quad k\in \mathbb{Z}.\end{equation*}$$

Integrating previous expression, substituting x_n = nL/N, and regrouping terms results in

$$\begin{equation*}{\hat{{f}^{N}}}_{\hspace{-2.5pt}k}=\sum\limits _{n=0}^{N-1}{f}_{\hspace{-1.5pt}n}\frac{L}{N}\enspace {\text{e}}^{-\text{i}\left(\frac{2\pi k}{N}\right)\left(n+\frac{1}{2}\right)}\frac{1}{\text{i}\left(\frac{2\pi k}{N}\right)}\left[{\text{e}}^{\text{i}\left(\frac{2\pi k}{2N}\right)}-{\text{e}}^{-\text{i}\left(\frac{2\pi k}{2N}\right)}\right]k\in \mathbb{Z},\end{equation*}$$

which can be simplified using the definition of the sinc function

$$\begin{equation*}\text{sinc}(x)=\frac{\mathrm{sin}(x)}{x}=\frac{1}{\text{i}u}({\mathrm{e}}^{\mathrm{i}u/2}-{\mathrm{e}}^{\mathrm{i}u/2})\end{equation*}$$

to give

$$\begin{align}\mathcal{F}({f}^{N})& =\frac{L}{N}\enspace {\text{e}}^{-\text{i}\left(\frac{\pi k}{N}\right)}\text{sinc}\left(\frac{\pi k}{N}\right)\sum\limits _{n=0}^{N-1}{f}_{\hspace{-1.5pt}n}\enspace {\text{e}}^{-\text{i}\left(\frac{2\pi kn}{N}\right)}\\ & =\frac{L}{N}\enspace {\text{e}}^{-\text{i}\left(\frac{\pi k}{N}\right)}\text{sinc}\left(\frac{\pi k}{N}\right){\mathcal{F}}^{d}({f}_{\hspace{-1.5pt}n})\quad k\in \mathbb{Z}.\end{align}$$

Comparing this expression with the standard discrete approximation equation (15) it can be observed that here, in addition to considering infinite frequencies, the DFT is weighted by the sinc function. This exact expression of the Fourier transform of a piecewise constant function can be used to define periodic differential operators, which reduce the noise associated to Gibbs phenomena or aliasing effect (Eloh et al 2019), but in majority of the Fourier homogenization approaches, the DFT is directly used to compute the CFT, as done in equation (15).

At this point, it is important to note the difference between the Fourier series and the DFT. The Fourier series is the exact Fourier transform of a continuous periodic function and corresponds to an additive decomposition in an infinite number of frequencies k, each with its corresponding weight equations (10) and (11). The Fourier series is convergent i.e. Fourier coefficients, or weights, tend to zero for k → ∞ (Lejeune-Dirichlet 1829) and the original function is only recovered exactly when summing infinite terms. Truncating the Fourier series to a finite number of frequencies allows to obtain an approximation of the original function and the accuracy depends on the number of terms used in the series. The DFT (equation (12)) aims to express a discrete periodic function also as a decomposition in frequencies but using a finite number of them, where this number is equal to the number of points in the discrete function. Therefore, there is no truncation in the definition of the DFT and it corresponds to a bijective application such that

$$\begin{equation*}{f}_{n}={\mathcal{F}}^{-1}(\mathcal{F}({f}_{n})).\end{equation*}$$

Note that, for k = 1, the angle $\theta =\frac{2\pi }{N}n$ moves from 0 to 2π. Due to periodicity of the trigonometrical functions, to minimize aliasinig effect, the interval is usually translated to [−π, π]. In this case, x ∈ [−L/2, L/2] and ξ ∈ [−N/2, N/2]. Furthermore, in most numerical solvers, the frequencies are reordered in the following manner

$$\begin{equation}\mathbf{q}=\begin{cases}0,1,2,\dots ,\frac{N-1}{2},-\frac{N-1}{2},-\frac{N-1}{2}+1,\dots ,-1\quad \hfill & \text{if}\enspace N\enspace \text{even}\hfill \\ 0,1,2,\dots ,\frac{N}{2},-\frac{N}{2}+1,-\frac{N}{2}+2,\dots ,-1\quad \hfill & \text{if}\enspace N\enspace \text{odd}\hfill \end{cases}\end{equation} \tag{ 16 }$$

q_n = q[n]. In FFT packages, the frequencies in the DFT are usually defined between $\left[-\frac{1}{2},\frac{1}{2}\right]$ as

$$\begin{equation}{\xi }_{n}={q}_{n}\frac{1}{N}.\end{equation} \tag{ 17 }$$

Then, the DFT can be defined using the reordered frequencies as

$$\begin{equation}\hat{f}(\xi )={\hat{f}}_{\hspace{-2.5pt}k}=\sum\limits _{n=0}^{N-1}{f}_{\hspace{-1.5pt}n}\enspace {\mathrm{e}}^{-\text{i}\frac{2\pi }{N}{q}_{k}{q}_{n}}.\end{equation} \tag{ 18 }$$

The DFT by definition is computationally very expensive, $\mathcal{O}({n}^{2})$, and in practice it is never computed using the above equation; a deeper analysis of its computation is provided in section 2.1.4.

Properties of the DFT

The properties of the DFT that are interesting for solid mechanics applications are

• Periodicity
  $$\begin{equation*}{\hat{f}}_{k+N}={\hat{f}}_{\hspace{-2.5pt}k}.\end{equation*}$$
• Transform of a real function. If f is the array representing a real periodic function, then
  $$\begin{equation*}\hat{f}(\xi )={\hat{f}}^{\ast }(-\xi )\end{equation*}$$
  so, for even N, only (N − 1)/2 + 1 terms have to be computed.
• Linearity
  $$\begin{equation*}\mathcal{F}(a\boldsymbol{f}+b\boldsymbol{g})=a\hat{\boldsymbol{f}}+b\hat{\boldsymbol{g}}.\end{equation*}$$
• Integration
  $$\begin{equation*}\hat{f}(0)=\sum\limits _{n=0}^{N}{f}_{n}.\end{equation*}$$
• Differentiation. If f is defined in a periodic interval L
  $$\begin{equation}{\hat{{f}^{\prime }}}_{\hspace{-2.5pt}k}=\frac{\mathrm{d}{\hat{f}}_{\hspace{-2.5pt}k}}{\mathrm{d}x}=\frac{\mathrm{d}{\hat{f}}_{\hspace{-2.5pt}k}}{\mathrm{d}n}\frac{\mathrm{d}n}{\mathrm{d}x}=\text{i}2\pi \frac{{q}_{k}}{N}\frac{N}{L}{\hat{f}}_{\hspace{-2.5pt}k}=\text{i}2\pi {\xi }_{k}\frac{N}{L}{\hat{f}}_{\hspace{-2.5pt}k}\end{equation} \tag{ 19 }$$
  with ${\xi }_{k}\in [-\frac{1}{2},\frac{1}{2}]$.
• (Circular) Convolution
  $$\begin{equation}{(\hspace{2.3pt}\boldsymbol{f}\ast \boldsymbol{g})}_{n}=\sum\limits _{j=0}^{N-1}{f}_{n-j}\cdot {g}_{j}\end{equation} \tag{ 20 }$$
  $$\begin{equation}{\hat{\mathbf{f}\ast \mathbf{g}}}_{k}={\hat{f}}_{\hspace{-2.5pt}k}\cdot {\hat{g}}_{k}.\end{equation} \tag{ 21 }$$
• Shift theorem
  $$\begin{equation*}{\hat{f}}_{n+{n}_{0}}={\text{e}}^{-\text{i}\frac{2\pi }{N}k{n}_{0}}{\hat{f}}_{n}.\end{equation*}$$

Matrix representation of the DFT

As stated before, the DFT (equation (18)) is a linear transformation of a discrete function (an array of values), therefore, it can be written in matrix form as

$$\begin{equation}\hat{\mathbf{f}}=\mathbf{F}\cdot \mathbf{f}=\left[\begin{matrix}{\omega }^{0\cdot 0}& {\omega }^{0\cdot 1}& \dots & {\omega }^{0\cdot N-1}\\ {\omega }^{1\cdot 0}& {\omega }^{1\cdot 1}& \dots & {\omega }^{1\cdot N-1}\\ {\vdots}& {\vdots}& {\vdots}& {\vdots}\\ {\omega }^{N-1\cdot 0}& {\omega }^{N-1\cdot 1}& \dots & {\omega }^{N-1\cdot N-1}\end{matrix}\right]\left.\begin{cases}{f}_{\hspace{-2pt}0}\\ {f}_{\hspace{-2pt}1}\\ {\vdots}\\ {f}_{\hspace{-1.5pt}n}\end{cases}\right\},\end{equation} \tag{ 22 }$$

where

$$\begin{equation*}{\omega }^{k\cdot n}={\left({\mathrm{e}}^{-\text{i}2\pi \frac{1}{N}}\right)}^{k\cdot n}.\end{equation*}$$

The resulting matrix F is a Vandermonde square matrix, fully dense and invertible. In addition, the periodicity of the Euler formula,

$$\begin{equation*}{\omega }^{k\cdot n}={\omega }^{\mathrm{mod}(k\cdot n/N)},\end{equation*}$$

where mod() is the remainder of the operation, allows to simplify the transformation matrix as

$$\begin{equation}\mathbf{F}(\omega )=\left[\begin{array}{lllll} 1 & 1 & 1 & \dots & 1 \\ 1 & {\omega }^{1} & {\omega }^{2} & \dots & {\omega }^{N-1} \\ 1 & {\omega }^{2} & {\omega }^{4} & \dots & {\omega }^{2(N-1)} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ 1 & {\omega }^{N-1} & {\omega }^{2(N-1)} & \dots & {\omega }^{(N-1)(N-1)} \end{array}\right].\end{equation} \tag{ 23 }$$

The inverse DFT in its matrix form is then given by

$$\begin{equation*}{\mathbf{F}}^{-1}=\frac{1}{N}\mathbf{F}(-\omega ).\end{equation*}$$

As previously mentioned, the DFT is a very expensive algorithm ($\mathcal{O}({n}^{2})$) since it requires n operations for each of the n Fourier components, which is also the cost of the product of the DFT matrix by an array. However, this cost can be drastically reduced using the symmetries of the transformation matrix, as it was first proposed in Cooley and Tukey (1965), giving rise to the so-called FFT algorithm.

In order to present the FFT algorithm, consider the following example. Let N = 4 and let f be a discrete function with four values f₀, f₁, f₂ and f₃. Following the notation of previous section ω corresponds to

$$\begin{equation*}\omega ={\mathrm{e}}^{-\text{i}2\pi /4}={\mathrm{e}}^{-\text{i}\pi /2}.\end{equation*}$$

Note that due to periodicity, all powers of ω greater than 2 can be expressed as the negative of lower powers, in particular

$$\begin{equation*}{\omega }^{3}=-{\omega }^{1},\qquad {\omega }^{2}=-{\omega }^{2}=-1,\qquad {\omega }^{4}=1.\end{equation*}$$

The full transformation is given by

$$\begin{equation*} \left.\left\{\begin{matrix} {\hat{f}}_{\hspace{-2.5pt}0} \\ {\hat{f}}_{\hspace{-2.5pt}1} \\ {\hat{f}}_{\hspace{-2.5pt}2} \\ {\hat{f}}_{\hspace{-2.5pt}3} \end{matrix}\right.\right\}=\left[\begin{matrix} 1 & 1 & 1 & 1 \\ 1 & \omega & {\omega }^{2} & {\omega }^{3} \\ 1 & {\omega }^{2} & {\omega }^{4} & {\omega }^{6} \\ 1 & {\omega }^{3} & {\omega }^{6} & {\omega }^{9} \end{matrix}\right]\left.\left\{\begin{matrix} {f}_{\hspace{-2pt}0} \\ {f}_{\hspace{-2pt}1} \\ {f}_{\hspace{-2pt}2} \\ {f}_{\hspace{-2pt}3} \end{matrix}\right.\right\}=\left[\begin{matrix} 1 & 1 & 1 & 1 \\ 1 & \omega & {\omega }^{2} & -\omega \\ 1 & {\omega }^{2} & 1 & {\omega }^{2} \\ 1 & -\omega & {\omega }^{2} & \omega \end{matrix}\right]\left.\left\{\begin{matrix} {f}_{\hspace{-2pt}0} \\ {f}_{\hspace{-2pt}1} \\ {f}_{\hspace{-2pt}2} \\ {f}_{\hspace{-2pt}3} \end{matrix}\right.\right\}.\end{equation*}$$

The key idea is splitting the data to reduce the number of operations. f is split into two arrays, f_even = [f₀, f₂] and f_odd = [f₁, f₃]. Then, these arrays are transformed by DFTs of size N/2 = 2

$$\begin{equation*}{\hat{\mathbf{f}}}_{\text{even}}=\left.\left\{\begin{array}{c}{\hat{f}}_{\hspace{-2.5pt}\mathrm{e}\mathrm{v},0}\\ {\hat{f}}_{\hspace{-2.5pt}\mathrm{e}\mathrm{v},1}\end{array}\right.\right\}=\left[\begin{array}{cccc}1& 1\\ 1& {\omega }_{2}\end{array}\right]\left.\left\{\begin{array}{c}{f}_{\hspace{-2pt}0}\\ {f}_{\hspace{-2pt}2}\end{array}\right.\right\}=\left.\left\{\begin{array}{c}{f}_{\hspace{-2pt}0}+{f}_{\hspace{-2pt}}\\ {f}_{\hspace{-2pt}0}+{\omega }_{2}{f}_{\hspace{-2pt}2}\end{array}\right.\right\}\end{equation*}$$

$$\begin{equation*}{\hat{\mathbf{f}}}_{\text{odd}}=\left.\left\{\begin{matrix} {\hat{f}}_{\hspace{-2.5pt}\mathrm{o}\mathrm{d}\mathrm{d},0} \\ {\hat{f}}_{\hspace{-2.5pt}\mathrm{o}\mathrm{d}\mathrm{d},1} \end{matrix}\right.\right\}=\left[\begin{matrix} 1 & 1 \\ 1 & {\omega }_{2} \end{matrix}\right]\left.\left\{\begin{matrix} {f}_{\hspace{-2pt}1} \\ {f}_{\hspace{-2pt}3} \end{matrix}\right.\right\}=\left.\left\{\begin{matrix} {f}_{\hspace{-2pt}1}+{f}_{\hspace{-2pt}3} \\ {f}_{\hspace{-2pt}1}+{\omega }_{2}{f}_{\hspace{-2pt}3} \end{matrix}\right.\right\}\end{equation*}$$

being ω₂ = ω² = e^−iπ = −1. Now the objective is obtaining $\hat{\mathbf{f}}$ from ${\hat{\mathbf{f}}}_{\text{even}}$ and ${\hat{\mathbf{f}}}_{\text{odd}}$

$$\begin{equation*}{\hat{f}}_{\hspace{-2.5pt}0}={\omega }^{0\cdot 0}{f}_{\hspace{-2pt}0}+{\omega }^{0\cdot 1}{f}_{\hspace{-2pt}1}+{\omega }^{0\cdot 2}{f}_{\hspace{-2pt}2}+{\omega }^{0\cdot 3}{f}_{\hspace{-2pt}3}={\hat{f}}_{\mathrm{e}\mathrm{v},0}+{\hat{f}}_{\hspace{-2.5pt}\mathrm{o}\mathrm{d}\mathrm{d},0}={\hat{f}}_{\hspace{-2.5pt}\mathrm{e}\mathrm{v},0}+{\omega }^{0}{\hat{f}}_{\hspace{-2.5pt}\mathrm{o}\mathrm{d}\mathrm{d},0}\end{equation*}$$

$$\begin{align*} {\hat{f}}_{\hspace{-2.5pt}1}& ={\omega }^{1\cdot 0}{f}_{\hspace{-2pt}0}+{\omega }^{1\cdot 1}{f}_{\hspace{-2pt}1}+{\omega }^{1\cdot 2}{f}_{\hspace{-2pt}2}+{\omega }^{1\cdot 3}{f}_{\hspace{-2pt}3}=({f}_{\hspace{-2pt}0}+{\omega }^{2}{f}_{\hspace{-2pt}2})+(\omega {f}_{\hspace{-2pt}1}+{\omega }^{1\cdot 3}{f}_{\hspace{-2pt}3}) \\ & =({f}_{\hspace{-2pt}0}+{\omega }_{2}{f}_{\hspace{-2pt}2})+\omega ({f}_{\hspace{-2pt}1}+{\omega }^{2}{f}_{\hspace{-2pt}3})=({f}_{\hspace{-2pt}0}+{\omega }_{2}{f}_{\hspace{-2pt}2})+\omega ({f}_{\hspace{-2pt}1}+{\omega }_{2}{f}_{\hspace{-2pt}3})={\hat{f}}_{\hspace{-2.5pt}\mathrm{e}\mathrm{v},1}+{\omega }^{1}{\hat{f}}_{\hspace{-2.5pt}\mathrm{o}\mathrm{d}\mathrm{d},1} \end{align*}$$

$$\begin{align*}{\hat{f}}_{\hspace{-2.5pt}2}& ={\omega }^{2\cdot 0}{f}_{\hspace{-2pt}0}+{\omega }^{2\cdot 1}{f}_{\hspace{-2pt}1}+{\omega }^{2\cdot 2}{f}_{\hspace{-2pt}2}+{\omega }^{2\cdot 3}{f}_{\hspace{-2pt}3}=({f}_{\hspace{-2pt}0}+{\omega }^{4}{f}_{\hspace{-2pt}2})+({\omega }^{2}{f}_{\hspace{-2pt}1}+{\omega }^{6}{f}_{\hspace{-2pt}3})\\ & =({f}_{\hspace{-2pt}0}+{f}_{\hspace{-2pt}2})+{\omega }^{2}({f}_{\hspace{-2pt}1}+{\omega }_{2}^{4}{f}_{\hspace{-2pt}3})=({f}_{\hspace{-2pt}0}+{f}_{\hspace{-2pt}2})-1({f}_{\hspace{-2pt}1}+{f}_{\hspace{-2pt}3})={\hat{f}}_{\hspace{-2.5pt}\mathrm{e}\mathrm{v},0}-{\omega }^{0}{\hat{f}}_{\hspace{-2.5pt}\mathrm{o}\mathrm{d}\mathrm{d},0}\end{align*}$$

$$\begin{align*}{\hat{f}}_{\hspace{-2.5pt}3}& ={\omega }^{3\cdot 0}{f}_{\hspace{-2pt}0}+{\omega }^{3\cdot 1}{f}_{\hspace{-2pt}1}+{\omega }^{3\cdot 2}{f}_{\hspace{-2pt}2}+{\omega }^{3\cdot 3}{f}_{\hspace{-2pt}3}=({f}_{\hspace{-2pt}0}+{\omega }^{6}{f}_{\hspace{-2pt}2})+({\omega }^{3}{f}_{\hspace{-2pt}1}+{\omega }^{9}{f}_{\hspace{-2pt}3})\\ & =({f}_{\hspace{-2pt}0}+{\omega }^{2}{f}_{\hspace{-2pt}2})+{\omega }^{3}({f}_{\hspace{-2pt}1}+{\omega }^{6}{f}_{\hspace{-2pt}3})=({f}_{\hspace{-2pt}0}+{\omega }_{2}{f}_{\hspace{-2pt}2})\omega ({f}_{\hspace{-2pt}1}+{\omega }_{2}{f}_{\hspace{-2pt}3})={\hat{f}}_{\hspace{-2.5pt}\mathrm{e}\mathrm{v},1}-{\omega }^{1}{\hat{f}}_{\hspace{-2.5pt}\mathrm{o}\mathrm{d}\mathrm{d},1}.\end{align*}$$

The DFT of size N = 4 is then obtained combining the two smaller DFTs. The computational cost of the direct DFT was 4^*4 = 16 operations. The FFT approach implies 2 times 2^*2 operations instead.

Generalizing the result

$$\begin{equation}{\hat{f}}_{\hspace{-2.5pt}k}={\hat{f}}_{\hspace{-2.5pt}\mathrm{e}\mathrm{v},k}+{\omega }_{N}^{k}{\hat{f}}_{\hspace{-2pt}\mathrm{o}\mathrm{d}\mathrm{d},k}\quad k=0,1,\dots ,N/2\end{equation} \tag{ 24 }$$

$$\begin{equation}{\hat{f}}_{k+N/2}={\hat{f}}_{\hspace{-2.5pt}\mathrm{e}\mathrm{v},k}-{\omega }_{N}^{k}{\hat{f}}_{\hspace{-2pt}\mathrm{o}\mathrm{d}\mathrm{d},k}\quad k=0,1,\dots ,N/2.\end{equation} \tag{ 25 }$$

This idea can be exploited recursively where each DFT of size N is computed by dividing it into two smaller DFTs of size (N/2), until N = 2.

In summary, the DFT computed using the definition in equation (12) or the DFT matrix in equation (22) is computationally very expensive: the cost grows as $\mathcal{O}({N}^{2})$, but the FFT algorithm (Cooley and Tukey 1965) allows to compute the DFT of an array with a computational cost of $\mathcal{O}(N\enspace {\mathrm{log}}_{2}\enspace N)$, where N is usually taken as a power of 2. Finally, note that FFT is just a clever way of computing the DFT, not an alternative transformation, and it does not imply any approximation in the computation of the DFT. Nevertheless, DFT without the FFT algorithm would not be useful for numerical methods and therefore, it is the FFT algorithm, and not the DFT, that gives name to the numerical approaches.

Thus far, all the transformations and algorithms were restricted to 1D. However, if we want to use the FFT approach to solve PDEs or to perform homogenization, then we need to define the transformations in 2-dimensional (2D) or 3-dimensional (3D) space because the functions involved are defined in 2D or 3D spaces. In the following, the definitions and developments will be done directly for the case of interest: real periodic functions in a discrete form.

The 2D DFT is defined as a sequential DFT on each axis. Let f be a real function defined in a rectangular domain $\left[\frac{-{L}_{x}}{2},\frac{{L}_{x}}{2}\right]\times \left[\frac{-{L}_{y}}{2},\frac{{L}_{y}}{2}\right]$

$$\begin{equation*}f:\left[-\frac{{L}_{x}}{2},\frac{{L}_{x}}{2}\right]\times \left[-\frac{{L}_{y}}{2},\frac{{L}_{y}}{2}\right]\to \mathbb{R}.\end{equation*}$$

The discrete form of f can be defined as the value of the function in a regular grid with N_x equi-spaced points in x and N_y points in y,

$$\begin{equation*}{f}_{{n}_{x},{n}_{y}}=f[x({n}_{x}),y({n}_{y})]\end{equation*}$$

and 0 ⩽ n_x < N_x , 0 ⩽ n_y < N_y .

The multidimensional DFT is then defined as

$$\begin{equation}\boxed{{\hat{f}}_{{k}_{x},{k}_{y}}=\sum\limits _{{n}_{y}=0}^{{N}_{y}-1}\left(\sum\limits _{{n}_{x}=0}^{{N}_{x}-1}{f}_{{n}_{x},{n}_{y}}\enspace {\mathrm{e}}^{-\text{i}2\pi {q}_{{k}_{x}}{q}_{{n}_{x}}\frac{1}{{N}_{x}}}\right){\mathrm{e}}^{-\text{i}2\pi {q}_{{k}_{y}}{q}_{{n}_{y}}\frac{1}{{N}_{y}}}} \end{equation} \tag{ 26 }$$

with 0 ⩽ k_x < N_x , 0 ⩽ ky < Ny. Equivalent to the 1D case, an N-dimensional frequency array ξ can be defined,

$$\begin{equation*}{\boldsymbol{\xi }}_{{k}_{x},{k}_{y}}=\boldsymbol{\xi }\left[{k}_{x},{k}_{y}\right]=\left[\frac{1}{{N}_{x}}{q}^{{N}_{x}}[{k}_{x}],\frac{1}{{L}_{y}}{q}^{{N}_{y}}[{k}_{y}]\right]\end{equation*}$$

and

$$\begin{equation} {q}^{N}=\begin{cases}\left[0,1,2,\dots ,\frac{N-1}{2},-\frac{N-1}{2},-\frac{N-1}{2}+1,\dots ,-1\right]\quad & \text{if}\enspace N\enspace \text{even} \\ \left[0,1,2,\dots ,\frac{N}{2},-\frac{N}{2}+1,-\frac{N}{2}+2,\dots ,-1\right]\quad & \text{if}\enspace N\enspace \text{odd} \end{cases}\end{equation} \tag{ 27 }$$

the values of this array for all the k_x and k_y values, (N_x ⋅ N_y discrete values) correspond to the Fourier points.

The multi-dimensional FFT fulfills the same properties than its 1D counterpart. In particular, is interesting to recall the partial derivative in a direction x

$$\begin{equation}\mathcal{F}{\left(\frac{\partial f}{\partial x}\right)}_{{k}_{x},{k}_{y}}=\text{i}2\pi \frac{1}{{L}_{i}}{q}^{{N}_{x}}[{k}_{x}]{\hat{f}}_{{k}_{x},{k}_{y}}=\text{i}2\pi {\xi }_{{k}_{x}}\frac{{N}_{x}}{{L}_{x}}{\hat{f}}_{{k}_{x},{k}_{y}}\end{equation} \tag{ 28 }$$

that can be directly extended to three dimensions and combined to form the Fourier form of differential operators as the grad, div or rot.

FFT packages

FFT is one of the most used algorithms in a wide range of fields including image analysis, signal processing and resolution of PDEs. Consequently, several algorithms have been developed that either improve the performance of the Cooley and Tukey algorithm or extend its application to odd and non-power-of-two discrete functions. A review of these algorithms is out of the scope of this paper; an interested reader could refer to the work of Duhamel and Vetterli (1990).

From a practical viewpoint, several implementations of the FFT algorithms can be found; most of these implementations are open-source. The algorithms provided in these packages are highly optimized in terms of memory and number of operations and also take advantage of parallel programming using either multi-threading, distributed MPI or GPU based parallelization. Among all the packages available, probably the most extended one is the so-called FFTW (abbreviation for 'the fastest Fourier transform in the West') http://fftw.org. Developed by Frigo and Johnson (Frigo and Johnson 2005), FFTW is a C library and includes DFT functions in 1D, 2D, 3D or more dimensions and of both real and complex data. Classically, FFTW was parallelized in threads and provides a very good scaling in shared memory systems. From version 3.3.1 onwards, the code also includes MPI parallel distribution. Another package of interest is fftMPI (Plimpton 2018) developed in Sandia Laboratory by Plimpton and co-workers https://fftmpi.sandia.gov. The package performs FFTs in parallel for distributed memory systems, where the FFT grid is distributed across processors. Another package for massive parallelization of three dimensional FFTs using MPI is the P3DFFT, https://p3dfft.net. This project was initiated at San Diego Supercomputer Center (SDSC) at UC San Diego by Pekurovsky (Pekurovsky 2012). The library include codes developed both in Fortran and C++. Finally, GPU acceleration has also been applied to the FFT algorithms since their structure perfectly fits the requirements for an efficient implementation. In this line, CuFFT (NVidia 2021) is an FFT package developed for Cuda https://docs.nvidia.com/cuda/cufft/index.html, the parallel programming platform of Nvidia. The package allows to compute the FFT of almost general input functions with dimension given by product of powers of the first prime numbers—best performance as usual is obtained using powers of 2—using both single and double precision.

The Green's function of the linear elastic problem in a homogeneous medium appears in almost every homogenization framework because it is the basis of Eshelby's approach (Eshelby 1957), which is at the core of any micromechanical homogenization approach. In the case of FFT based homogenization, Green's function is also behind the formulation of the method, as it will be presented in section 3.

The particularization of the Green's function method for solving the elastic problem with body-forces is presented now. Let a homoge(neous medium characterized by the stiffness tensor $\mathbb{C}$ be subjected to a body force per unit volume field b(x). Then, the stress equilibrium PDE comprising the deformation compatibility condition is written as

$$\begin{equation*}\nabla \cdot \mathbb{C}:{\nabla }^{s}\mathbf{u}(\mathbf{x})+\mathbf{b}(\mathbf{x})=0,\end{equation*}$$

where u(x) is the displacement field. The linear operator acting on the function u(x) corresponds to

$$\begin{equation*}\mathcal{L}(\cdot )=\nabla \cdot \mathbb{C}:{\nabla }^{s}(\cdot )\end{equation*}$$

and the source term f(x) in equation (30) corresponds to −b(x).

The Green's function for this problem is a second order tensor function which provides the displacement caused by a force F located at x' in a point x,

$$\begin{equation}\mathbf{u}(\mathbf{x})=\mathbf{G}({\mathbf{x},\mathbf{x}}^{\prime })\cdot \mathbf{F}.\end{equation} \tag{ 35 }$$

In the case of a boundary value problem in a domain ${\Omega}\in {\mathbb{R}}^{3}$ and a body force field b(x), the solution of the problem is

$$\begin{equation}\mathbf{u}(\mathbf{x})={\int }_{{\Omega}}\mathbf{G}(\mathbf{x},{\mathbf{x}}^{\prime })\cdot \mathbf{b}({x}^{\prime })\mathrm{d}{\mathbf{x}}^{\prime }=\mathbf{G}\ast \mathbf{b}\end{equation} \tag{ 36 }$$

where the Green's function are particularized for the body shape and boundary conditions. An sketch of the procedure for obtaining u(x) using Green's function under a discrete force field F(x₁), F(x₂) etc, or a continuous force density field b(x) is depicted in figure 3.

Zoom In Zoom Out Reset image size

In particular, for a translation invariant problem where the effect of a unit force depends only on the relative position between x and x', the Green's function is simplified to

$$\begin{equation*}G(x,{x}^{\prime })=G(x-{x}^{\prime }).\end{equation*}$$

This case corresponds to a homogeneous material and an infinite or periodic domain. In this case, analytical expressions of the function can be derived (Mura 1987) to solve the PDE that results in particularizing the definition of G equation (30) to the elastic PDE. The Green's function G_km is the solution of this PDE

$$\begin{equation}{C}_{ijkl}{G}_{km,lj}(\mathbf{x}-{\mathbf{x}}^{\prime })+{\delta }_{im}\delta (\mathbf{x}-{\mathbf{x}}^{\prime })=0.\end{equation} \tag{ 37 }$$

Equation (37) can be easily solved in Fourier space for a periodic medium, providing a rather simple closed-form expression of the Green's function. Also closed-form expressions can be found in real space for an elastic isotropic infinite body, expressed as a function of the Lamé coefficients λ and μ (Mura 1987).

The basic scheme proposed in Moulinec and Suquet (1994, 1998) solves the implicit Lippmann–Schwinger equation (50) using the fixed-point iterative method (Picard's iterations). For the numerical resolution of the Lippmann–Schwinger equation which implies the calculation of direct and inverse Fourier transforms and a Green's operator, the spatial and frequency domain are discretized, allowing the use of DFTs. The spatial and frequency domain discretization can be found in section 2.1.3, equation (16).

For the basic scheme, the discretization of the problem relies on the DFTs that may be interpreted as a trigonometric collocation discretization of the fields, where the trigonometric polynomial of order N (N is the discretization) are the interpolation functions (Nguyen et al 2011, Schneider 2015). In the case of an odd number of voxels, to prevent the loss of symmetry, the value of the gamma operator is redefined as the inverse of the reference stiffness (Moulinec and Suquet 1998).

The resulting algorithm is given in algorithm 1 where the definition of a reference medium ${\mathbb{C}}^{0}$ is required as input data. This value is a numerical parameter, that does not affect the result but has a strong effect on the convergence rate. Moulinec and Suquet (Moulinec and Suquet 1998) proposed the following Lamé constants of an isotropic reference medium as the best choice for fastest convergence:

$$\begin{equation}{\lambda }^{0}=\frac{1}{2}\left(\underset{\mathbf{x}\in {\Omega}}{\mathrm{inf}}\enspace \lambda \left(\mathbf{x}\right)+\underset{\mathbf{x}\in {\Omega}}{\mathrm{sup}}\enspace \lambda \left(\mathbf{x}\right)\right)\quad \text{and}\quad {\mu }^{0}=\frac{1}{2}\left(\underset{\mathbf{x}\in {\Omega}}{\mathrm{inf}}\enspace \mu \left(\mathbf{x}\right)+\underset{\mathbf{x}\in {\Omega}}{\mathrm{sup}}\enspace \mu \left(\mathbf{x}\right)\right),\end{equation} \tag{ 55 }$$

where $\lambda \left(\mathbf{x}\right)$ and $\mu \left(\mathbf{x}\right)$ are the elastic properties of the phases contained in the domain under consideration.

Algorithm 1.  Basic scheme

Data: ${\mathbb{C}}^{0}$, $\mathbb{C}(\mathbf{x})$, tol,$\bar{\boldsymbol{\varepsilon }}$

Result: ɛ (x)

${\boldsymbol{\varepsilon }}^{0}(\mathbf{x})=\bar{\boldsymbol{\varepsilon }}$

while $\frac{{\Vert}\nabla \cdot {\boldsymbol{\sigma }}^{i}{{\Vert}}_{L2}}{{\Vert}\left\langle {\mathbf{\Omega }\mathbf{\sigma }}^{i}\right\rangle {\Vert}}\hspace{-1.5pt} > \mathrm{t}\mathrm{o}\mathrm{l}$ do

${\boldsymbol{\tau }}^{i}(\mathbf{x})=\left[\mathbb{C}(\mathbf{x})-{\mathbb{C}}^{0}\right]:{\boldsymbol{\varepsilon }}^{i}(\mathbf{x})$

${\hat{\boldsymbol{\tau }}}^{i}(\boldsymbol{\xi })=\mathcal{F}({\boldsymbol{\tau }}^{i}(\mathbf{x}))$

${\hat{\tilde{\boldsymbol{\varepsilon }}}}^{i+1}(\boldsymbol{\xi })=-{\hat{\mathbb{\Gamma }}}^{0}(\boldsymbol{\xi }):{\hat{\boldsymbol{\tau }}}^{i}(\boldsymbol{\xi })$

${\boldsymbol{\varepsilon }}^{i+1}(\mathbf{x})={\mathcal{F}}^{-1}({\hat{\tilde{\boldsymbol{\varepsilon }}}}^{i+1}(\boldsymbol{\xi }))+\bar{\boldsymbol{\varepsilon }}$

end

The convergence (algorithm stopping) criterion originally proposed for the basic scheme reads as

$$\begin{equation}\frac{{\Vert}\nabla \cdot \boldsymbol{\sigma }{{\Vert}}_{{L}_{2}}}{{\Omega}{\Vert}< \hspace{-1.5pt}\boldsymbol{\sigma }\hspace{-1.5pt} > {\Vert}}=\frac{\sqrt{{\sum }_{\boldsymbol{\xi }}\vert \boldsymbol{\xi }\cdot \hat{\boldsymbol{\sigma }}{\vert }^{2}}}{{\Vert}\hat{\boldsymbol{\sigma }}(0){\Vert}}< \hspace{-1.5pt}\mathrm{t}\mathrm{o}\mathrm{l},\end{equation} \tag{ 56 }$$

where ${\Vert}\hspace{-1.5pt}{\bullet}\hspace{-1.5pt}{{\Vert}}_{{L}_{2}}$ represents the L₂-norm of the vector field and ||•|| is the Frobenius norm of a tensor. Invoking Parseval's theorem, the expression can be easily evaluated in Fourier space. Alternative convergence criteria can be found based on the compatibility condition or the imposed macroscopic conditions (Monchiet and Bonnet 2012, Moulinec and Silva 2014).

Only a few years after the development of the basic scheme, a family of accelerated solvers were developed (Eyre and Milton 1999, Michel et al 2000, 2001) with the aim to improve upon the convergence rate offered by the basic scheme. These methods (Eyre and Milton 1999, Michel et al 2000, 2001) and others proposed later (Monchiet and Bonnet 2012, Moulinec and Silva 2014, Schneider et al 2019) share the same approach, which is to introduce a dual strain and/or dual stress based formulation as well as the modification of the polarization term in the Lippmann–Schwinger equation to include strain compatibility and stress equilibrium with a prescribed macroscopic load.

The first accelerated solver was developed by Eyre and Milton (1999) to solve the periodic heat conduction problem. Michel et al (2001) developed an analogous scheme for the mechanical problem, called the accelerated scheme. Years later, Monchiet and Bonnet (2012) presented a more general class of algorithms, called polarization-based schemes, deriving motivation from the work of Eyre and Milton (1999) and Michel et al (2001). It was later shown by Moulinec and Silva (2014) that all these algorithms, although being developed from very different ideas, converge into a class of generalized schemes. In this section, the generalization proposed in Moulinec and Silva (2014) to derive all the methods from a single set of equations will be followed and this group of methods will be referred to as the polarization methods. However, the augmented Lagrangian scheme (Michel et al 2000, 2001), which is also part of this group, will be presented separately.

$\bar{\boldsymbol{\tau }}$. An implicit equation is derived by pre-multiplying the Lippmann–Schwinger equation by ${\left[{\mathbb{C}}^{0}\right]}^{-1}:\left[\mathbb{C}(\mathbf{x})-{\mathbb{C}}^{0}\right]$ and splitting it into two terms that are factorized with a coefficient α, which is a user-specified coefficient that affects the convergence rate. Rearranging terms and expressing the result as a function of the polarization field τ , which becomes the unknown of the problem, the implicit polarization scheme equation in its expanded form reads as

$$\begin{align}\hat{\boldsymbol{\tau }}\left(\boldsymbol{\xi }\right)& =\bar{\boldsymbol{\tau }}-\alpha {\mathbb{C}}^{0}:{\hat{\mathbb{\Gamma }}}^{0}\left(\boldsymbol{\xi }\right):\left({\mathbb{C}}^{0}:\mathcal{F}\left\{{\left[\mathbb{C}(\mathbf{x})-{\mathbb{C}}^{0}\right]}^{-1}:\boldsymbol{\tau }(\mathbf{x})\right\}+\hat{\boldsymbol{\tau }}(\boldsymbol{\xi })\right)\\ & \quad +\alpha \left[{\mathbb{C}}^{0}-{\mathbb{C}}^{0}:{\hat{\mathbb{\Gamma }}}^{0}\left(\boldsymbol{\xi }\right):{\mathbb{C}}^{0}\right]:\mathcal{F}\left\{{\left[\mathbb{C}(\mathbf{x})-{\mathbb{C}}^{0}\right]}^{-1}:\boldsymbol{\tau }(\mathbf{x})\right\}.\end{align} \tag{ 57 }$$

where α is a user-specified coefficient α = 2 and α = 1 yield the accelerated scheme (Eyre and Milton 1999, Michel et al 2000) and the augmented Lagrangian scheme (Michel et al 2000, 2001), respectively. A general expression of the polarization scheme that includes the definition of more coefficients affecting the convergence rate is provided in Monchiet and Bonnet (2012), Monchiet (2015).

As in the basic scheme, the discretization of the polarization implicit equation is performed using the DFT approach. The algorithm for the polarization scheme is shown in algorithm 2, where equation (57) is solved iteratively via Picard's iterations (fixed point iterations); a reference medium needs to be defined and it will affect the convergence rate. The optimal choice for the elastic constants of an isotropic reference medium proposed are

$$\begin{equation}{\lambda }^{0}=\sqrt{\underset{\mathbf{x}\in {\Omega}}{\mathrm{inf}}\enspace \lambda \left(\mathbf{x}\right)\cdot \underset{\mathbf{x}\in {\Omega}}{\mathrm{sup}}\enspace \lambda \left(\mathbf{x}\right)}\quad \text{and}\quad {\mu }^{0}=\sqrt{\underset{\mathbf{x}\in {\Omega}}{\mathrm{inf}}\enspace \mu \left(\mathbf{x}\right)\cdot \underset{\mathbf{x}\in {\Omega}}{\mathrm{sup}}\enspace \mu \left(\mathbf{x}\right)}.\end{equation} \tag{ 58 }$$

An exhaustive study of the choice of the numerical parameters of the method has been performed by To et al (2017).

Algorithm 2. Polarization scheme

Data: ${\mathbb{C}}^{0}$, $\mathbb{C}(\mathbf{x})$, tol,$\bar{\boldsymbol{\tau }}$

Result: ɛ (x)

${\boldsymbol{\tau }}^{0}(\mathbf{x})=\bar{\boldsymbol{\tau }}$

while $\frac{{\Vert}{\boldsymbol{\tau }}^{i+1}(\mathbf{x})-{\boldsymbol{\tau }}^{\boldsymbol{i}}(\mathbf{x}){{\Vert}}_{{L}_{2}}}{{\Vert}{\boldsymbol{\tau }}^{i+1}{{\Vert}}_{{L}_{2}}}\hspace{-1.5pt} > \mathrm{t}\mathrm{o}\mathrm{l}$ do

${\boldsymbol{\tau }}^{i}(\mathbf{x})={\mathcal{F}}^{-1}({\hat{\boldsymbol{\tau }}}^{i}(\boldsymbol{\xi }))$

${\boldsymbol{\varepsilon }}^{i}(\mathbf{x})={\left[\mathbb{C}(\mathbf{x})-{\mathbb{C}}^{0}\right]}^{-1}:{\boldsymbol{\tau }}^{i}(\mathbf{x})$

${\hat{\boldsymbol{\varepsilon }}}^{i}(\boldsymbol{\xi })=\mathcal{F}({\boldsymbol{\varepsilon }}^{i}(\mathbf{x}))$

${\hat{\boldsymbol{\sigma }}}^{i}(\boldsymbol{\xi })={\mathbb{C}}^{0}:{\hat{\boldsymbol{\varepsilon }}}^{i}(\boldsymbol{\xi })+{\hat{\boldsymbol{\tau }}}^{i}(\boldsymbol{\xi })$

${\hat{\boldsymbol{\tau }}}^{i+1}(\boldsymbol{\xi })={\hat{\boldsymbol{\tau }}}^{i}(\boldsymbol{\xi })-\alpha {\mathbb{C}}^{0}:{\hat{{\Gamma}}}^{0}(\boldsymbol{\xi }):{\hat{\boldsymbol{\sigma }}}^{i}(\boldsymbol{\xi })+\alpha \left[{\mathbb{C}}^{0}-{\mathbb{C}}^{0}:{\hat{{\Gamma}}}^{0}\left(\boldsymbol{\xi }\right):{\mathbb{C}}^{0}\right]:{\hat{\boldsymbol{\varepsilon }}}^{i}\left(\boldsymbol{\xi }\right)$

end

The following convergence criterion for the ith iteration was prescribed for algorithm 2 by Eyre and Milton (1999):

$$\begin{equation}\frac{{\Vert}{\boldsymbol{\tau }}^{i+1}(\mathbf{x})-{\boldsymbol{\tau }}^{i}(\mathbf{x}){{\Vert}}_{{L}_{2}}}{{\Vert}{\boldsymbol{\tau }}^{i+1}{{\Vert}}_{{L}_{2}}}< \hspace{-1.5pt}\mathrm{t}\mathrm{o}\mathrm{l}.\end{equation} \tag{ 59 }$$

Other convergence criteria have also been proposed (Moulinec and Silva 2014).

The augmented Lagrangian scheme

Michel et al (Michel et al 2000, 2001) originally developed the so-called augmented Lagrangian scheme to tackle problems with infinite phase contrast. This scheme is optimal for composites with very soft inclusions. In addition to the standard strain and stress fields, ɛ, σ , the method introduces two auxiliary fields $\mathbf{e}\left(\mathbf{x}\right)$, $\boldsymbol{\lambda }\left(\mathbf{x}\right)$, and forces their compatibility and equilibrium, respectively, by minimization with the augmented Lagrangian method. It was later shown by Moulinec and Silva (2014) that the algorithm becomes a particular choice of the polarization scheme with α = 1.

The augmented Lagrangian scheme involves solving a constrained optimization problem that reads as

$$\begin{equation}\underset{\mathbf{e}\left(\mathbf{x}\right)}{\mathrm{I}\mathrm{n}\mathrm{f}}\left\{\underset{\boldsymbol{\varepsilon }\left(\mathbf{x}\right)}{\mathrm{I}\mathrm{n}\mathrm{f}}\left\{{\int }_{{\Omega}}\frac{1}{2}\mathbf{e}\left(\mathbf{x}\right):\mathbb{C}\left(\mathbf{x}\right):\mathbf{e}\left(\mathbf{x}\right)\mathrm{d}{\Omega}\right\}\right\},\end{equation} \tag{ 60 }$$

where the constraints affecting the problem (compatibility condition) are prescribed as

$$\begin{equation}\boldsymbol{\varepsilon }\left(\mathbf{x}\right)-\mathbf{e}\left(\mathbf{x}\right)=0,\end{equation} \tag{ 61 }$$

with $\boldsymbol{\varepsilon }\left(\mathbf{x}\right)$ being periodic, compatible and with average $\bar{\boldsymbol{\varepsilon }}$ (imposed). The resulting optimization yields the augmented Lagrangian expression that reads as

$$\begin{equation} {\mathcal{L}}_{{\mathbb{C}}^{0}}\left(\boldsymbol{\varepsilon },\mathbf{e},\boldsymbol{\lambda }\right)={\int }_{{\Omega}}\left(\frac{1}{2}\mathbf{e}:\mathbb{C}(\mathbf{x}):\mathbf{e}+\boldsymbol{\lambda }:\left(\boldsymbol{\varepsilon }-\mathbf{e}\right)+\frac{1}{2}\left(\boldsymbol{\varepsilon }-\mathbf{e}\right):{\mathbb{C}}^{0}:\left(\boldsymbol{\varepsilon }-\mathbf{e}\right)\right)\mathrm{d}{\Omega},\end{equation} \tag{ 62 }$$

where $\boldsymbol{\lambda }\left(\mathbf{x}\right)$ denotes the Lagrange multiplier associated with the constraint in equation (61). The first term of equation (62) is the strain energy, the second one is the Lagrange multiplier term and the third one is the penalty term associated with the constraint (61).

The constrained optimization problem in equation (60) now becomes a saddle-point problem of ${\mathcal{L}}_{{\mathbb{C}}^{0}}$ for the field e as follows:

$$\begin{equation}\underset{\mathbf{e}}{\mathrm{I}\mathrm{n}\mathrm{f}}\enspace \underset{\boldsymbol{\varepsilon }}{\mathrm{I}\mathrm{n}\mathrm{f}}\enspace \underset{\boldsymbol{\lambda }}{\mathrm{Sup}}\enspace {\mathcal{L}}_{{\mathbb{C}}^{0}}\left(\boldsymbol{\varepsilon },\mathbf{e},\boldsymbol{\lambda }\right)=\underset{\boldsymbol{\lambda }}{\mathrm{Sup}}\enspace \underset{\mathbf{e}}{\mathrm{I}\mathrm{n}\mathrm{f}}\enspace \underset{\boldsymbol{\varepsilon }}{\mathrm{I}\mathrm{n}\mathrm{f}}\enspace {\mathcal{L}}_{{\mathbb{C}}^{0}}\left(\boldsymbol{\varepsilon },\mathbf{e},\boldsymbol{\lambda }\right),\end{equation} \tag{ 63 }$$

where the interchange between Sup and Inf can be justified when the energy function defined in equation (60) is convex and has sufficient growth at infinity.

The saddle-point can be reached by means of Uzawa's iterative algorithm (Uzawa and Arrow 1989, Glowinski and Le Tallec 1989). During Uzawa's iterations, for a given value of λ ⁱ⁻¹ and eⁱ⁻¹, an elastic problem has to be solved. It consists of finding the compatible strain field ɛ ⁱ in a homogeneous reference medium ${\mathbb{C}}^{0}$ by solving the following problem,

$$\begin{equation}\nabla \cdot {\mathbb{C}}^{0}:{\boldsymbol{\varepsilon }}^{i}=-\nabla \cdot \left({\boldsymbol{\lambda }}^{i-1}-{\mathbb{C}}^{0}:{\mathbf{e}}^{i-1}+{\mathbb{C}}^{0}:\bar{\boldsymbol{\varepsilon }}\right)=-\nabla \cdot \boldsymbol{\tau }(\mathbf{x}),\end{equation} \tag{ 64 }$$

where the right-hand side is the divergence of a polarization term τ (x) in which all the fields involved are known from previous iterations. The solution of this problem, ${\boldsymbol{\varepsilon }}^{i}\left(\mathbf{x}\right)$, can be explicitly computed using the periodic Green's operator (equations (53) and (54)) associated with the reference medium. Then, the auxiliary $\mathbf{e}\left(\mathbf{x}\right)$ field is obtained by solving the following linear equation at each point,

$$\begin{equation}{\mathbf{e}}^{i}={\boldsymbol{\varepsilon }}^{i}+{\left(\mathbb{C}+{\mathbb{C}}^{0}\right)}^{-1}:\left({\boldsymbol{\lambda }}^{i-1}-\mathbb{C}:{\boldsymbol{\varepsilon }}^{i}\right).\end{equation} \tag{ 65 }$$

Finally the Lagrange multipliers are updated as follows:

$$\begin{equation}{\boldsymbol{\lambda }}^{i}={\boldsymbol{\lambda }}^{i-1}+{\mathbb{C}}^{0}:\left({\boldsymbol{\varepsilon }}^{i}-{\mathbf{e}}^{i}\right).\end{equation} \tag{ 66 }$$

Once convergence has been reached, e is (almost) equal to ɛ and λ is (almost) equal to σ . The resulting algorithm is shown in algorithm 3. The optimal choice of the reference medium coincides with that for the polarization scheme in equation (58). Michel et al (2000) proposed the following convergence criteria involving the auxiliary fields

$$\begin{equation} \frac{{\Vert}{\boldsymbol{\varepsilon }}^{i+1}(x)-{\mathbf{e}}^{i+1}(x){{\Vert}}_{{L}_{2}}}{{\Omega}{\Vert}\bar{\boldsymbol{\varepsilon }}{\Vert}}\hspace{-1.5pt} > \mathrm{t}\mathrm{o}\mathrm{l}\quad \text{and}\quad \frac{{\Vert}\mathbb{C}(\mathbf{x}):{\boldsymbol{\varepsilon }}^{i+1}(x)-{\boldsymbol{\lambda }}^{i+1}(x){{\Vert}}_{{L}_{2}}}{{\Omega}{\Vert}{\mathbb{C}}^{0}:\bar{\boldsymbol{\varepsilon }}{\Vert}}\hspace{-1.5pt} > \mathrm{t}\mathrm{o}\mathrm{l}.\end{equation} \tag{ 67 }$$

Algorithm 3. Augmented Lagrangian scheme

Data: ${\mathbb{C}}^{0}$, $\mathbb{C}(\mathbf{x})$, tol,$\bar{\boldsymbol{\varepsilon }}$

Result: ɛ (x)

${\mathbf{e}}^{0}(x)=\bar{\boldsymbol{\varepsilon }};{\boldsymbol{\lambda }}^{0}(x)=0$

while $\frac{{\Vert}{\boldsymbol{\varepsilon }}^{i+1}(x)-{\mathbf{e}}^{i+1}(x){{\Vert}}_{L2}}{{\Omega}{\Vert}\bar{\boldsymbol{\varepsilon }}{\Vert}}\hspace{-1.5pt} > \mathrm{t}\mathrm{o}\mathrm{l}\;\;\mathbf{\text{and}}\;\;\frac{{\Vert}\mathbb{C}(\mathbf{x}):{\boldsymbol{\varepsilon }}^{i+1}(x)-{\boldsymbol{\lambda }}^{i+1}(x){{\Vert}}_{L2}}{{\Omega}{\Vert}{\mathbb{C}}^{0}:\bar{\boldsymbol{\varepsilon }}{\Vert}}\hspace{-1.5pt} > \mathrm{t}\mathrm{o}\mathrm{l}$ do

${\boldsymbol{\tau }}^{i}(\mathbf{x})={\boldsymbol{\lambda }}^{i}(x)-{\mathbb{C}}^{0}:{\mathbf{e}}^{i}(\mathbf{x})$

${\hat{\boldsymbol{\tau }}}^{i}(\boldsymbol{\xi })=\mathcal{F}({\boldsymbol{\tau }}^{i}(x))$

${\hat{\tilde{\boldsymbol{\varepsilon }}}}^{i+1}(\boldsymbol{\xi })=-{\hat{\mathbb{\Gamma }}}^{0}(\boldsymbol{\xi }):{\hat{\boldsymbol{\tau }}}^{i}(\boldsymbol{\xi })$

${\boldsymbol{\varepsilon }}^{i+1}(\mathbf{x})={\mathcal{F}}^{-1}({\hat{\tilde{\boldsymbol{\varepsilon }}}}^{i+1}(\boldsymbol{\xi }))+\bar{\boldsymbol{\varepsilon }}$

${\mathbf{e}}^{i+1}(\mathbf{x})={\boldsymbol{\varepsilon }}^{i+1}(x)+{\left(\mathbb{C}(\mathbf{x})+{\mathbb{C}}^{0}\right)}^{-1}:({\mathbb{C}}^{0}{\boldsymbol{\varepsilon }}^{i+1}(\mathbf{x})+{\boldsymbol{\lambda }}^{i}(\mathbf{x}))$

${\boldsymbol{\lambda }}^{i+1}(\mathbf{x})={\boldsymbol{\lambda }}^{i}(x)+{\mathbb{C}}^{0}:({\boldsymbol{\varepsilon }}^{i+1}-{\mathbf{e}}^{i+1})\left.\right)$

end

Moulinec and Silva (2014) have shown that obtaining an accurate solution may require additional convergence criteria on compatibility and equilibrium of the dual fields e and λ , respectively.

An alternate approach to the accelerated schemes for improved convergence rates over the basic scheme was proposed by Zeman et al (2010) and Vondřejc et al (2012) for the homogenization of a vector field in the context of electrostatics. In this approach, the integral Lippmann–Schwinger is discretized using the trigonometric collocation method in which trigonometric polynomials are used to interpolate the value of the fields in real space. After projecting the Lippmann–Schwinger equation to this discrete functional space, the resulting expression in linear elastic materials can be considered as a linear system of equations in which the unknown is the value of the strain at the grid points. This equation can be solved efficiently using a Krylov method such as the conjugate gradient method or the biconjugate gradient method (Zeman et al 2010). The linear equation obtained from the Lippmann–Schwinger equation reads as

$$\begin{equation}\boldsymbol{\varepsilon }(\mathbf{x})+{\mathcal{F}}^{-1}\left\{{\hat{\mathbb{\Gamma }}}^{0}\left(\boldsymbol{\xi }\right):\mathcal{F}\left\{\left[\mathbb{C}(\mathbf{x})-{\mathbb{C}}^{0}\right]:\boldsymbol{\varepsilon }(\mathbf{x})\right\}\right\}=\bar{\boldsymbol{\varepsilon }},\end{equation} \tag{ 68 }$$

where ɛ (x) is the discrete (nodal based) value of the strain ɛ . Note that the left-hand side of equation (68) corresponds to a linear operator applied to the field ɛ (x), which can be directly used in the Krylov solver. The algorithm for the Krylov-based scheme is shown in algorithm 4.

Algorithm 4. Krylov-based scheme

Data: ${\mathbb{C}}^{0}$, $\mathbb{C}(\mathbf{x})$, tol, $\bar{\boldsymbol{\varepsilon }}$

Result: ɛ (x)

ɛ (x) = 0

Solve $\mathcal{A}(\boldsymbol{\varepsilon }(\mathbf{x}))=\bar{\boldsymbol{\varepsilon }}$ by conjugate gradients with a tolerance of tol with$\mathcal{A}(\boldsymbol{\varepsilon }(\mathbf{x})){:=}\boldsymbol{\varepsilon }(\mathbf{x})+{\mathcal{F}}^{-1}\left\{{\hat{\mathbb{\Gamma }}}^{0}\left(\boldsymbol{\xi }\right):\mathcal{F}\left\{\left[\mathbb{C}(\mathbf{x})-{\mathbb{C}}^{0}\right]:\boldsymbol{\varepsilon }(\mathbf{x})\right\}\right\}$

As in the case of the basic scheme and the polarization approaches, it is demonstrated that the solution is insensitive to the choice of the auxiliary reference medium (Zeman et al 2010). In summary, this scheme shows a higher convergence rate compared to the basic scheme, but still presents some limitation in terms of stiffness contrast of the phases, since the matrix associated to the linear system of equations becomes ill-posed (with a large condition number) for very large contrast.

In parallel to the collocation approach proposed in Zeman et al (2010) and Vondřejc et al (2012), Brisard et al (Brisard and Dormieux 2010, 2012, Brisard and Legoll 2014) adopted the Hill–Mandel principle to define a direct, point-wise, discretization of the integral Lippmann–Schwinger equation where a subsequent truncation of underlying infinite Fourier series is required for the use of the DFT (see section 2.1.3 for more details). As a by-product, it provides an energetically consistent rule for the homogenization of a boundary value problem. The polarization field is taken as the unknown variable, which allows to solve infinite phase contrasts. However, the method requires pre-computing a Green's operator, ${\hat{\mathbb{\Gamma }}}^{0c}(\boldsymbol{\xi })$, consistent with their variational approach whose computational cost is very high and is therefore replaced by alternative filtering strategies. In a discretized domain the resulting linear equation for this method reads as

$$\begin{equation}{\left[\mathbb{C}(\mathbf{x})-{\mathbb{C}}^{0}\right]}^{-1}:\boldsymbol{\tau }(\mathbf{x})+{\mathcal{F}}^{-1}\left\{{\hat{\mathbb{\Gamma }}}^{0c}(\boldsymbol{\xi }):\hat{\boldsymbol{\tau }}(\boldsymbol{\xi })\right\}=\bar{\boldsymbol{\varepsilon }}.\end{equation} \tag{ 69 }$$

Alternative approaches of a consistent Green's operator in the Lippmann–Schwinger equation have also been proposed in Eloh et al (2019).

In a non-linear geometry framework in FFT based homogenization, the problem is generally considered using a total Lagrangian approach in which X represent a material point in the reference configuration. As in the linear case, the RVE is formed by different phases, however, each phase can have its own constitutive non-linear behavior, arranged in a reference periodic volume Ω₀ (see figure 4) where each point X ∈ Ω₀ belongs to one of those phases.

The objective of the homogenization problem is almost identical to the small strain setting. For a given macroscopic deformation gradient, find the fluctuation in the periodic deformation gradient and first Piola–Kirchhoff stress fields in the reference configuration; these fields have to fulfill compatibility and equilibrium conditions.

The governing equations of the non-linear geometry problem defined in a reference RVE Ω₀ in the absence of body forces are:

$$\begin{equation}\begin{cases}{\nabla }_{0}\cdot \mathbf{P}(\mathbf{X})=0\quad \\ \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\enspace {< \hspace{-1.5pt}\mathbf{F}(\mathbf{X})\hspace{-1.5pt}\hspace{-1.5pt} > }_{{{\Omega}}_{0}}=\bar{\mathbf{F}}\quad \\ \mathbf{F}(\mathbf{X})\enspace \mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{d}\mathrm{i}\mathrm{c}\quad \\ \tilde{\mathbf{P}}(\mathbf{X})\cdot \mathbf{N}\enspace \mathrm{a}\mathrm{n}\mathrm{t}\mathrm{i}-\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{d}\mathrm{i}\mathrm{c}\end{cases},\enspace ,\end{equation} \tag{ 88 }$$

where ∇₀⋅ is the divergence operator in the reference configuration, P(X) is the local first Piola–Kirchhoff stress tensor field and F(X) is the local deformation gradient. The latter field is decomposed in its reference volume averaged value in the domain $\bar{\mathbf{F}}$ and the fluctuating field $\tilde{\mathbf{F}}(\mathbf{X})$.

Irrespective of the use of a finite or a small deformation framework, material non-linearities occur when some phase in the domain is described by a non-linear constitutive behavior. In the non-linear setting, the constitutive behavior involves the first Piola–Kirchhoff stress as a function of the deformation gradient: P(F, X). In addition, the non-linear response can be deformation history and/or rate dependent. In the most general case, the constitutive laws for linear and non-linear geometry can be defined by a non linear relationship that may depend on the strain measure, deformation rate and/or internal variables, which are collectively represented using α .

$$\begin{equation}\boldsymbol{\sigma }=\boldsymbol{\sigma }\left(\boldsymbol{\varepsilon },\boldsymbol{\alpha },\dot{\boldsymbol{\varepsilon }},\dots \enspace \right)\enspace \;\text{or}\;\enspace \mathbf{P}=\mathbf{P}\left(\mathbf{F},\boldsymbol{\alpha },\mathbf{L},\dots \enspace \right).\end{equation} \tag{ 89 }$$

For simplicity but without the loss of generality, internal variables will be omitted in the algorithms presented below.

Some non-linear approaches rely on the use of tangent stiffness operators $\mathbb{C}\left(\mathbf{x}\right)$. In a small strain framework with discrete time increments, these operators are defined as the partial derivative of the stress increment with respect to the strain increment $\mathbb{C}\left(\mathbf{x}\right)=\frac{\partial \mathbf{\Delta }\mathbf{\sigma }}{\partial \mathbf{\Delta }\mathbf{\varepsilon }}$. In a finite strain setting described using a total Lagrangian approach, the material tangent required is $\mathbb{K}\left(\mathbf{x}\right)=\frac{\partial \mathbf{P}}{\partial \mathbf{F}}$.

For a non-linear simulation, the loading path is usually discretized into regular or adaptive load increments. When the time variable does not play a role on the constitutive behavior, the loading path is usually divided into pseudo-time increments, whereas when the time variable is considered, loading path is divided into real time increments. At each increment, which corresponds to a (pseudo-)time increment, the macroscopic prescribed strain/stress tensors are updated proportionally and a new problem is solved using the previous converged increment as starting point. In the case of non-linear simulations for time-independent and path-independent constitutive models, as in the case of hyperelastic solids, this discretization is only used to obtain intermediate solutions. Otherwise, the (pseudo-)time discretization is used to impose a particular loading rate and path.

For a time t, the set of variables which are stored are σ _t , ɛ _t and α _t in the case of small strains and P_t , F_t and α _t for finite strains. For a given time increment Δt with an associated increment of prescribed macroscopic strain ${\Delta}\bar{\boldsymbol{\varepsilon }}$ or ${\Delta}\bar{\mathbf{F}}$, then the problem becomes incremental and consists in solving σ _t+Δt , ɛ _t+Δt , α _t+Δt , or in the case of non linear geometry, P_t+Δt , F_t+Δt , α _t+Δt such that equilibrium is satisfied, equations (38) or (88), under the updated prescribed macroscopic strain ${\bar{\boldsymbol{\varepsilon }}}_{t+{\Delta}t}$ or ${\bar{\mathbf{F}}}_{t+{\Delta}t}$.

In the case of small strains and non-linear material models, the Lippmann–Schwinger equation is reformulated by the replacing the definition of the polarization as

$$\begin{equation}\boldsymbol{\tau }(\mathbf{x})=\boldsymbol{\sigma }(\mathbf{x})-{\mathbb{C}}^{0}:\boldsymbol{\varepsilon }(\mathbf{x}),\end{equation} \tag{ 90 }$$

where the local stresses are computed directly from the definition of non-linear constitutive laws. Following the original form of equation (50) and using equation (90), the final expression problem yields into a non-linear Lippmann–Schwinger equation, which in Fourier space reads as

$$\begin{equation}\hat{\boldsymbol{\varepsilon }}(\boldsymbol{\xi })=\bar{\boldsymbol{\varepsilon }}-{\hat{\mathbb{\Gamma }}}^{0}(\boldsymbol{\xi }):\left(\hat{\boldsymbol{\sigma }}(\boldsymbol{\xi })-{\mathbb{C}}^{0}:\hat{\boldsymbol{\varepsilon }}(\boldsymbol{\xi })\right).\end{equation} \tag{ 91 }$$

The extension of the Lippmann–Schwinger to finite strains using a pure Lagrangian framework, the most common approach in the literature, was proposed by Lahellec et al (2003). As stated in the introduction of this section, under this framework the equilibrium is expressed in the reference configuration using the first Piola–Kirchhoff stress while the deformation gradient becomes the unknown in the problem. The reformulation of the Lippmann–Schwinger equation in this case corresponds to

$$\begin{equation}\hat{\mathbf{F}}(\boldsymbol{\xi })=\bar{\mathbf{F}}-{\hat{\mathbb{\Gamma }}}^{0f}(\boldsymbol{\xi }):\left(\hat{\mathbf{P}}\left(\boldsymbol{\xi }\right)-{\mathbb{K}}^{0}:\hat{\mathbf{F}}(\boldsymbol{\xi })\right),\end{equation} \tag{ 92 }$$

where ${\hat{\mathbb{\Gamma }}}^{0f}$ is the fourth order gamma function in Fourier space, defined as the derivative of the Green's function of the reference medium and given by

$$\begin{equation}{\hat{\mathbb{\Gamma }}}^{0f}={\hat{{\Gamma}}}_{ijkl}^{0f}={\left.{\hat{{G}^{0}}}_{ik}{\xi }_{j}{\xi }_{l}\right\vert }_{(ik)(jh)},\end{equation} \tag{ 93 }$$

where $\hat{{\mathbb{G}}^{0}}$ stands for the Green's function in Fourier space (acoustic tensor, equation (52)) and the symbol (ik)(jh) denotes symmetrization with respect to the indices i, k and j, h only. Note that, contrary to the gamma operator in the small strain setting, ${\mathbb{\Gamma }}^{0f}$ only has major symmetries.

Basic scheme

The basic scheme proposed in Moulinec and Suquet (1994, 1998) was also used, from its original introduction, to solve small strain non-linear homogenization problems (equation (91)) using the fixed-point iterative procedure. This procedure consists in solving equation (91) using the result of the solution obtained from the previous iteration to compute the right-hand side of the equation. The resulting algorithm, for a single time increment, has been given in algorithm 7. A reference medium ${\mathbb{C}}^{0}$ is required as in the linear case and its optimum can be calculated from the equivalent isotropic elastic properties of the non-linear behaviors using equation (55).

Algorithm 7. Non-linear small strain basic scheme.

Data: ${\mathbb{C}}^{0}$, σ ( ɛ (x)), tol, $\bar{\boldsymbol{\varepsilon }}$

Result: ɛ (x)

${\boldsymbol{\varepsilon }}^{0}(\mathbf{x})=\bar{\boldsymbol{\varepsilon }}$

while $\frac{{\Vert}\nabla \cdot {\boldsymbol{\sigma }}^{i}{{\Vert}}_{L2}}{{\Vert}\left\langle {\boldsymbol{\sigma }}^{i}\right\rangle {\Vert}}\hspace{-1.5pt} > \mathrm{t}\mathrm{o}\mathrm{l}$ do

${\boldsymbol{\tau }}^{i}(\mathbf{x})=\boldsymbol{\sigma }({\boldsymbol{\varepsilon }}^{i}(\mathbf{x}))-{\mathbb{C}}^{0}:{\boldsymbol{\varepsilon }}^{i}(\mathbf{x})$

${\hat{\boldsymbol{\tau }}}^{i}(\boldsymbol{\xi })=\mathcal{F}({\boldsymbol{\tau }}^{i}(\mathbf{x}))$

${\hat{\tilde{\boldsymbol{\varepsilon }}}}^{i+1}(\boldsymbol{\xi })=-{\hat{{\Gamma}}}^{0}(\boldsymbol{\xi }):{\hat{\boldsymbol{\tau }}}^{i}(\boldsymbol{\xi })$

${\boldsymbol{\varepsilon }}^{i+1}(\mathbf{x})={\mathcal{F}}^{-1}({\hat{\tilde{\boldsymbol{\varepsilon }}}}^{i+1}(\boldsymbol{\xi }))+\bar{\boldsymbol{\varepsilon }}$

end

This direct extension of the basic scheme to non-linear material problems adds non-linearity to the original implicit equation and does not make use of material tangents. As a consequence, the resulting convergence rate becomes even poorer than its linear counterpart.

In the case of the non-linear basic scheme for finite strains, Lahellec et al (2003) proposed a linearization of equation (92) to yield into a Newton–Raphson procedure, which makes use of the material tangent operators to improve the convergence rate. The linearized Lippmann–Schwinger equation is solved using the basic scheme at each Newton–Raphson iteration using the stress field and consistent tangents obtained in the previous iteration.

The deformation gradient field is also obtained in an iterative manner and its value for iteration i + 1 is given by Fⁱ⁺¹ = Fⁱ + δ F. The stress field is linearized around the deformation gradient by

$$\begin{equation}{\mathbf{P}}^{i+1}={\mathbf{P}}^{i}+{\left.\frac{\partial \mathbf{P}}{\partial \mathbf{F}}\right\vert }_{\mathbf{F}={\mathbf{F}}^{i}}:\delta \mathbf{F}={\mathbf{P}}^{i}+{\mathbb{K}}^{i}:\delta \mathbf{F},\end{equation} \tag{ 94 }$$

where ${\mathbb{K}}^{i}$ is the (non-symmetric) material tangent computed for the ith iteration. The fixed point form of the Lippmann–Schwinger equation is simplified then to

$$\begin{equation}{\hat{\mathbf{F}}}^{i+1}(\boldsymbol{\xi })={\hat{\mathbf{F}}}^{i}(\boldsymbol{\xi })-{\hat{\mathbb{\Gamma }}}^{0f}(\boldsymbol{\xi }):\left[\mathcal{F}\left\{{\mathbb{K}}^{i}(\mathbf{X}):\delta \mathbf{F}(\mathbf{X})\right\}+{\hat{\mathbf{P}}}^{i}(\mathbf{F}(\mathbf{X}))\right].\end{equation} \tag{ 95 }$$

Note that the term ${\hat{\mathbb{\Gamma }}}^{0f}(\boldsymbol{\xi }):{\mathbb{C}}^{0}:{\hat{\mathbf{F}}}^{i}={\hat{\tilde{\mathbf{F}}}}^{i}$ has been introduced in ${\hat{\bar{\mathbf{F}}}}^{i}$ to yield into ${\hat{\mathbf{F}}}^{i}$. The computational cost of this linearized version of the basic scheme is smaller, since the number of evaluations of the constitutive equations is much lower with respect to fixed point procedure, as a trade off providing the material consistent tangent.

These non-linear extensions for the basic scheme can be easily introduced into a (pseudo-)time incremental approach by computing the equilibrium state at the end of each discretized time increment. The resulting algorithm for finite strains is described in algorithm 8.

Algorithm 8. Non linear time incremental basic scheme for finite strains.

Data: ${\Delta}\bar{\mathbf{F}},{\Delta}t,{t}_{\text{final}},{\mathrm{t}\mathrm{o}\mathrm{l}}_{\text{lin}},{\mathrm{t}\mathrm{o}\mathrm{l}}_{nw}$

Result: Ft (X)

t = 0

Ft = I

while t ⩽ tfinal do

${\mathbf{F}}^{0}={\mathbf{F}}_{t}+{\Delta}{\bar{\mathbf{F}}}_{t+{\Delta}t}$

while ${\Vert}\boldsymbol{\xi }\cdot {\hat{\mathbf{P}}}^{i}(\boldsymbol{\xi }){{\Vert}}_{\infty }/{\Vert}{\hat{\mathbf{P}}}^{i}(\mathbf{0}){\Vert}\geqslant {\mathrm{t}\mathrm{o}\mathrm{l}}_{nw}$ do

Pi = P(Fi ); ${\mathbb{K}}^{i}={\left.\frac{\partial \mathbf{P}}{\partial \mathbf{F}}\right\vert }_{\mathbf{F}={\mathbf{F}}^{i}}$

Solve δ F by fixed point iteration method (algorithm 7) with a tolerance of tollin for: ${\hat{\delta \mathbf{F}}}^{i+1}(\boldsymbol{\xi })=-{\hat{\mathbb{\Gamma }}}^{0f}(\boldsymbol{\xi }):\left[\mathcal{F}\left\{{\mathbb{K}}^{i}(\mathbf{X}):\delta \mathbf{F}(\mathbf{X})\right\}+{\hat{\mathbf{P}}}^{i}(\mathbf{F}(\mathbf{X}))\right]$

Fi+1 = Fi + δ F

end

Ft+Δt = Fi+1

t = t + Δt

end

Polarization-based schemes

In Michel et al (1999, 2000, 2001), the accelerated and augmented Lagrangian schemes, belonging to the group of polarization schemes, were directly proposed for non-linear constitutive equations. The polarization based schemes were extended by Vinogradov and Milton (2008) for a non-linear geometry setting for thermoelastic problems; a Newton–Raphson procedure was also used in this extension. The polarization scheme is used at each Newton iteration to solve the linearized problem, using the stress and the local consistent tangent fields computed during the evaluation of the constitutive equations at the previous Newton iteration. In the small strain framework, the implicit expression to be solved is

$$\begin{equation}{\hat{\boldsymbol{\tau }}}^{i+1}(\boldsymbol{\xi })={\hat{\boldsymbol{\tau }}}^{i}-\alpha {\mathbb{C}}^{0}:{\hat{\mathbb{\Gamma }}}^{0}(\boldsymbol{\xi }):{\hat{\boldsymbol{\sigma }}}^{i}(\boldsymbol{\xi })+\alpha \left[{\mathbb{C}}^{0}-{\mathbb{C}}^{0}:{\hat{\mathbb{\Gamma }}}^{0}\left(\boldsymbol{\xi }\right):{\mathbb{C}}^{0}\right]:{\hat{\boldsymbol{\varepsilon }}}^{i}(\boldsymbol{\xi }),\end{equation} \tag{ 96 }$$

where the material consistent tangent is used inside the definition of polarization tensor at previous iteration ${\hat{\boldsymbol{\tau }}}^{i}$. The finite strains framework of the Newton–Raphson method combined with the polarization scheme was presented in Kabel et al (2014) by using the definitions described in equation (94). Regarding this method, the optimal choice of algorithmic parameters and alternatives to reduce the number of evaluations of the constitutive equations during the iterative procedure were studied in detail in Schneider et al (2019). Analogously to the basic scheme, the polarization scheme can be easily posed in a (pseudo-)time incremental algorithm similar to algorithm 8.

Krylov-based schemes

Gélébart and Mondon-Cancel (2013) first used the Krylov-based methods to solve non-linear problems. The method proposed is based on a Newton–Raphson algorithm and can be applied to various kinds of behaviors (time dependant or independent, with or without internal variables) through a conventional integration procedure as used in finite element codes. Similar to the basic scheme version (Lahellec et al 2003), the Newton–Raphson procedure consists in applying corrections to the unknown field ɛ until the convergence is reached.

$$\begin{equation}{\boldsymbol{\varepsilon }}^{i+1}={\boldsymbol{\varepsilon }}^{i}+\delta \boldsymbol{\varepsilon }\quad \text{and}\quad {\boldsymbol{\sigma }}^{i+1}={\boldsymbol{\sigma }}^{i}+{\left.\frac{\partial \boldsymbol{\sigma }}{\partial \boldsymbol{\varepsilon }}\right\vert }_{\boldsymbol{\varepsilon }={\boldsymbol{\varepsilon }}^{i}}:\delta \boldsymbol{\varepsilon }={\boldsymbol{\sigma }}^{i}+{\mathbb{C}}^{i}:\delta \boldsymbol{\varepsilon }.\end{equation} \tag{ 97 }$$

The strain correction field, δ ɛ , is obtained using a Krylov iterative algorithm (section 3.2.3). At the beginning of each time increment, the increment of the macroscopic strain is added to the (pseudo-)time discretization, $\bar{\boldsymbol{\varepsilon }}$ to the state of the strain field ɛ ⁱ by defining ${\boldsymbol{\varepsilon }}^{0}={\boldsymbol{\varepsilon }}_{t}+{\Delta}\bar{\boldsymbol{\varepsilon }}$. Then, the linear equation corresponds to

$$\begin{equation} \delta \boldsymbol{\varepsilon }+{\mathcal{F}}^{-1}\left\{{\hat{\mathbb{\Gamma }}}^{0}\left(\boldsymbol{\xi }\right):\mathcal{F}\left\{\left[{\mathbb{C}}^{i}(\mathbf{x})-{\mathbb{C}}^{0}\right]:\delta \boldsymbol{\varepsilon }(\mathbf{x})\right\}\right\}=-{\boldsymbol{\varepsilon }}^{i}+{\mathcal{F}}^{-1}\left\{{\hat{\mathbb{\Gamma }}}^{0}\left(\boldsymbol{\xi }\right):{\hat{\boldsymbol{\tau }}}^{i}(\mathbf{x})\right\}+{\bar{\boldsymbol{\varepsilon }}}^{i}.\end{equation} \tag{ 98 }$$

The small strain non-linear Krylov approach was extended to a non-linear geometrical setting by Kabel et al (2014) using the finite strains counterparts of stresses and strains defined in equation (94). In addition, the (pseudo-)time incremental approach for Krylov-based schemes consists on replacing the fixed point iterative solver in algorithm 8 by a Kryov solver.

Non-linear quasi-Newton approaches and non-linear conjugate gradient methods

Some schemes for the non-linear Lippmann–Schwinger equation, alternative to the standard Newton–Raphson, have been proposed recently by Schneider. In Schneider (2017), a Nesterov's method was proposed that leads to a significant speed up compared to the basic scheme for linear problems with moderate contrast, and compares favorably to the (Newton-)conjugate gradient method for some problems. In Schneider (2019), the Barzilai–Borwein implementation of the basic scheme was proposed, which is a fast and robust alternative algorithm that avoids the manual choice of algorithmic parameters.

An exhaustive comparison of different quasi-Newton methods for FFT homogenization has been presented in Wicht et al (2020). For example, Shanthraj et al (2015), Chen et al (2019b), Wicht et al (2021a) investigate the non-linear extension of Lippmann–Schwinger methods by Anderson acceleration and demonstrates that this combination leads to robust and fast general-purpose solvers.

Yet another alternative to Newton–Raphson to solve the non-linear Lippmann–Schwinger equation was explored in Schneider (2020a) who used non-linear conjugate gradient, Fletcher–Reeves algorithm, as a dynamical system with state-dependent nonlinear damping. It was demonstrated by numerical experiments that this approach may represent a competitive, memory-efficient and parameter-choice free solution method.

The non-linear extension of Fourier–Galerkin approaches was proposed in Zeman et al (2017), de Geus et al (2017). As in the linear case (section 3.3), the weak formulation of the equilibrium is derived using the virtual work principle. The problem statement, in its finite strain form, consists in finding, for a given macroscopic deformation gradient history $\bar{\mathbf{F}}(t)$, the deformation gradient field F(x) for every time t that fulfills

$$\begin{equation}\mathbb{G}\ast \mathbf{P}={\mathcal{F}}^{-1}\left\{\hat{\mathbb{G}}:\mathcal{F}\left\{\mathbf{P}\right\}\right\}=0,\end{equation} \tag{ 99 }$$

where the projection operator $\mathbb{G}$ (defined in equation (74)) is used similarly to the small strain linear version to project any arbitrary tensor field into its compatible part.

In the case of non-linear material behavior, the discrete expression of the weak form of equilibrium given by equation (99) defines a non-linear system of algebraic equations. This equation is discretized in time and then, for each time increment, is solved iteratively by the Newton–Raphson method using a linearization. Combining the equilibrium equation (99) with the linearization of deformation gradient and stress equation (94) the next equation is derived for the iteration i and time t with a target deformation gradient ${\bar{\mathbf{F}}}^{t}$.

$$\begin{equation}{\mathcal{F}}^{-1}\left\{\hat{\mathbb{G}}:\mathcal{F}\left\{{\mathbb{K}}^{i}:\delta \mathbf{F}\right\}\right\}=-{\mathcal{F}}^{-1}\left\{\hat{\mathbb{G}}:\mathcal{F}\left\{\mathbf{P}({\mathbf{F}}^{i})\right\}\right\}.\end{equation} \tag{ 100 }$$

The expression in equation (100) is a linear system of equations in which the unknown is the correction at the iteration i of the deformation gradient δ F. In that equation, Fⁱ is obtained in a previous iteration being ${\mathbf{F}}^{0}={\bar{\mathbf{F}}}^{t}$. This linear system, as it happened to its linear version, is rank deficient and the use of Krylov solvers is also fundamental here. The iterative Newton process stops when the equilibrium is achieved. This happens when either the right-hand side of the equation $\mathcal{G}(\mathbf{P})$ or the correction of the Newton–Raphson δ F become sufficiently small. The solving sequence is summarized in algorithm 9.

Algorithm 9. Non linear time incremental Fourier–Galerkin scheme for finite strains.

Data: ${\Delta}\bar{\mathbf{F}},{\Delta}t,{t}_{\text{final}},{\mathrm{t}\mathrm{o}\mathrm{l}}_{\text{lin}},{\mathrm{t}\mathrm{o}\mathrm{l}}_{nw}$

Result: Ft (X)

t = 0

Ft = I

while t ⩽ tfinal do

${\mathbf{F}}^{0}={\mathbf{F}}_{t}+{\Delta}{\bar{\mathbf{F}}}_{t+{\Delta}t}$

while ${\Vert}\delta \mathbf{F}{\Vert}/{\Vert}{\bar{\mathbf{F}}}^{i}{\Vert}\geqslant {\mathrm{t}\mathrm{o}\mathrm{l}}_{nw}$ do

Pi = P(Fi ); ${\mathbb{K}}^{i}={\left.\frac{\partial \mathbf{P}}{\partial \mathbf{F}}\right\vert }_{\mathbf{F}={\mathbf{F}}^{i}}$

Solve δ F by conjugate gradient with a tolerance of tollin for:${\mathcal{F}}^{-1}\left\{\hat{\mathbb{G}}:\mathcal{F}\left\{{\mathbb{K}}^{i}:\delta \mathbf{F}\right\}\right\}=-{\mathcal{F}}^{-1}\left\{\hat{\mathbb{G}}:\mathcal{F}\left\{{\mathbf{P}}^{i}\right\}\right\}$

Fi+1 = Fi + δ F

end

Ft+Δt = Fi+1

t = t + Δt

end

The extension of FFT methods based on the displacement field to finite strains and/or non-linear behavior was proposed in Schneider et al (2017), Lucarini and Segurado (2019b). In the case of finite strains, the equilibrium is posed in the reference configuration. As in the linear case, the variable in which the problem is defined is the displacement field, in particular the displacement fluctuations $\tilde{\mathbf{u}}$. The non-linear differential equation defined in the equilibrium equation is then solved iteratively by the Newton–Raphson method.

For the DBFFT approach developed by Lucarini and Segurado (2019b), the displacement fluctuation field and the first Piola stress are linearized around the last iteration i of the displacement field fluctuations ${\tilde{\mathbf{u}}}^{i}$ leading to

$$\begin{equation}{\tilde{\mathbf{u}}}^{i+1}={\tilde{\mathbf{u}}}^{i}+\delta \tilde{\mathbf{u}}\end{equation} \tag{ 101 }$$

$$\begin{equation}\mathbf{P}(\mathbf{I}+{\nabla }_{0}({\mathbf{u}}^{k+1}\left(\mathbf{x}\right)))\approx \mathbf{P}({\bar{\mathbf{F}}}_{\bar{U}}^{k+1}+{\nabla }_{0}{\tilde{\mathbf{u}}}^{i})+{\mathbb{K}}^{i}:{\nabla }_{0}\delta \tilde{\mathbf{u}}\enspace \;\text{with}\;\enspace {\mathbb{K}}^{i}=\frac{\partial \mathbf{P}}{\partial {\nabla }_{0}\delta \tilde{\mathbf{u}}}. \end{equation} \tag{ 102 }$$

Combining the equilibrium equation equation (88) with the linearization of the stress equation (102), the following equation is obtained

$$\begin{equation}{\nabla }_{0}\cdot {\mathbb{K}}^{i}:{\nabla }_{0}\delta \tilde{\mathbf{u}}=-{\nabla }_{0}\cdot \mathbf{P}\left(\mathbf{I}+{\nabla }_{0}{\mathbf{u}}^{i}\right)\end{equation} \tag{ 103 }$$

which is a linear differential equation whose unknown is $\delta \tilde{\mathbf{u}}$. In equation (103), the initialization of the variables are adapted from the averaged incremental boundary conditions ${\nabla }_{0}{\mathbf{u}}^{0}={\bar{\mathbf{F}}}_{t+{\Delta}t}-\mathbf{I}+{\nabla }_{0}{\tilde{\mathbf{u}}}_{t}$ and ${\mathbb{K}}^{i}$ denotes the material consistent tangent evaluated in the ith iteration.

Then, following the same approach as in the linear case (see section 3.4), the linear differential equation is transformed to Fourier space and the fields are replaced by their discrete counterparts. The result is a linear system of complex numbers

$$\begin{equation}\hat{\mathbb{d}}:\mathcal{F}\left({\mathbb{K}}^{i}:\left({\mathcal{F}}^{-1}\left(\hat{\mathbb{g}}\cdot \delta \hat{\tilde{\mathbf{u}}}\right)\right)\right)=-\hat{\mathbb{d}}:\mathcal{F}\left(\mathbf{P}\left({\nabla }_{0}{\mathbf{u}}^{i}+\mathbf{I}\right)\right)\end{equation} \tag{ 104 }$$

in which the unknown vector to solve is $\delta \tilde{\mathbf{u}}$. The linear equation for each Newton iteration (104) is formally identical to the linear equation obtained for the case of linear elasticity (equation (86)). The difference is that in equation (104) the gradient operator $\hat{\mathbb{g}}$ is no longer symmetric being this operator defined in Fourier space as

$$\begin{equation}{\left[\hat{\mathbb{g}}\left(\boldsymbol{\xi }\right)\right]}_{ijk}={\delta }_{ik}{\xi }_{j}.\end{equation} \tag{ 105 }$$

The resulting system of equations in complex numbers becomes fully determined removing the real Fourier transform symmetries and the zero frequency. As in the other cases, the linear equation can be solved by any iterative solver. The Newton iterations finish when the norm of correction ${\nabla }_{0}\delta \tilde{\mathbf{u}}$ at the last linear iteration is below a tolerance tol_nw .

In the case of a finite strain non-linear analysis, the definition of the preconditioner M is very similar to the linear version, being replaced in equation (87) the average stiffness matrix by the average material tangent $\bar{\mathbb{K}}$. This preconditioner is computed at the beginning of every Newton iteration. The resulting algorithm is detailed below (algorithm 10).

Algorithm 10. Non linear time incremental DBFFT method for finite strains.

Data: ${\Delta}\bar{\mathbf{F}},{\Delta}t,{t}_{\text{final}},{\mathrm{t}\mathrm{o}\mathrm{l}}_{\text{lin}},{\mathrm{t}\mathrm{o}\mathrm{l}}_{nw}$

Result: ut (X)

t = 0

ut = 0

while t ⩽ tfinal do

${\mathbf{u}}^{0}={\mathbf{u}}_{t}+{\Delta}\bar{\mathbf{F}}\cdot \mathbf{X}$

while ${\Vert}{\nabla }_{0}\delta \tilde{\mathbf{u}}{\Vert}/{\Vert}\bar{{\nabla }_{0}{\mathbf{u}}^{i}}+\mathbf{I}{\Vert}\geqslant {\mathrm{t}\mathrm{o}\mathrm{l}}_{nw}$ do

Pi = P(∇0 ui + I); ${\mathbb{K}}^{i}={\left.\frac{\partial \mathbf{P}}{\partial \mathbf{F}}\right\vert }_{\mathbf{F}={\nabla }_{0}{\mathbf{u}}^{i}+\mathbf{I}}$

Solve $\delta \hat{\tilde{\mathbf{u}}}$ by conjugate gradient with a tolerance of tollin for:$\mathbf{M}\cdot \hat{\mathbb{d}}:\mathcal{F}\left\{{\mathbb{K}}^{i}:\left({\mathcal{F}}^{-1}\left\{\hat{\mathbb{g}}\cdot \delta \hat{\tilde{\mathbf{u}}}\right)\right\}\right\}=-\mathbf{M}\cdot \hat{\mathbb{d}}:\mathcal{F}\left\{{\mathbf{P}}^{i}\right\}$

${\mathbf{u}}^{i+1}={\mathbf{u}}^{i}+\delta \tilde{\mathbf{u}}$

end

ut+Δt = ui+1

t = t + Δt

end

In these classical schemes, load is introduced imposing the whole macroscopic strain tensor. When the stress tensor is the macroscopic data, the dual problem in a straightforward manner to overcome it. In this case, the input data for the algorithms are the fourth order compliance tensors and the fluctuating stress tensor becomes unknown of the problem as stated in Bhattacharya and Suquet (2005), Monchiet and Bonnet (2012), Monchiet (2015).

When mixed control is required (e.g. uniaxial tension), iterative approaches are the usual solution (Michel et al 1999). Under this procedure, the algorithm still imposes the macroscopic strain tensor as input of the algorithm, but this tensor is corrected during the iterative process in order to obtain, at the end of the iteration process, the macroscopic stress tensor components imposed as input.

In Michel et al (1999) it is presented an algorithm for small strains that includes corrections as a residual formulation in each iteration into the prescribed macroscopic strain via

$$\begin{equation} {\bar{\varepsilon }}_{IJ}^{i+1}={\bar{\varepsilon }}_{IJ}^{i}+{\left[{< \hspace{-1.5pt}\mathbb{C}\hspace{-1.5pt}\hspace{-1.5pt} > }^{-1}:\left(\frac{{\bar{\boldsymbol{\varepsilon }}}^{i}:\bar{\boldsymbol{\sigma }}+\bar{\boldsymbol{\sigma }}:\left({< \hspace{-1.5pt}\mathbb{C}\hspace{-1.5pt} > }^{-1}:< \hspace{-1.5pt}{\boldsymbol{\sigma }}^{i}\hspace{-1.5pt} > -{\bar{\boldsymbol{\varepsilon }}}^{i}\right)}{\bar{\boldsymbol{\sigma }}:{< \hspace{-1.5pt}\mathbb{C}\hspace{-1.5pt} > }^{-1}:\bar{\boldsymbol{\sigma }}}\bar{\boldsymbol{\sigma }}-< \hspace{-1.5pt}{\boldsymbol{\sigma }}^{i}\hspace{-1.5pt} > \right)\right]}_{IJ},\end{equation} \tag{ 106 }$$

where $\left\langle \enspace \right\rangle$ denotes the spatial average of the field inside and the convergence is reached when ${\bar{\varepsilon }}_{IJ}^{i+1}={\bar{\varepsilon }}_{IJ}^{i}$. An alternative simplified version was proposed in Lebensohn and Cazacu (2012), Eisenlohr et al (2013) where, in the case of finite strains, the imposed macroscopic deformation gradient is corrected as

$$\begin{equation}{\bar{F}}_{IJ}^{i+1}={\bar{F}}_{IJ}^{i}-{\left[{< \hspace{-1.5pt}{\mathbb{K}}^{i}\hspace{-1.5pt} > }^{-1}:\left(< \hspace{-1.5pt}{\mathbf{P}}^{i}\hspace{-1.5pt} > -\bar{\mathbf{P}}\right)\right]}_{IJ},\end{equation} \tag{ 107 }$$

where the terms of the averaged fourth order tensor $\mathbb{K}$ where the stress is not imposed are removed, then inverted and refilled.

As a result of the corrections to the applied strain during iterations, the number of iterations per step significantly increases with respect to strain control versions.

An alternative approach for classic schemes under mixed control was proposed by Kabel et al (2016) under the finite strains formalism. The method consists in modifying the Krylov-based Lippmann–Schwinger algorithm through the introduction of projection tensors to account for the load control and in solving the problem using an iterative Krylov solver. This method for stress control is more efficient than the iterative approaches. The fourth order projector tensors, $\mathbb{Q}$ and $\mathbb{M}={\left[\mathbb{Q}:{\mathbb{C}}^{0}:\mathbb{Q}\right]}^{-1}$ act on the null frequency and $\mathbb{Q}$ accounts for the components and directions where the stresses are imposed. The Lippmann–Schwinger equation in its implicit form reads as

$$\begin{equation} \mathbf{F}=\bar{\mathbf{F}}+\mathbb{M}:\left(\bar{\mathbf{P}}-\mathbb{Q}:{\mathbb{C}}^{0}:\bar{\mathbf{F}}\right)-\mathbb{M}:\mathbb{Q}:{\left\langle \mathbf{P}\left(\mathbf{F}\right)-{\mathbb{C}}^{0}:\mathbf{F}\right\rangle }_{{{\Omega}}_{0}}-{{\Gamma}}^{0}:\left(\mathbf{P}\left(\mathbf{F}\right)-{\mathbb{C}}^{0}:\mathbf{F}\right).\end{equation} \tag{ 108 }$$

A similar idea was exploited by Lucarini and Segurado (2019a) for the Fourier–Galerkin method but, instead of modifying the equation, changes are performed only on the Fourier projection operator, preserving the simple structure of the resulting equations. This algorithm modifies the null frequency of the projector operator, replacing the components where the stress is imposed by the fourth order identity tensor. This allows the introduction of an average imposed stress (in the form of the macroscopic first Piola stress tensor) into the equilibrium equation reading as

$$\begin{equation}{\mathcal{G}}^{\ast }\left(\mathbf{P}\left(\mathbf{F}\right)-\bar{\mathbf{P}}\right)=0,\end{equation} \tag{ 109 }$$

where $\bar{\mathbf{P}}$ is the target average stress tensor. Note that the only components of this tensor that will have an effect the solution will be the components IJ which are stress controlled and which are the ones in which the projector operator has been modified. It is demonstrated that the method does not imply any extra computational cost compares with the strain controlled version.

$$\begin{equation}{\hat{{G}^{\ast }}}_{ijkl}=\begin{cases}{\delta }_{ik}{\delta }_{jl} & \text{if}\enspace \mathbf{\xi }=0\;\;\text{for}\;\;\text{stress}\;\;\text{controled}\;\;IJ\;\;\text{terms} \\ {0}_{ijkl} & \text{if}\enspace \mathbf{\xi }=0\;\;\text{for}\;\;\text{strain}\;\;\text{controlled}\;\;ij\;\;\text{terms} \\ {0}_{ijkl} & \text{for}\;\;\text{Nyquist}\;\;\text{frequencies} \\ {\delta }_{ik}\frac{{\xi }_{j}{\xi }_{l}}{\boldsymbol{\xi }\cdot \boldsymbol{\xi }}\quad & \text{for}\enspace \mathbf{\xi }\ne 0 \end{cases}.\end{equation} \tag{ 110 }$$

In the case of the displacement based method (Lucarini and Segurado 2019b) the deformation gradient field will be split into the components due to macroscopic imposed Piola stress ${\bar{\mathbf{F}}}_{\bar{f}}$ and the ones due to prescribed deformation gradient ${\bar{\mathbf{F}}}_{\bar{U}}$.

$$\begin{equation}\mathbf{F}\left(\mathbf{x}\right)=\mathbf{I}+{\left[{\bar{\mathbf{F}}}_{\bar{U}}-\mathbf{I}\right]}_{ij}+{\left[{\bar{\mathbf{F}}}_{\bar{f}}-\mathbf{I}\right]}_{IJ}+{\nabla }_{0}\tilde{\mathbf{u}},\end{equation} \tag{ 111 }$$

where ${\bar{\mathbf{F}}}_{\bar{f}}$ is unknown and is only defined for the IJ components where Piola stress is imposed. Solving the equilibrium equation (104) at each time increment by an iterative Newton approach, the linear system to solve at each Newton iteration i is

$$\begin{align}\begin{aligned}\hat{\mathbb{d}}:\mathcal{F}\left({\mathbb{K}}^{i}:\left({\mathcal{F}}^{-1}\left(\hat{\mathbb{g}}\cdot \delta \hat{\tilde{\mathbf{u}}}\right)+\delta {\bar{\mathbf{F}}}_{\bar{f}}\right)\right)& =-\hat{\mathbb{d}}:\mathcal{F}\left(\mathbf{P}\left({\nabla }_{0}{\mathbf{u}}^{i}+\mathbf{I}\right)\right)\\ {\left[\mathcal{F}\left({\mathbb{K}}^{i}:\left({\mathcal{F}}^{-1}\left(\hat{\mathbb{g}}\cdot \delta \hat{\tilde{\mathbf{u}}}\right)+\delta {\bar{\mathbf{F}}}_{\bar{f}}\right)\right)\left(0\right)\right]}_{IJ}& ={\left[{\bar{\mathbf{P}}}_{t+{\Delta}t}\right]}_{IJ}-{\left[\mathcal{F}\left(\mathbf{P}\left({\nabla }_{0}{\mathbf{u}}^{i}+\mathbf{I}\right)\right)\left(0\right)\right]}_{IJ}\end{aligned},\end{align} \tag{ 112 }$$

being ${\nabla }_{0}{\mathbf{u}}^{i}={\left[{\bar{\mathbf{F}}}_{\bar{U}\enspace t+{\Delta}t}-\mathbf{I}\right]}_{ij}+{\left[{\bar{\mathbf{F}}}_{\bar{f}}^{i-1}-\mathbf{I}\right]}_{IJ}+{\nabla }_{0}{\tilde{\mathbf{u}}}^{i-1}$ and starting the algorithm with ${\nabla }_{0}{\tilde{\mathbf{u}}}^{0}={\nabla }_{0}{\tilde{\mathbf{u}}}_{t}$ and ${\bar{\mathbf{F}}}_{\bar{f}}^{0}={\bar{\mathbf{F}}}_{\bar{f}\enspace t}$. In (112), the unknown to be solved is the variable $\left\{\hat{\tilde{\mathbf{u}}}\vert {\bar{\mathbf{F}}}_{\bar{f}}\right\}$ and, if the zero frequency and doubled terms due to the symmetries of the real Fourier transform are removed in the first equation, the system becomes fully determined and Hermitian.

The first application of the FFT approach to study polycrystals was proposed by Lebensohn (2001), who developed the viscoplastic FFT (VP-FFT), a fixed-point iterative algorithm to compute the local response of elastic and pure viscoplastic anisotropic 3D polycrystals under the action of a macroscopic velocity gradient. The VP-FFT formulation is closely related with the viscoplastic self consistent scheme (VPSC, Lebensohn and Tomé (1993)) using a pure visco-plastic small strain formulation accounting for the crystal rotations and using the same tangent linearization. The local behavior at a material point x inside a grain is represented with a crystal viscoplastic model in which ${\dot{\boldsymbol{\varepsilon }}}^{p}$ is the plastic strain rate tensor, which is a function of the deviatoric Cauchy stress σ ' as

$$\begin{equation}{\dot{\boldsymbol{\varepsilon }}}^{p}(\mathbf{x})=\sum\limits _{s}\text{sym}\left({\mathbf{m}}^{s}(\mathbf{x})\right){\dot{\gamma }}^{s}=\sum\limits _{s}\text{sym}\left({\mathbf{m}}^{s}(\mathbf{x})\right){\dot{\gamma }}_{0}^{s}{\left(\frac{\mathbf{m}(\mathbf{x}):{\boldsymbol{\sigma }}^{\prime }(\mathbf{x})}{{\tau }_{\mathrm{c}}^{s}(\mathbf{x})}\right)}^{n},\end{equation} \tag{ 114 }$$

where m^s is the Schmid tensor for slip system s, ${\dot{\gamma }}^{s}$ is the slip rate, ${\tau }_{\mathrm{c}}^{s}$ is the critical resolved shear stress for slip system s, ${\dot{\gamma }}_{0}^{s}$ is the reference shear rate and n is the power law exponent, which is the inverse of the rate sensitivity parameter. It is important to remark that this original model (Lebensohn 2001) does not consider elastic strains and, therefore, the total strain is isochoric. The viscoplastic response of each material point is incorporated using a tangent approximation,

$$\begin{equation}{\boldsymbol{\sigma }}^{\prime }(\mathbf{x})={\mathbb{M}}^{tg-1}(\mathbf{x}):\mathbf{d}(\mathbf{x})+{\mathbf{S}}^{0}(\mathbf{x}),\end{equation} \tag{ 115 }$$

where d(x) is the strain rate, S⁰(x) is the back-extrapolated stress and ${\mathbb{M}}^{tg}(\mathbf{x})$ is the tangent compliance modulus, defined by deriving equation (114) with respect to the deviatoric Cauchy stress, and given by

$$\begin{equation}{\mathbb{M}}^{tg}(\mathbf{x})=\frac{\partial {\dot{\boldsymbol{\varepsilon }}}^{p}(\mathbf{x})}{\partial {\boldsymbol{\sigma }}^{\prime }(\mathbf{x})}=n\sum\limits _{s}{\dot{\gamma }}_{0}^{s}\frac{{\boldsymbol{m}}^{s}(\mathbf{x}):{\boldsymbol{m}}^{s}(\mathbf{x})}{{\tau }_{\mathrm{c}}^{s}(\mathbf{x})}{\left(\frac{{\boldsymbol{m}}^{s}(\mathbf{x}):{\boldsymbol{\sigma }}^{\prime }(\mathbf{x})}{{\tau }_{\mathrm{c}}^{s}(\mathbf{x})}\right)}^{n-1}.\end{equation} \tag{ 116 }$$

Following the basic scheme, a homogeneous reference medium is defined. In this case the medium follows a pure visco-plastic linear behavior which relates macroscopic stress and strain rates (Σ' and D) using a tangent approach as

$$\begin{equation}{\mathbf{\Sigma }}^{\prime }={\mathbb{L}}^{0}:\mathbf{D}+{\mathbf{S}}^{\infty },\end{equation} \tag{ 117 }$$

where ${\mathbb{L}}^{0}$ is the macroscopic tangent stiffness modulus and S^∞ is the back-extrapolated stress of the reference medium. The local perturbation field in the deviatoric stress, φ ' is given by

$$\begin{equation}{\boldsymbol{\varphi }}^{\prime }(\mathbf{x})={\tilde{\boldsymbol{\sigma }}}^{\prime }(\mathbf{x})+{\mathbb{L}}^{0}:\tilde{\mathbf{d}}(\mathbf{x}),\end{equation} \tag{ 118 }$$

where $\tilde{\mathbf{d}}(\mathbf{x})=\mathbf{d}(\mathbf{x})-\mathbf{D}$ and ${\tilde{\boldsymbol{\sigma }}}^{\prime }(\mathbf{x})={\boldsymbol{\sigma }}^{\prime }(\mathbf{x})-{\mathbf{\Sigma }}^{\prime }$ are the local fluctuations in the strain rate and deviatoric stress. Assuming incompressibility, the local system of PDEs for the VP-FFT problem is derived from the stress equilibrium as

$$\begin{equation}\begin{aligned}\nabla \cdot ({\mathbb{L}}^{0}:{\nabla }^{s}\mathbf{v}(\mathbf{x}))+{\nabla }^{s}\cdot {\boldsymbol{\varphi }}^{\prime }(\mathbf{x})-\nabla p(\mathbf{x})& =0,\quad \forall \enspace \mathbf{x}\in V,\\ \nabla \cdot \mathbf{v}& =0,\quad \forall \enspace \mathbf{x}\in V,\end{aligned}\end{equation} \tag{ 119 }$$

where v is the velocity field and p is the pressure.

Using the error in equilibrium vs number of iterations as a benchmark, Lebensohn performed a series of convergence tests of the VP-FFT model in a polycrystalline RVE by varying spatial resolution in 3D, rate sensitivity, reference medium stiffness (Voigt and Reuss averages, and tangent modulus obtained from the VPSC estimate (Lebensohn and Tomé 1993)). Two crystal phases were included in the simulations, hard and soft phases, having the RVEs different volume fraction of the hard phase. It was shown that lowest error in equilibrium is obtained for the finer grids and the convergence is achieved faster for lower rate sensitivity, Voigt reference medium, single phase polycrystals, and polycrystals with lower volume fraction of the hard phase. Lebensohn then studied the predictive capabilities of the VP-FFT model and compared the texture evolution due to rolling of copper up to 50% thickness reduction predicted by the VPSC model, which models each grain as a viscoplastic homogeneous inclusion with homogeneous values of the stress and strain rates, with respect the VP-FFT model in which regions within a grain can evolve differently from the rest of the grain (see figure 13). The VP-FFT predicted texture was found to be the closest to the experimental texture.

Zoom In Zoom Out Reset image size

Since its conception, the VP-FFT model has been applied to many different problems including but not limited to comparisons with homogenization estimates for the effective behavior and statistical fluctuations of stress and strain-rate fields in 2D and 3D viscoplastic polycrystals (Lebensohn et al 2004, 2005, 2008). The VP-FFT model allows taking direct input from orientation imaging microscopy images of polycrystals, which is a meaningful endeavor because the VP-FFT model can account for the effect of grain topology and grain neighborhood on the local and macroscopic response.

Zoom In Zoom Out Reset image size

An elastic-viscoplastic FFT (EVP-FFT) model was first proposed by Lebensohn et al (Lebensohn et al 2012) within a small strain formulation using the augmented Lagrangian scheme for accelerated convergence. In this model, both hydrostatic and deviatoric parts of the Cauchy stress tensor were taken into account such that

$$\begin{equation}\boldsymbol{\sigma }(\mathbf{x})=\mathbb{C}(\mathbf{x}):\left(\boldsymbol{\varepsilon }(\mathbf{x})-{\boldsymbol{\varepsilon }}^{p}(\mathbf{x})\right),\end{equation} \tag{ 120 }$$

where $\mathbb{C}(\mathbf{x})$ is the local elastic stiffness, ɛ ^p (x) is the plastic strain whose evolution is given by equation (114) and ɛ (x) is the local total strain, which is the symmetric part of the displacement gradient u_k,l (x). Then the FFT basic scheme (section 3.2.1) is applied combined with an iterative technique for mixed macroscopic control (see section 3.6.1), based on the method proposed in Michel et al (2001).

Lebensohn et al (2012) applied the EVP-FFT model to study the interplay between elastic and plastic anisotropy, and its effects on the macroscopic and local response during elastic–plastic transition by comparing the response of the Cu polycrystal with that of an artificial polycrystalline material whose single crystal elastic properties are chosen such that it gives the same macroscopic elastic response as the polycrystal Cu. Furthermore, the same plastic properties were used for both materials. The comparison showed the superiority of the EVP-FFT model over the VP-FFT model predictions that had been previously performed by Rollett et al (2010) for a similar configuration (figure 14).

(Eisenlohr et al (2013) generalized the small strain EVP-FFT model (Lebensohn et al 2012) to a finite strain formulation to get the first finite strain EVP-FFT model. The finite strain approach followed here uses a total Lagrangian approach in which deformation gradient is used as kinematic variable similar to the adaptation of the basic scheme for finite strains described in section 3.5.3.

In Eisenlohr et al (2013) the finite strain EVP-FFT predictions were compared with those from finite strain EVP-FEM simulations showing that both methods matched very well at macroscopic and microscopic level (in figure 15). In addition, it is demonstrated that the finite strain EVP-FFT could simulate microstructures that the FEM code could not handle for such fine resolution. The FFT method in this application showed a small but noticeable high-frequency fluctuation (Gibbs oscillations). It is reported that the strain/stress fields spatial variability in the FFT simulation resulted significantly larger than that of the corresponding FEM at a given mesh/grid resolution stating that this observation is a consequence of the solutions and not the spurious oscillations.

Zoom In Zoom Out Reset image size

The method proposed in Eisenlohr et al (2013) was further improved by Shanthraj et al (2015) from the numerical view point. In this work, the basic-scheme is replaced by different non-linear efficient approaches including the augmented Lagrangian scheme, the polarization scheme and several non-linear solvers for systems of equations as non-linear GMREs. The result was a very efficient approach for large simulations of polycrystals. A more recent alternative FFT implementation of finite strain crystal plasticity was proposed in Lucarini and Segurado (2019c). The model followed the same finite strain description as Eisenlohr et al (2013), Shanthraj et al (2015) but relied on the Fourier–Galerkin approach for solving the microfields. The main difference with respect to other approaches was the absence of reference medium, which alleviates the need of recomputing the average tangent response at each iteration. A subsequent improvement of the CP-FFT implementation in Lucarini and Segurado (2019c) was the development of an efficient non-iterative control scheme (Lucarini and Segurado 2019a) for general mixed boundary conditions.

Lebensohn and Needleman (Lebensohn and Needleman 2016) extended the small strain EVP-FFT model to propose the first FFT-based numerical implementation of the non-local (strain gradient) plasticity theory of Gurtin (2000, 2002). In addition to the advantages of the local EVP-FFT model, the non-local EVP-FFT model allows accounting for slip system level back-stresses generated by the presence and accumulation of dislocations during plasticity. However, in addition to the standard stress and strain-rate boundary conditions in the conventional EVP-FFT model, one also needs to impose periodicity of plastic shear γ^s at the slip-system level, provide micro-stress ( ξ ^s —work-conjugate to ∇γ^s ) conditions at internal boundaries and ensure that the average microstress is consistent with the state prescribed. The governing PDEs of the non-local EVP-FFT model are

$$\begin{equation}\begin{aligned}& \nabla \cdot \boldsymbol{\sigma }(\mathbf{x})=0,\quad \forall \enspace \mathbf{x}\in V\\ & {\pi }^{s}(x)-{\tau }^{s}(x)-\nabla \cdot {\boldsymbol{\xi }}^{s}(\mathbf{x})=0,\quad \forall \;s\in {N}^{s};\enspace \forall \enspace \mathbf{x}\in V\end{aligned},\end{equation} \tag{ 121 }$$

where π^s is the work conjugate to γ^s and τ^s (x) = m^s : σ is the resolved shear stress. In the non-local EVP-FFT approach, the following constitutive and kinematic relationships are respected at each material point x ∈ V:

$$\begin{equation}\begin{aligned}\boldsymbol{\sigma }& =\mathbb{C}:(\boldsymbol{\varepsilon }-{\boldsymbol{\varepsilon }}^{p})\\ {\dot{\boldsymbol{\varepsilon }}}^{p}& =\sum\limits _{s}{\mathbf{m}}^{s}{\dot{\gamma }}_{0}^{s}{\left(\frac{{\pi }^{s}}{{\tau }_{\mathrm{c}}^{s}}\right)}^{n}\mathrm{sign}\left({\pi }^{s}\right)=\sum\limits _{s}{\mathbf{m}}^{s}{\dot{\gamma }}_{0}^{s}{\left(\frac{{\tau }^{s}+\nabla \cdot {\boldsymbol{\xi }}^{s}}{{\tau }_{\mathrm{c}}^{s}}\right)}^{n}\mathrm{sign}\left({\tau }^{s}+\nabla \cdot {\boldsymbol{\xi }}^{s}\right)\\ {\boldsymbol{\xi }}^{s}& =-{l}^{2}{\pi }_{0}\mathbf{X}:\left[\left(\nabla \times {\boldsymbol{\beta }}^{p}\right)\cdot {\mathbf{m}}^{s}\right]\end{aligned},\end{equation} \tag{ 122 }$$

where β ^p is the plastic distortion (a second order tensor), whose symmetric part is the plastic strain, l and π₀ are material parameters that have the dimensions of length and stress, and X is the third order Levi-Civita permutation tensor with components in Einstein summation notation as e_ijk ; the tensorial operations used in the equation ξ ^s are adopted from (Upadhyay et al 2013).

The parameter l can be tuned with respect to the grain size d in order to increase or decrease the non-local effect on the local and macroscopic mechanical response, as can be seen in figure 16.

Zoom In Zoom Out Reset image size

A crucial step in the non-local EVP-FFT model is to compute ∇ ⋅ ξ ^s , which involves the partial derivatives of the spatial distribution of the Schmidt tensor, m^s (x), and plastic distortion β ^p . It is well known that partial derivatives are easily computed in Fourier space, however, using standard DFT based on the continuum Fourier transform can result in the formation of Gibbs oscillations (as reviewed in section 2.1.3). A workaround to this problem is to compute the partial derivatives using some finite difference scheme, which corresponds to the use of some modified frequencies in Fourier space Neumann et al (2001), Müller (2002), Press et al (2002). In this work the central difference scheme is used, resulting in a strong reduction of the spurious oscillations in the derivatives.

More recently, Marano et al (2021) also proposed an FFT implementation of the non-local strain gradient plasticity model of Gurtin (2000, 2002) in a manner similar to Lebensohn and Needleman (2016). Figure 17 shows the formation of kink-bands under a single slip condition as predicted by the local EVP-FFT model and the strain gradient FFT model of Marano et al (2021). The strain gradient FFT predicted bands were found to be larger in number and intensity in comparison to the local EVP-FFT predictions.

Zoom In Zoom Out Reset image size

Upadhyay and co-workers (Upadhyay 2014, Upadhyay et al 2016a) extended the classical continuum based EVP-FFT model to a couple stress continuum based EVP-FFT (CSEVP-FFT) approach in order to capture the role of grain boundaries on the local and bulk mechanical response of nanocrystalline materials. The model takes into account the role of local changes in the elastic and plastic orientations induced during deformation and their role on the Cauchy stress response via a strong coupling between the symmetric Cauchy stress tensor σ and the second-order deviatoric couple stress tensor M^D through the higher order equilibrium equation (Upadhyay et al 2013, Upadhyay 2014, Upadhyay et al 2016a):

$$\begin{equation}\nabla \cdot \boldsymbol{\sigma }-\frac{1}{2}\nabla \cdot \left[\mathbf{X}\cdot \left(\nabla \cdot {\mathbf{M}}^{D}\right)\right]=0.\end{equation} \tag{ 123 }$$

The couple stress tensor has the units N m⁻¹ and it is the work-conjugate of the second-order elastic curvature tensor κ ^e , which is the gradient of the elastic rotation vector ω ^e such that ${\mathbf{M}}^{D}=\mathbb{A}:{\boldsymbol{\kappa }}^{e}$, where $\mathbb{A}$ is a 4th order tensor with units N (Upadhyay et al 2013). Development of the CSEVP-FFT model has led to the formulation of an extended periodic Lippmann–Schwinger type equation arising from the higher order equilibrium equation (123). The modified Green's function for this extended periodic Lippmann–Schwinger type equation is obtained as a solution to the following equation:

$$\begin{equation} \nabla \cdot \left({\mathbb{C}}^{0}:\nabla \mathbf{G}(\boldsymbol{x}-{\boldsymbol{x}}^{\prime })\right)+\frac{1}{4}\nabla \cdot \left[\mathbf{X}\cdot \left(\nabla \cdot \left({\mathbb{A}}^{0}:\nabla \left(\nabla \times \mathbf{G}(\boldsymbol{x}-{\boldsymbol{x}}^{\prime })\right)\right)\right)\right]+1\delta (\boldsymbol{x}-{\boldsymbol{x}}^{\prime })=0,\end{equation} \tag{ 124 }$$

where ${\mathbb{C}}^{0}$ and ${\mathbb{A}}^{0}$ are elastic moduli of a reference medium and $\mathbf{G}\left({\boldsymbol{x}-\boldsymbol{x}}^{\prime }\right)$ is the Green's tensor, which for the CSEVP problem is expressed in Fourier space using index notation as

$$\begin{equation}{\hat{G}}_{ki}(\xi )={\left({\xi }_{l}{\xi }_{j}{C}_{ijkl}^{0}-\frac{1}{4}{\xi }_{o}{\xi }_{n}{\xi }_{l}{\xi }_{j}{e}_{mok}{e}_{ijp\hspace{.5pt}}{A}_{plmn}^{0}\right)}^{-1}.\end{equation} \tag{ 125 }$$

The solution to (124) gives the compatible total strain and curvature tensors as

$$\begin{equation}\begin{aligned}\boldsymbol{\varepsilon }(\boldsymbol{x})& =\boldsymbol{E}+\frac{1}{2}{\mathcal{F}}^{-1}\left[i\left(\hat{\mathbf{G}}(\boldsymbol{\xi })\cdot \hat{\boldsymbol{f}}(\boldsymbol{\xi })+\hat{\boldsymbol{f}}(\boldsymbol{\xi })\cdot {\hat{\mathbf{G}}}^{\mathrm{T}}(\boldsymbol{\xi })\right)\otimes \boldsymbol{\xi }\right]\\ \boldsymbol{\kappa }(\boldsymbol{x})& =\boldsymbol{K}-\frac{1}{2}{\mathcal{F}}^{-1}\left[\left\{\mathbf{X}:\left[\boldsymbol{\xi }\otimes \left(\hat{\mathbf{G}}(\boldsymbol{\xi })\cdot \hat{\boldsymbol{f}}(\boldsymbol{\xi })\right)\right]\right\}\otimes \boldsymbol{\xi }\right]\end{aligned},\end{equation} \tag{ 126 }$$

where $i=\sqrt{-1}$ and $\boldsymbol{f}(\boldsymbol{x})=\nabla \cdot \boldsymbol{\tau }(\boldsymbol{x})+\frac{1}{2}\nabla \cdot \left[\mathbf{X}\cdot \left(\nabla \cdot \boldsymbol{\mu }(\boldsymbol{x})\right)\right]$.

In order to avoid Gibbs' oscillations due to higher order derivatives involved in the computation of equation set (126), these equations are solved using the finite difference based approach to compute the frequency multipliers in Fourier space corresponding to partial derivatives (Neumann et al 2001, Müller 2002, Press et al 2002); henceforth, this approach is known as the discrete differentiation approach.

The CSEVP-FFT model (Upadhyay 2014, Upadhyay et al 2016a) requires imposing macroscopic couple stress and curvature rate conditions in addition to the conventional Cauchy stress and strain rate. These additional conditions allow simulating bending and torsion of the RVEs. The implementation of the couple stress continuum model takes into account the size dependent response. Grain boundary interfaces are characterized using elastic curvatures that are representative of their structure and defect content. Depending on the magnitude and distribution of these curvatures, local Cauchy stresses are generated in the grain boundary neighborhood. These stresses result in the activation of slip systems besides those predicted from the Schmid law. At the macro scale, this results in a strain rate dependent 'softening' or the inverse Hall–Petch effect (figure 18), which is typical for nanocrystalline materials with grain sizes less than a few tens of nanometers.

Zoom In Zoom Out Reset image size

The CSEVP-FFT model Upadhyay et al (2016a) also predicts the formation and evolution of geometrically necessary dislocations via the temporal evolution of the polar dislocation density tensor $\boldsymbol{\alpha }=-\nabla \times {\boldsymbol{\varepsilon }}^{p}+{\boldsymbol{\kappa }}^{p\mathrm{T}}-\text{tr}({\boldsymbol{\kappa }}^{p})1$. In addition, to dislocations, the model also predicts the formation of disclinations, which are rotational type of line crystal defects, through the evolution of their polar density θ = −∇ × κ ^p . The evolutions of ${\Vert}\boldsymbol{\alpha }{\Vert}$ and ${\Vert}\boldsymbol{\theta }{\Vert}$ were studied and they were found to be the highest in the vicinity of grain boundaries.

Haouala et al (2020) proposed a so-called lower-order strain gradient plasticity FFT model that is based on a physically based crystal plasticity model in finite strains which enriches the viscoplastic constitutive relationship (114) by accounting for the contribution from both statistically stored and geometrically necessary dislocation densities, ρ_SSD and ρ_GND, respectively, on the critical resolved shear stress ${\tau }_{\mathrm{c}}^{s}$. The evolution of ρ_SSD is given by the balance between the rate of accumulation and annihilation of dislocations according to the widely used Kocks–Mecking law Mecking and Kocks (1981). The evolution of ρ_GND was obtained by solving a constrained minimization problem based on the L₂ minimization approach proposed in the work of Arsenlis and Parks (1999). The constitutive equations for ${\tau }_{\mathrm{c}}^{s}$ of this model (Haouala et al 2020) are:

$$\begin{equation}\begin{aligned}{\tau }_{\mathrm{c}}^{s}& =\mu {b}^{s}\sqrt{\sum\limits _{{s}^{\prime }}{h}^{s{s}^{\prime }}{\rho }^{{s}^{\prime }}}\\ {\dot{\rho }}_{\text{SSD}}^{s}& =\frac{1}{{b}^{s}}\left(\frac{1}{{l}^{s}}-2{y}_{\mathrm{c}}{\rho }_{\text{SSD}}\right)\left\vert {\dot{\gamma }}^{s}\right\vert \\ {l}^{s}& =\frac{K}{\sqrt{{\sum }_{{s}^{\prime }\ne s}\left({\rho }_{\text{SSD}}^{{s}^{\prime }}+{\rho }_{\text{GND}}^{{s}^{\prime }}\right)}}\\ {\bar{\rho }}_{\text{GND}}& ={\text{arg min}}_{\bar{\rho }}\left\{{\bar{\rho }}^{\mathrm{T}}\bar{\rho }+\lambda \left(\mathbf{A}\bar{\rho }-\bar{\alpha }\right)\right\}={\left({\mathbf{A}}^{\mathrm{T}}\mathbf{A}\right)}^{-1}{\mathbf{A}}^{\mathrm{T}}\bar{\alpha },\end{aligned}\end{equation} \tag{ 127 }$$

where μ is the shear modulus, b^s is the magnitude of the Burgers vector of dislocations on a slip system s, y_c is the effective dislocation annihilating distance, l^s is the dislocation mean-free path, K is the similitude coefficient Sauzay and Kubin (2011), ${\bar{\rho }}_{\text{GND}}$ is a vector of dimension N_s (the total number of slip systems) containing the GND density on each system s, $\bar{\alpha }$ is the vector representation of the polar dislocation density tensor, A is a matric of dimensions 9 × N_s , and λ are Lagrange multipliers. From the algorithmic view point, the model relies on using a Fourier–Galerkin approach de Geus et al (2017) for the mechanical solver. The GNDs from equation (127) are computed at load step every increment using the curl operator and, similar to all the previous higher order crystal plasticty models reviewed here, a discrete differentiation rule is used to reduce noise.

The lower order model (Haouala et al 2020) does not involve introducing any ad hoc length scale parameter being the non-local effect controlled by the length of the burgers vector. Therefore, the model present a simple structure as that of the local EVP-FFT model and yet it is able to capture the size-dependent plasticity response as can be seen in figure 19.

Zoom In Zoom Out Reset image size

FFT approaches present several benefits in the estimation of the fatigue response of polycrystalline materials since the very good performance for these methods allow to solve very detailed microstructures or to perform many simulations for a statistical analysis. The first use of FFT approaches for fatigue was presented by Rovinelli et al (Rovinelli et al 2015) who used the small-strain based EVP-FFT model (Lebensohn et al 2012) to analyze the microstructure variability in short crack growth behavior. The authors analyze different fatigue indicator parameters (FIPs) as possible driving forces for short crack growth and based on their values postulated potential crack paths. This work was further extended to analyze real experiments performed in a synchrotron using machine learning techniques (Rovinelli et al 2018a, 2018b).

An alternative way was followed by Lucarini and Segurado (Lucarini and Segurado 2019c, 2020) who propose the use of FFT homogenization to estimate fatigue life as function of cyclic loading and polycrystalline microstructure, through a microstructure sensitive fatigue life prediction approach (McDowell and Dunne 2010). The model by Lucarini and Segurado (2019c) is a finite strain crystal elasto-visco plastic approach based on the Fourier Galerkin formulation (de Geus et al 2017, Zeman et al 2017). In order to account for kinematic hardening and softening that could occur during cyclic loading, a modified version of the viscoplastic power law (114) was used to model the crystal behavior

$$\begin{equation}{\dot{\gamma }}^{s}={\dot{\gamma }}_{0}{\left(\frac{\vert {\tau }^{s}-{\chi }^{s}\vert }{{g}_{\mathrm{c}}^{s}+{g}_{s}^{s}}\right)}^{\frac{1}{m}}\enspace \mathrm{sign}({\tau }^{s}-{\chi }^{s}).\end{equation} \tag{ 136 }$$

In equation (136), the critical resolved shear stress on the s slip system is the denominator and has two contributions, a monotonic term ${g}_{\mathrm{c}}^{\alpha }$ and a cyclic softening term ${g}_{s}^{s}$. The flow rule also includes a dependency of the slip rate with a back stress, χ^s , that determines the kinematic hardening contribution and which evolution follows the modified version of the Ohno and Wang model proposed in (Cruzado et al 2017).

Based on the microscopic fields obtained in the simulation of a cyclic deformation of a polycrystalline RVE, the approach computes a series of fatigue indicator parameters to estimate the fatigue life of the polycrystal. In particular two FIPs are considered: (i) the grain averaged accumulated plastic slip per cycle ${P}_{\text{cyc}}^{g}$ and (ii) the crystallographic plastic strain energy ${W}_{\text{cyc}}^{g}$. The two FIPs ${P}_{\text{cyc}}^{g,\mathrm{max}}\enspace$ and ${W}_{\text{cyc}}^{b,\mathrm{max}}\enspace$ are respectively associated with two integration regions, (i) the volume occupied by a grain and (ii) the volume defined, for each integration point, by the bands of a given width contained in the grain and parallel to the slip planes (Castelluccio and McDowell 2015) such that,

$$\begin{equation}\begin{aligned}{P}_{\text{cyc}}^{g,\mathrm{max}}\enspace & =\underset{g=1,{N}^{g}}{\mathrm{max}}\frac{1}{{V}^{g}}\underset{{V}^{g}}{\int }{P}_{\text{cyc}}(x)\mathrm{d}{V}^{g}=\underset{g=1,{N}^{g}}{\mathrm{max}}\frac{1}{{V}^{g}}\underset{{V}^{g}}{\int }\underset{\text{cyc}}{\int }{\mathbf{L}}_{\mathbf{p}}(\mathbf{x}):{\mathbf{L}}_{\mathbf{p}}(\mathbf{x})\mathrm{d}t\enspace \mathrm{d}{V}^{g}\\\ {W}_{\text{cyc}}^{b,\mathrm{max}}\enspace & =\underset{b=1,{N}^{b}}{\mathrm{max}}\left\{\underset{s}{\mathrm{max}}\frac{1}{{V}^{b}}\underset{{V}^{b}}{\int }{W}_{\text{cyc}}^{s}(\mathbf{x})\mathrm{d}{V}^{b}\right\}\\ & =\underset{b=1,{N}^{b}}{\mathrm{max}}\left\{\underset{s}{\mathrm{max}}\frac{1}{{V}^{b}}\underset{{V}^{b}}{\int }\underset{\text{cyc}}{\int }{\tau }^{s}{\dot{\gamma }}^{s}(\mathbf{x})\mathrm{d}t\enspace \mathrm{d}{V}^{b}\right\},\end{aligned}\end{equation} \tag{ 137 }$$

where P_cyc(x) is the local value of the accumulated plastic slip per cycle, ${W}_{\text{cyc}}^{s}(x)$ is the local value of the energy density dissipated per cycle for each slip system s, ${\mathbf{L}}_{\mathbf{p}}(x)={\sum }_{s}{\;\mathbf{m}}^{s}{\dot{\gamma }}^{s}$ is the plastic velocity gradient, V^g is the volume of the gth grain and N^g is the total number of grains, V^b is the volume of the bth band and N^b is the total number of bands in the RVE.

Two FFT based implementations of the finite strain FIP model were proposed by Lucarini and Segurado (Lucarini and Segurado 2019c): (i) classical CFT based FFT and (ii) discrete differentiation approach implementation to correct the Gibbs phenomena by introducing new projection operators based on discrete derivatives (Willot et al 2014). Both implementations were validated by comparing their predictions with FEM predicted macroscopic response, microscopic fields and probability to nucleate cracks in an fcc RVE under different cyclic loading histories. The material simulated was Inconel 718. For the FFT simulations, the polycrystalline RVEs are subjected to a nearly uniaxial-stress cyclic strain history. Figure 26(a) shows the comparison between P_cyc(x) predicted from the FEM, FFT and discrete differentiation based DFT models for a slice of this RVE after three cycles. Qualitatively, the FIP patterns obtained with the three numerical implementations are very close. Next, the histograms of ${P}_{\text{cyc}}^{g,\mathrm{max}}\enspace$ and ${W}_{\text{cyc}}^{b,\mathrm{max}}$ are plotted for the three models as shown in figures 26(b) and (c), respectively. The discrete operator result was more near than the standard FFT approach to the FEM results. In order to further highlight the different between the standard and discrete FFT implementations, a distribution of stiff spherical particles (ten times stiffer than the Inconel 718 crystals) were introduced into the RVE and their effect on the cyclic response was studied. Figure 26 shows a slice of the RVE with stiffer second phase particles and the line plot of the stress component σ_xx vs the position along the white line shown in the slice. The DFT prediction matches very well with the FEM prediction, whereas the FFT prediction shows strong Gibbs oscillations in the vicinity of the interfaces between the primary and secondary phases.

Zoom In Zoom Out Reset image size

This modeling approach was further improved into a statistical approach which also account for specimen size effects (Lucarini and Segurado 2020). In this work, a set of random but statistically equivalent polycrystalline RVEs are used to determine the statistical distribution of the FIP corresponding to a given RVE size. Then, the article proposes an upscaling law, based on extreme value distributions, which is able to estimate the FIP distribution of a larger RVE. The result is a fatigue life model which incorporates the specimen size in the life prediction.

The simulation of fracture and damage has a lot of interest in micromechanics in order to capture the effect of the microstructure on the nucleation and development of damage in the material. Fracture is a localized phenomena which involves the development of new free surfaces in the body. Its integration with the FEM implies either effectively creating new surfaces by redefining the domain for each crack propagation step (linear elastic fracture mechanics with crack propagation, (Bittencourt et al 1996)), incorporating special purpose interface elements (cohesive fracture, (Dugdale 1960, Barenblatt 1962)) or enhancing the simulation with strong discontinuities embedded in the element (Sancho et al 2007) or using an auxiliar nodal field as done in the X-FEM approach (Moës et al 1999). The direct extension of these approaches to FFT homogenization is not easy due to the nature of the FFT method which does not alloy to include directly internal surfaces. Nevertheless, some recent works propose some alternatives to incorporate cohesive type approaches in FFT to consider interphase damage and decohesion (Sharma et al 2018).

An alternative view point in fracture problems is to smooth out the crack as a continuous field ϕ concentrated in a thin volume defined by a characteristic length ℓ. The original fracture mechanics problem is recovered when ℓ → 0. These type of models can be classified into the category of variational approaches (Bourdin et al 2008) and they are usually referred to as phase-field fracture models (Bourdin et al 2008, Miehe et al 2010b, 2015). The phase-field fracture approach was developed initially for finite elements but its formulation allows simple implementation in other boundary value numerical methods, as FFT based approaches. Indeed, spectral solvers arise as an ideal framework to implement the model since, as the crack path is defined by a phase field ϕ which evolves through the simulation domain, it does not require special elements or remeshing techniques and a regular grid performs ideally in this case.

The weak form of the phase-field fracture model is obtained applying variational principles and can be written in its most simple form as

$$\begin{equation} {\int }_{\partial {{\Omega}}_{F}}\mathbf{t}\cdot \dot{\mathbf{u}}\enspace \mathrm{d}A={\int }_{{\Omega}}{g}^{\prime }(\phi ){U}_{0}(\boldsymbol{{\epsilon}}(\mathbf{u}))\dot{\phi }+R\left(\frac{1}{\ell }\phi \dot{\phi }+\ell \nabla \phi \cdot \nabla \dot{\phi }\right)-\nabla \cdot \boldsymbol{\sigma }\cdot \dot{\mathbf{u}}\enspace \mathrm{d}{\Omega},\end{equation} \tag{ 138 }$$

where U₀ is the elastic energy and R is the fracture toughness (critical energy released by the creation of a unit crack surface) (Miehe et al 2015) and g(ϕ) is a softening function which takes into account the effect of the diffuse crack on the elastic energy. The function g(ϕ) can have different forms, and the most simple expression is g(ϕ) = (1 − ϕ)² + k with k being a small positive parameter for stabilizing the numerical algorithm. Note that these expressions do not include a history field to ensure irreversibility (Miehe et al 2015) and are therefore valid only for monotonic loading.

The strong form resulting of previous integral is a system of coupled PDEs

$$\begin{equation}\begin{cases}\nabla \cdot \boldsymbol{\sigma }(\mathbf{u},\phi )=0\quad \\ {g}^{\prime }(\phi ){U}_{0}(\boldsymbol{{\epsilon}}(\mathbf{u}))+R\left(\frac{1}{\ell }\phi -\ell {\nabla }^{2}\phi \right)=0\quad \end{cases}\end{equation} \tag{ 139 }$$

with boundary conditions t = σ ⋅ n on ∂Ω_F and ∇ϕ ⋅ n = 0 for ∂Ω_∇ϕ .

Very recently, several adaptations of the phase-field fracture model to FFT approaches have appeared in a short time period (Chen et al 2019b, Ernesti et al 2020, Ma and Sun 2020, Cao et al 2020), illustrating the strong interest of the community in this problem. All the works are focused on brittle fracture as the original phase-field formulation. The use of staggered approaches and discrete Green operators or derivatives are the common numerical receipts for all the implementations. Discrete derivatives are introduced to alleviate the Gibbs effect near shape interfaces or jump conditions. Staggered schemes are chosen since they are the natural integration of multifield problems in FFT and also because it is well known that staggered approaches might provide a more robust way of controlling crack propagation (Ambati et al 2015).

In the first work (Chen et al 2019b) the phase-field fracture model of Miehe et al (Miehe et al 2010b) was adopted and the staggered algorithm in Miehe et al (2010a) was followed to separately solve the mechanical and phase-field problems. Chen et al (Chen et al 2019b) applied the model to study different problems in 2D and 3D including single edge notched tensile specimen with inhomogeneous ℓ or R, crack body under shear loads, symmetric and asymmetric double edge notched tensile specimens, crack branching and coalescence, continuous fiber reinforced composite with void. In the following, the example of crack branching and coalescence studied by Chen et al (Chen et al 2019b) is presented. Figure 27(a) shows the geometry of the unit cell and the material parameters for the two solid parts used for this test. The upper part is 10 times stiffer and tougher than the lower part. Macroscopic tension in the X-axis is applied until failure. Figure 27(b) shows the predicted macroscopic response and figure 27(c) shows the microstructural snapshots at four different strain levels depicted in figure 27(b). At point 1, a crack initiates from the notch, which results in a slope change in the stress–strain curve. Then, at point 2, the crack bifurcates into two branches in the vicinity of the interface of the two solid parts. As the load continues to increase, the two branches penetrate into the upper part and coalesce (point 3), leading to the final failure (point 4). The simulation successfully captures both crack branching and coalescence.

Zoom In Zoom Out Reset image size

Chen et al (2019b) and Cao et al (2020) use the basic scheme to solve the mechanical problem. In Chen et al 2019b a fixed point iteration scheme is also exploited for the phase field evolution and their results clearly show the potential of this approach in a highly parallelized FFT code, solving notched samples and then elastic long fiber composites. In Ma and Sun (2020) the focus is made on simulating the fracture propagation in strongly anisotropic materials as three-dimensional elastic polycrystals. Ernesti et al (2020) focus on implicit solvers and propose the heavy ball algorithm to solve the mechanical problem. This work is applied to the study matrix brittle damage of composites.

Continuum damage is an interesting alternative to fracture models since, contrary to fracture mechanics, it is based a priori in a smeared damage which is represented by a continuous damage field. Continuum damage models are specially interesting for ductile failure due to the disperse nature of the damage within the material, and several continuum damage models have been developed in the last 40 years following more or less phenomenological paths (Gurson 1977, Tvergaard and Needleman 1984, Lemaitre 1985). Nevertheless, the numerical solution of boundary value problems with this class of constitutive laws results in a pathological discretization dependence due to loss of ellipticity of the problem after strain softening. Different regularization techniques have been proposed to overcome this limitation. Among the different approaches, higher-order or implicit gradient approaches enhance the constitutive equations considering a non-local version of some relevant internal variable (equivalent plastic strain ${\bar{\varepsilon }}_{\mathrm{e}\mathrm{q}}$ or porosity in porous plasticity) Peerlings et al (1996). The non-local field is obtained solving the Helmholtz-type PDE using as source function the corresponding local field ɛ_eq,

$$\begin{equation}{\bar{\varepsilon }}_{\mathrm{e}\mathrm{q}}-{\ell }^{2}{\nabla }^{2}{\bar{\varepsilon }}_{\mathrm{e}\mathrm{q}}={\varepsilon }_{\mathrm{e}\mathrm{q}},\end{equation} \tag{ 140 }$$

being ℓ a characteristic length which is related with the damage band width.

The resulting formulation consists of a system of two coupled PDEs: the equation for the mechanical equilibrium—i.e. ∇ ⋅ σ = 0 and the additional equation for computing the non-local field, i.e. equation (140). It is worth nothing that, from a mathematical view point, the implicit gradient regularization resembles phase-field fracture models as a particular case of averaging equation (de Borst and Verhoosel (2016), Steinke et al (2017)).

Spectral solvers become an ideal modeling framework for implicit gradient damage models for the same reasons than phase-field fracture approaches. Nevertheless, only a very small number of studies can be found in the literature which use non-local ductile damage and FFT solvers (Boeff et al 2015, Magri et al 2021). In the work of Boeff et al (2015), an iterative algorithm is proposed for the solution of an implicit gradient regularization of a simple damage model. Their approach was based on a staggered scheme in which the mechanical equation is solved using the basic scheme (Moulinec and Suquet 1994) and the smoothening PDE for solving the damage variable was simplified for an homogenous model. The resulting model was able to prevent the pathological discretization dependency but, due to the poor convergence rate of the basic scheme, the simulations were limited to relatively simple geometries in a two-dimensional setting. The model proposed by Magri et al (2021) also relies in a staggered approach, but the mechanical problem is solved using a Galerkin FFT approach (Vondřejc et al 2014) and the Helmholtz smoothening considers heterogeneous characteristic length, result of having different phases with different damage response. This model uses both GTN (Tvergaard and Needleman 1984) and Lemaitre's (Lemaitre 1985) damage models as local continuum approaches. The resulting framework was able to regularize the macroscopic response and presented a very good numerical performance allowing to simulate very large three dimensional RVEs (figure 28). Moreover, the resolution of the Helmholtz approach in a heterogeneous media allowed to emulate the Neumann-free boundary conditions of the non-local variables in the internal interphases.

Zoom In Zoom Out Reset image size

Focusing the attention on interface damage, Sharma et al (2018) studied, also through gradient damage models, the interface decohesion. Instead of using a sharp interface approach, an interfacial band of finite thickness was modeled under the condition that the total dissipation resulting from the volumetric damage process is equal to the dissipation that would have resulted if a sharp interface (cohesive zone model) model were to be used. Sharma et al (2018) modeled the contribution of damage strain using a non-local damage field, which was regularized in order to minimize the effect of grid spacing. Two regularization techniques were used: integral based averaging (Baz̆ant and Pijaudier-Cabot 1988) and gradient damage based regularization (Peerlings et al 1996). The study conducted by Sharma et al (2018) concluded that the gradient damage approach resulted in the mechanical response and dissipation to be sensitive to the accuracy with which flux boundary conditions were enforced at the edges of the interface band. Meanwhile, the integral approach provided more promising results and it was reported to be straightforward to implement.

Chen et al (2015) proposed an FFT-based staggered algorithm that coupled the VP-FFT model (Lebensohn 2001) and a phase-field recrystallization model based on strain-induced grain boundary migration (SIBM), within which the migration of a pre-existing boundary is initiated via a bulging of a part of the boundary followed by the propagation into a higher dislocation density region leaving behind a dislocation free region (Humphreys et al 2017). At a time step t, the staggered algorithm involves first computing plastic deformation using the VP-FFT model, then updating the grain orientations based on the local rotation rate followed by nucleation of recrystallized grains using the SIBM criterion ${\varepsilon }^{p}(VM)\hspace{-1.5pt} > {\varepsilon }_{\text{crit}}^{p}$ (with ${\varepsilon }_{\text{crit}}^{p}$ being a critical plastic strain value), solving stress equilibrium using the EL-FFT method with the plastic strain introduced as an known eigenstrain, and finally updating the phase field functions using the model of Chen and co-workers Hu and Chen (2001), Chen (2002), Krill and Chen (2002).

Chen et al (Chen et al 2015) first validated their model by studying the deformation of a single crystal and compared with the analytical solution of recrystallization kinetics that comes from the so-called Avrami or JMAK static recrystallization model (Avrami 1939, 1940). Next, stress equilibrium was validated by comparing the von Mises local elastic strain field after nucleation of recrystallization between their proposed FFT-based solver and the theoretical model of Budiansky and Wu Budiansky and Wu (1961). Details on the model parameters and RVE discretization can be found in (Chen et al 2015). The boundary conditions correspond to uniaxial tension along x₃ with an applied strain rate component along the tensile axis ${\dot{E}}_{33}=1{\;\text{s}}^{-1}$. Figure 29 shows microstructural snapshots during the recrystallization simulations at an applied strain of 0.1. The deformed grains are shown in orange, the recrystallized grains in white, and the grain boundaries in black. The recrystallized grains start to grow at sites with relatively high plastic strain, and spread into relatively low plastic strain regions. In particular, the 2D sections show early recrystallization in grains with relatively high stored energy (upper left and middle right) and delayed recrystallization in regions of relatively low stored energy (e.g. lower left grains).

Zoom In Zoom Out Reset image size

Over the years, several thermo-mechanical FFT models have been proposed. One of the first models is the locally linear thermo-elastic model by Vinogradov and Milton (Vinogradov and Milton 2008), implemented via the accelerated scheme of Eyre and Milton (Eyre and Milton 1999). The model was applied to study thermo-elastic (analogously, thermo-electric) composites. In this approach, the stress (or current) polarization tensor has a contribution from the heterogeneous temperature field in addition to the contribution coming from heterogeneous elasticity and strain (charge). The stress polarization tensor reads:

$$\begin{equation}\boldsymbol{\tau }=\left(\mathbb{C}-{\mathbb{C}}^{0}\right):\boldsymbol{\varepsilon }-\boldsymbol{\beta }(T-{T}_{\text{ref}}),\end{equation} \tag{ 141 }$$

where β is the second order thermal moduli tensor, T is the local temperature, T_ref is the reference temperature. The rest of the computation follows similar steps as in the Eyre and Milton (Eyre and Milton 1999) implementation. A basic scheme implementation is also proposed and the convergence analysis reveals that indeed the accelerated implementation converges more rapidly than the basic scheme. The numerical examples however focused on studying the local and macroscopic response of conductive non-linear composites.

Anglin et al (Anglin et al 2014) proposed a heterogeneous thermo-elastic FFT model by introducing the thermal strain as an eigenstrain field in the elastic constitutive response and as a temperature polarization into the expression of the stress polarization field, as done by Vinogradov and Milton Vinogradov and Milton (2008). Anglin et al (2014) implement the model using the basic scheme and validated it against the analytical solutions Eshelby (1957), Mura (1987), Sherer (1986) for spherical and cylindrical inclusions embedded in an infinite medium with same and different elastic properties.

These FFT based thermo-mechanical implementations do not allow imposing standard temperature or heat-flux boundary conditions, which makes it difficult to simulate heat conduction processes. However, it is possible to introduce a homogeneous or heterogeneous distribution of temperature through the eigenstrain approach mentioned above. In the heterogeneous case, however, one must be careful about the periodic boundary conditions.

Ambos et al (Ambos et al 2015) use the thermo-elastic FFT scheme of Vinogradov and Milton (Vinogradov and Milton 2008) to simulate thermal loading by imposing a uniform temperature change everywhere in their RVE. In one of their case studies, a uniform purely thermal loading ΔT = 1 and a purely hydrostatic macroscopic loading $\left\langle {\varepsilon }_{m}\right\rangle =\left\langle \text{tr}(\boldsymbol{\varepsilon })\right\rangle \hspace{-1.5pt}/3$ were applied on two separate instances of the RVE. Figure 30 shows 2D RVE slices in which the mean stress, equivalent strain and stress fields are represented. The highest values for the equivalent strain field occur along the grain boundaries, in particular, near corners or junctions. The equivalent stress field shows the same trend, however, the mean stresses are less localized.

Zoom In Zoom Out Reset image size

More recently, Wicht et al (2021b) proposed a fully-coupled small-strain thermo-mechanical FFT model based on the asymptotic homogenization of thermo-mechanically coupled generalized standard materials of Chatzigeorgiou et al (Chatzigeorgiou et al 2016) and applied it study the thermo-mechanical response of fiber-reinforced composites.

Vidyasagar et al (Vidyasagar et al 2017) were the first ones to propose a fully-coupled electro-mechanics FFT model. The fields implicated in this phenomenon are stress and strain tensors σ, from the mechanical side, and the vector fields of electric displacement d and polarization of the ferroelectric material p from the electric side. The governing equations of the model are summarized below, using index notation due to the high order (up to rank 6) of the tensors involved in the material response,

$$\begin{align}\boldsymbol{\sigma }& =\left(\mathbb{C}+2\left(\mathbb{F}\cdot \boldsymbol{p}\right)\cdot \boldsymbol{p}\right):\boldsymbol{\varepsilon }+\left(\mathbb{B}\cdot \boldsymbol{p}\right)\cdot \boldsymbol{p}+(((\mathbb{G}\cdot \boldsymbol{p})\cdot \boldsymbol{p})\cdot \boldsymbol{p})\cdot \boldsymbol{p}\nonumber\\ \boldsymbol{d}& =\boldsymbol{p}+{\kappa }_{0}\boldsymbol{e}\nonumber\\ \nabla \cdot \boldsymbol{\sigma }& =0\nonumber\\ \nabla \cdot \boldsymbol{d}& =0\nonumber\\ \eta \dot{\boldsymbol{p}}& =\boldsymbol{e}-\frac{\partial {{\Psi}}_{\text{pol}}}{\partial \boldsymbol{p}}-\frac{\partial {{\Psi}}_{\text{coupl}}}{\partial \boldsymbol{p}}+\nabla \cdot (\mathbb{K}:\nabla \boldsymbol{p})\quad \text{with}\nonumber\\ {{\Psi}}_{\text{pol}}& =\frac{1}{2}\boldsymbol{p}\cdot {\mathbb{A}}^{1}\cdot \boldsymbol{p}+\frac{1}{4}\boldsymbol{p}\cdot (\boldsymbol{p}\cdot {\mathbb{A}}^{2}\cdot \boldsymbol{p})\cdot \boldsymbol{p}+\frac{1}{6}\boldsymbol{p}\cdot (\boldsymbol{p}\cdot (\boldsymbol{p}\cdot {\mathbb{A}}^{3}\cdot \boldsymbol{p})\cdot \boldsymbol{p})\cdot \boldsymbol{p}\nonumber\\ & \quad +\frac{1}{8}\boldsymbol{p}\cdot (\boldsymbol{p}\cdot (\boldsymbol{p}\cdot (\boldsymbol{p}\cdot {\mathbb{A}}^{4}\cdot \boldsymbol{p})\cdot \boldsymbol{p})\cdot \boldsymbol{p})\cdot \boldsymbol{p}\quad \text{and}\nonumber\\ {{\Psi}}_{\text{coupl}}& =\boldsymbol{\varepsilon }:((\mathbb{B}\cdot \boldsymbol{p})\cdot \boldsymbol{p})+\boldsymbol{\varepsilon }:((\mathbb{F}\cdot \boldsymbol{p})\cdot \boldsymbol{p}):\boldsymbol{\varepsilon }+\boldsymbol{\varepsilon }:((((\mathbb{G}\cdot \boldsymbol{p})\cdot \boldsymbol{p})\cdot \boldsymbol{p})\cdot \boldsymbol{p}).\end{align} \tag{ 142 }$$

The second-order tensor ${\mathbb{A}}^{1}$, the fourth-order tensors $\mathbb{B}$, $\mathbb{K}$ and ${\mathbb{A}}^{2}$, the sixth-order tensors $\mathbb{F}$, $\mathbb{G}$, and ${\mathbb{A}}^{3}$, and the eight-order tensor ${\mathbb{A}}^{4}$ are material constants defining mechanical, electrical and coupled response. Ψ_pol is the generic multi-stable Landau–Devonshire energy density needed to enforce the symmetric equilibrium states of tetragonal perovskites, Ψ_coupl is the electro-mechanical coupling energy. Equations (142)(a) and (b) are energetic constitutive relationships, equation (142d) comes from Gauss' law, and equation (142e) is a dissipative constitutive relationship. Equations (142)(c)–(e) are solved in Fourier space; an explicit Euler time discretization scheme is used for (142)e. A staggered algorithm is implemented where at each time step the compatible strain fields ɛ_ij are first computed, followed by the electric field e_i and then the ferroelectric polarization field p_i . Gibbs oscillations were corrected using discrete differentiation approach. Then, after some benchmark tests, Vidyasagar et al (Vidyasagar et al 2017) applied the model to study ferroelectric switching. In one of the examples, electric cycling of an RVE is performed using a triangle-wave profile for the bias field at a quasi-static frequency of 0.04 Hz. Figure 31 shows the macroscopic and local fields of a lead zirconate titanate (PZT) polycrystalline RVE whose grain orientations are randomly assigned using a Gaussian profile with a mean at 0° around a mean orientation aligned with the applied electric field $\bar{\boldsymbol{e}}={e}_{2}{\boldsymbol{e}}_{2}$. It is found that the grain boundaries serve as nucleation and stress/electric field concentration sites and influence the hysteresis behavior of the effective response.

Zoom In Zoom Out Reset image size

Other multi-physics coupling FFT frameworks can be found in the literature. For example, Sixto-Camacho et al (2013) proposed a linear thermo-magneto-electro-elastic FFT model for heterogeneous media with periodic and rapidly oscillating coefficients. Rambausek et al (2019) introduced a multi-scale FE–FFT framework for magneto-active materials in a finite-strain setting. Cruzado et al (2021) developed an FFT framework to model the pseudo-elastic behavior and shape-memory effect in phase transforming materials. Meanwhile, Shanthraj et al (2019) proposed an FFT solver for multi-physics simulations (thermo-mechanics and damage) involving crystal plasticity.

High energy x-ray and neutron powder diffraction are important experimental techniques to understand the microstructural evolution of polycrystalline materials during deformation. Specifically, with these experiments one can study the deformation induced evolution of different groups of grains (known as grain families) within a polycrystal. A group of grains, or a grain family, is a set of grains that have one of their {hkl} lattice planes (defined with respect to the crystallographic axes) oriented in such a way that they satisfy the Bragg's condition for diffraction and result in the formation of diffraction patterns that are typically analyzed in the form of {hkl} intensity vs diffraction angle (2θ) peaks. By studying the evolution of different {hkl} peaks during deformation, one can deduce the evolution of different microstructural features. For example, in a deforming single-phase fcc polycrystal, one can follow the evolution of the average {hkl} peak position to deduce information on internal or lattice strains, the evolution of the average {hkl} full widths at half maximum to deduce the dislocation density evolution, and the evolution of {hkl} peak intensities to deduce the texture evolution. Generally, with neutron powder diffraction experiments one can study the microstructural evolution in a larger subset of the material volume than with synchrotron x-ray powder diffraction. Therefore, neutron diffraction is typically used to obtain a statistical understanding of the microstructural evolution, whereas synchrotron x-ray diffraction can be used to study the evolution of a very small subset of grains in the same kind of material or a large subset of grains in materials that have very small grains such as nanocrystalline materials.

At the end of section 4.7, we saw an application of the soft-coupled FE–FFT approach in synergy with biaxial mechanical testing experiments during in-situ neutron diffraction by Upadhyay and co-workers (Upadhyay et al 2016a, 2017a, 2019). However, the first application of an FFT solver being used in synergy with in-situ neutron diffraction experiments was done by Kanjarla et al (2012) who applied the EVP-FFT model in Lebensohn et al (2012) to perform a comprehensive study of the lattice strain evolution during uniaxial loading of dog-bone samples along the longitudinal and one of the transverse loading directions.

High-energy x-ray diffraction microscopy (HEDM) are imaging techniques that can provide information on the grain morphology, grain size, grain position, crystallographic orientation and intragranular strain distributions. Based on the type of x-ray beam used, placement of x-ray detectors and the type of information extracted, HEDM approaches can be classified into the following techniques (Renversade et al 2016): (i) diffraction contrast tomography (DCT) (Ludwig et al 2008, Johnson et al 2008, Ludwig et al 2009): a near-field (nf) technique where a wide box shaped monochromatic beam is used to illuminate a large portion of a (typically cylindrical shaped) sample. (ii) nf-HEDM (Suter et al 2006): a nf-technique where a planar monochromatic beam is used. (iii) Scanning 3D x-ray diffraction (scanning-3DXRD) (Hayashi et al 2015): a nf-technique that uses a pencil-shaped monochromatic beam. (iv) Differential-aperture x-ray microscopy (DAXM) (Larson et al 2002): a nf-technique that uses a polychromatic pencil beam. (v) Far-field 3DXRD (ff-3DXRD) (Poulsen et al 2001, Oddershede et al 2010, Bernier et al 2011): a monochromatic far field (ff) x-ray beam technique involving the use of a pencil beam to study internal strain and stress evolutions. (vi) High resolution reciprocal space mapping (HRRSM) (Jakobsen et al 2006): a ff-technique that uses a monochromatic x-ray pencil beam to study intragranular strain evolution in individual grains.

All these techniques have significant potential to be combined with FFT solvers to better understand the material response, however, only a few studies can be found which combine FFT with the different diffraction techniques. One of the first instances of this synergy was done by Lieberman et al (2016) to study microstructural effects on damage evolution in Cu polycrystals subjected to dynamic shock loading. The microstructure of a Cu polycrystal sample was characterized in its initial undeformed state and post shock loading via nf-HEDM. The initial undeformed state is used as microstructural input for the EVP-FFT model (as shown in figure 34) to study stress localization and identify sites for damage nucleation and compare these with the post shock loading microstructure. EVP-FFT simulations were performed with assistance from FE analysis to generate boundary conditions that were consistent with shock loading. The synergistic study revealed no correlation between damage nucleation and microstructural features and stress measures such as mean and deviatoric stress, stress triaxiality, surface traction or grain boundary inclination angle. However, correlation was found for differences in Taylor factor and accumulated plastic work across grain boundaries with the damage nucleation due to the incipient spall experiment. These insights would not have been possible only from an experimental analysis.

Zoom In Zoom Out Reset image size

A second group of papers which have made extensive use of the synergies between high energy diffraction and FFT solvers is the work performed in Rovinelli et al 2018a, 2018b. The authors used in-situ diffraction and phase contrast tomography to generate high-resolution 3D images of crack propagating into a polycrystalline aggregate. These images were later compared with the predictions of the small-strain based EVP-FFT model (Lebensohn et al 2012) combined with data-science techniques to provide trends on fatigue small crack propagation, as reviewed in section 4.5.1.

The Fourier continuation (FC) framework (Bruno and Lyon 2010, Lyon and Bruno 2010) employs a high-order spectral approach to the numerical analysis of PDEs by enabling fast and accurate interpolating Fourier series representations of non-periodic functions while avoiding the well-known Gibb's phenomenon. This allows high-order computation of spatial derivatives for use in a numerical PDE solver, resulting in algorithms of almost linear complexity that can be applied to general geometries and boundaries. Originally introduced by Bruno and Lyon as a foundation for computing such derivatives as found in general elliptic, hyperbolic and parabolic PDEs (Bruno and Lyon 2010, Lyon and Bruno 2010), the mathematical and computational principles of FC have been extended into a number of full-fledged solvers, including those for realistic solid mechanics applications in the time-domain by Amlani and others (Amlani and Bruno 2016, Amlani et al 2019a, Amlani and Pahlevan 2020). These solvers produce fast, high-order solutions with essentially no numerical dispersion; possess mild CFL constraints on the time step that scale linearly with spatial discretization sizes; and are easily parallelized in distributed-memory computing clusters.

Principles of FC for spatial differentiation

The general idea of FC is to construct an accurate Fourier expansion (see section 2.1.2) of a discretely non-periodic function in space. For example, considering a discretization x_i = i/(N − 1), i = 0, ..., N − 1 of size N for a general smooth function u defined on the unit interval [0, 1], FC algorithms append a small number of points d to the discretized function values u(x_i ) in order to form a (b = 1 + d)-periodic trigonometric polynomial u_cont(y) of the form

$$\begin{equation}{u}_{\text{cont}}(y)=\sum\limits _{k=-M}^{M}{b}_{k}{\mathrm{e}}^{\frac{2\pi \text{i}ky}{1+d}},\end{equation} \tag{ 143 }$$

where u_cont(x_i ) closely matches the original discrete values of u(x_i ). Spatial derivatives can then be computed exactly a term-wise differentiation of (143) as

$$\begin{equation}\frac{\partial u}{\partial y}({x}_{i})=\frac{\partial {u}_{\text{cont}}}{\partial y}({x}_{i})=\sum\limits _{k=-M}^{M}\left(\frac{2\pi \text{i}k}{1+d}\right){b}_{k}\enspace {\mathrm{e}}^{\frac{2\pi \text{i}k{x}_{i}}{L+d}}.\end{equation} \tag{ 144 }$$

Essentially, using a precomputed orthonormal polynomial basis procedure Albin and Bruno (2011), Amlani and Bruno (2016) that employs only a small number d_ℓ , d_r of points at the left and right boundaries of the original function f, FC algorithms add an additional C values to the original discretized function in order to form a periodic extension in [1, 1 + d] that transitions smoothly from u(1) back to u(0) (see figure 35 for a summary of the procedure). The resulting continued function can be viewed as a set of discrete values of a smooth and periodic function which can be approximated to very high-order on a slightly larger interval by a trigonometric polynomial. Once such a function is found, corresponding Fourier coefficients b_k in equation (143) can be obtained rapidly from an application of the FFT. A detailed prescription on accelerated construction of FCs can be found in Amlani and Bruno (2016).

Zoom In Zoom Out Reset image size

An application to ultrasonic non-destructive testing of materials

Ultrasonic non-destructive testing (NDT) is a powerful tool to study the integrity of a variety of plate-like and beam-like structures from aircraft wings to oil pipelines or bridges. Low-energy, high-frequency wave packets are introduced into a material to determine fundamental properties (e.g. elastic constants) or to detect defects (e.g. cracks or holes) by measuring and analyzing the propagation, reflection and attenuation of incident pulses. These pulses are excited by a certain combination of pressure and shear wave modes to enable propagation of a single guided wave. The subsequent dynamics and scattering patterns can then be used to extract important features—such as the positions, dimensions and orientations of defects—by solving the corresponding direct or inverse problems.

Recently, an FC-based solver for linear elastic 3D domains introduced in (Amlani and Bruno 2016) has been utilized in Amlani et al (2019b) to simulate ultrasonic NDT experiments on thin aluminum plates: a pulsed TV-holography (PTVH) system records the 2D acoustic field of the instantaneous out-of-plane displacement over the surface (López-Vázquez et al 2009, 2010). This enables measurements of both the spatial and temporal evolution of incident vibrations that subsequently scatter upon interaction with drilled holes of varying shapes and sizes. Physically-faithful numerical modeling of such configurations entails the use of fully dynamic equations governing elastic waves in a linear, isotropic, possibly heterogeneous medium contained in a general 3D domain Ω. The corresponding PDE is given by

$$\begin{equation}\begin{aligned} \rho (\boldsymbol{x})\frac{{\partial }^{2}\boldsymbol{u}}{\partial {t}^{2}}(\boldsymbol{x},t)& =\nabla \cdot \left[\mu (\boldsymbol{x})\left(\nabla \boldsymbol{u}(\boldsymbol{x},t)+\nabla {\boldsymbol{u}}^{\text{T}}(\boldsymbol{x},t)\right)+\lambda (\boldsymbol{x})\left(\nabla \cdot \boldsymbol{u}(\boldsymbol{x},t)\right)I\right]+\boldsymbol{f}(\boldsymbol{x},t), \\ \quad \boldsymbol{x}& \in {\Omega}\subset {\mathbb{R}}^{3},\quad t\geqslant {t}_{0} \end{aligned}\end{equation} \tag{ 145 }$$

for time-dependent position and displacement vectors x , u ; body force vector f ( x , t); and spatially-varying material properties specified by Lamé parameters μ( x ), λ( x ) and density ρ( x ). Traction (stress-free) boundary conditions are imposed on the surfaces of the hole and plate—an accurate representation of the experimental setup Amlani et al (2019b).

Considering the high-frequency and transient nature of the incident pulses (1 MHz to 10 MHz), this challenging problem has been recently addressed by the FC-based elastic solver of (Amlani and Bruno 2016) using parallel domain decomposition (left images of figure 36). This has enabled a systematic quantitative comparison between numerically simulated maps and filtered experimental displacement maps (Amlani et al 2019b). The results show very good agreement both in amplitude and phase (the L² errors in the boxes of figure 36 fall within 10% error in both amplitude and phase—indicating a very good match within experimental uncertainty Amlani et al (2019b)). This demonstrates the feasibility and applicability of FC-based solvers for the characterization of experimental transient scattering patterns measured with the PTVH technique. The subsequently obtained 3D simulated data, which is often computationally prohibitive to be obtained by dispersive finite element-based methods, has additionally provided a first-time look inside defect holes via numerical analysis of all three displacement components. This has enabled insight into a mode conversion during interaction in the thickness of the hole that is not observable by the experimental setup, nor can be captured by models that have been simplified in the interest of computational expense López-Vázquez et al (2009, 2010).

Zoom In Zoom Out Reset image size

A second group of works studying problems which does not fulfill the strict RVE periodicity have been recently introduced in FFT frameworks Zhou and Bhattacharya (2021), Segurado and Lebensohn (2021). These studies rely on the use Bloch boundary conditions, that were introduced by Bloch (Bloch 1929) to describe the quantum wave functions of electrons in a crystal lattice. The Bloch–Floquet formalism establishes that for any field u traveling as an harmonic wave in a periodic infinite medium with periodicity L, the value of the field at a point x ∈ (−∞, ∞) should follow

$$\begin{equation}u(x,t)=U(x){\mathrm{e}}^{\mathrm{i}kx-\mathrm{i}\omega t},\end{equation} \tag{ 146 }$$

where U(x) is a periodic function with periodicity L, U(x + nL) = U(x) for $n\in \mathbb{Z}$, and k and ω correspond to the wave number and frequency respectively. The key idea is introducing the Bloch condition (equation (146)) in the linear momentum balance to study the deformation of an infinite periodic medium with period L in which the periodicity in the displacement field is fulfilled over several cells, L_period > L. This periodicity is then prescribed by setting the wave number k in equation (146) as k = 2π/L_period. The problem was exploited in Zhou and Bhattacharya (2021) in an static setting, to find bifurcations in the deformation of a non-linear periodic medium. In Segurado and Lebensohn (2021), Bloch conditions were used to study acoustic wave propagation in elastic heterogeneous materials. In particular, an FFT based method was proposed to obtain the dispersion diagram of a periodic heterogeneous material. The accuracy of the method was assessed by comparing with analytical solutions and Fourier series based analysis. The numerical efficiency reached is dependent on the elastic phase contrast, as usual in FFT based approaches, being competitive with FEM for composites or polycrystals. The method was then used to study wave propagation in elastic polycrystals using large 3D RVEs containing hundreds of grains. The simulation allowed to account for the effect of grain size and texture in the wave group velocity.