# Numerical analysis of several FFT-based schemes for computational homogenization

We study the convergences of several FFT-based schemes that are widely applied in computational homogenization for deriving effective coefficients, and the term "convergence" here means the limiting behaviors as spatial resolutions going to infinity. Those schemes include Moulinec-Suquent's scheme [Comput Method Appl M, 157 (1998), pp. 69-94], Willot's scheme [Comptes Rendus Mécanique, 343 (2015), pp. 232-245], and the FEM scheme [Int J Numer Meth Eng, 109 (2017), pp. 1461-1489]. Under some reasonable assumptions, we prove that the effective coefficients obtained by those schemes are all convergent to the theoretical ones. Moreover, for the FEM scheme, we can present several convergence rate estimates under additional regularity assumptions.

• 4 publications
• 20 publications
02/22/2021

### Convergence error estimates at low regularity for time discretizations of KdV

We consider various filtered time discretizations of the periodic Kortew...
02/18/2020

### Convergence and error estimates for a finite difference scheme for the multi-dimensional compressible Navier-Stokes system

We prove convergence of a finite difference approximation of the compres...
07/12/2022

### Hybrid Physical-Neural ODEs for Fast N-body Simulations

We present a new scheme to compensate for the small-scales approximation...
01/25/2021

### L^p-Convergence Rate of Backward Euler Schemes for Monotone SDEs

We give a unified method to derive the strong convergence rate of the ba...
10/08/2019

### An Analysis of the Milstein Scheme for SPDEs without a Commutative Noise Condition

In order to approximate solutions of stochastic partial differential equ...
06/29/2022

### Implicit and fully discrete approximation of the supercooled Stefan problem in the presence of blow-ups

We consider two implicit approximation schemes of the one-dimensional su...
04/23/2022

### Unboundedness for Recursion Schemes: A Simpler Type System

Decidability of the problems of unboundedness and simultaneous unbounded...

## 1 Introduction and preliminaries

In an inhomogeneous material that obeys the linear elasticity law, the fourth-order stiffness tensor

varies at different . The central goal of computational homogenization methods is determining an effective constitute law, which may serve as a surrogate model for downstream applications Milton (2002); Zohdi and Wriggers (2008).

Let be a Representative Volume Element (RVE) where through the whole paper, and be the fourth-order stiffness tensor for any . We denote by the closure of smooth -periodic functions with respect to ordinary -norm Cioranescu and Donato (1999), and define accordingly. The notation will be replaced as or

for vector-valued functions, while

(, ) should be always understood as -norm (-norm, -norm) of on the domain regardless of whether is vector-valued or not. For a discrete set , let be the linear space of discrete functions which take inputs in and return values of , where may be , , , and (the set of symmetric matrices). The double dot operator “” and the single dot operator “” are defined conventionally for any order tensors (including vectors and matrices), and note that now the matrix-vector product is rewritten as . A popular method to derive the effective coefficients of is solving a periodic boundary value problem: find with such that for any

 ⟨∇s\vectorsymv:\tensorsymC(\vectorsymx):∇s\vectorsymu⟩=−⟨∇s\vectorsymv:\tensorsymC(\vectorsymx)⟩:\matrixsymE, (1)

where stands for , belongs to and . The effective stiffness tensor can be obtained by

 \tensorsymCeff:\matrixsymE=⟨\tensorsymC(\vectorsymx):(∇s\vectorsymu+\matrixsymE)⟩ (2)

via choosing different .

The starting point of FFT-based homogenization is rewriting the above variational form into the Lippmann-Schwinger equation. Taking as a linear elastic reference medium, which is a constant fourth-order tensor satisfying certain coercive conditions, we can introduce a Green function such that for any , the convolution is the solution to the variational form

 ⟨∇s\vectorsymv:\tensorsymC% ref:∇s(\tensorsymG0∗\matrixsymF)⟩=⟨∇s\vectorsymv:\matrixsymF⟩,∀\vectorsymv∈H1#(Y;Rd).

By setting and taking the Fourier series expansion of as , we can derive that for

 ˆ\tensorsymG0∗\matrixsymF[\vectorsymξ]=i2π∣∣\vectorsymξ∣∣2μ0⎡⎣(μ0+λ0)\vectorsymξ⊗\vectorsymξ(2μ0+λ0)∣∣\vectorsymξ∣∣2−\matrixsymId⎤⎦⋅ˆ\matrixsymF⋅\vectorsymξ,

where is the identity matrix and is the Kronecker product. Moreover, it is easy to show that eq. 1 is equivalent to

 \vectorsymu=−\tensorsymG0∗(δ\tensorsymC:(∇s\vectorsymu+\matrixsymE)+\tensorsymC% ref:\matrixsymE)=−\tensorsymG0∗(δ\tensorsymC:(∇s\vectorsymu+\matrixsymE)),

where and vanishes due to that is a constant matrix. If taking and , we can immediately read the above equation as

 \matrixsymε+\tensorsymΓ0∗(δ\tensorsymC:\matrixsymε)=\matrixsymE (3)

the so-called Lippmann-Schwinger equation Lippmann and Schwinger (1950). The alternative form of eq. 3 in the Fourier space is

Specifically, when , we can derive a closed form of as

 mnpq[\vectorsymξ]= −(λ0+μ0)\vectorsymξm\vectorsymξn\vectorsymξp\vectorsymξqμ0(λ0+2μ0)∣∣\vectorsymξ∣∣4

Now it comes to discretization schemes. For simplicity, we equally split into parts in every dimension. For any -dimensional index with , we denote by

 Y\vectorsymI=(\vectorsymI1N,\vectorsymI1+1N)×(\vectorsymI2N,\vectorsymI2+1N)×(\vectorsymI3N,\vectorsymI3+1N)

an element in the language of Finite Element Methods (FEM), which could be interpreted as a pixel (for 2D problems) or voxel (for 3D problems). Correspondingly, the frequency domain is denoted by

. We always postulate that an image-like (i.e., is constant in every ) rather that the ground truth is provided, and introduce a notation where is the geometric center of . The basic scheme of FFT-based homogenization was proposed by Moulinec and Suquent in Moulinec and Suquet (1995, 1998), and can be stated as algorithm 1. Depending on the type of a function’s domain, the Fourier coefficients are defined differently, i.e. if , then is determined by the formula

 f(\vectorsymx)=∑\vectorsymξ∈Zdˆf[\vectorsymξ]exp(2πi\vectorsymξ⋅\vectorsymx),

and if , then satisfies

 f⋆[\vectorsymI]=∑\vectorsymξ∈FNˆf⋆[\vectorsymξ]exp(2πi\vectorsymξ⋅\vectorsymIN).

Note that if is labeled with a superscript “”, the function should be understood as a discrete one with as its domain.

Moulinec-Suquent’s scheme (we also call it the basic scheme) gains considerable popularity. The reasons can be summarized as follows: 1) microstructures in heterogeneous materials may be rather complex, which causes that generating high quality meshes dominates overall computational overheads Hughes et al. (2005); 2) the primary information of microstructures is usually provided by modern digital volume imaging techniques Larson et al. (2002); Poulsen (2004); Landis and Keane (2010), and we may not be capable to retrieve the original geometric descriptions and to perform a preprocessing for FEMs; 3) the easy implementation of those schemes and highly optimized FFT packages (e.g., Intel®MKL, FFTW Frigo and Johnson (2005)) secure the global efficiency; 4) generally, it is the effective coefficients we are more interested in rather than the local displacement or the local stress, and unfitted structured meshes may have less influence on desired comparing to and . Moreover, those methods are not limited in linear elasticity, versatile applications such as nonlinear elasticity Moulinec and Suquet (1998); Michel et al. (2001), piezoelectricity Brenner (2009), damage and fracture mechanics Zhu and Yvonnet (2015); Chai et al. (2020), polycrystalline materials Segurado et al. (2018) can be found in literature. We strongly recommend a recent review article Schneider (2021) as an exhaustive reference for historical developments and the current state of the art of FFT-base homogenization methods.

Willot’s scheme is another popular discretization method of the Lippmann-Schwinger equation Willot (2015). The main idea in Willot’s scheme is replacing the gradient operator in the Fourier space with

 \vectorsymkN[\vectorsymξ]=iN4d∏m=1{exp(2πi\vectorsymξmN)+1}[tan(π\vectorsymξ1N),tan(π\vectorsymξ2N),tan(π\vectorsymξ3N)] (4)

that is derived from a finite difference stencil. If is isotropic with the Lamé coefficients , the closed form of is

 mnpq[\vectorsymξ]= 14μ0∣∣\vectorsymkN∣∣2{[\vectorsymkN]m[\widebar\vectorsymkN]qδnp+[\vectorsymkN]n[\widebar\vectorsymkN]qδmp +[\vectorsymkN]m[\widebar\vectorsymkN]pδnq+[\vectorsymkN]n[\widebar\vectorsymkN]pδmq} −(λ0+μ0)[\vectorsymkN]m[\widebar\vectorsymkN]n[\vectorsymkN]p[\widebar\vectorsymkN]qμ0(λ0+2μ0)∣∣\vectorsymkN∣∣4,

and this is the major difference from the basic scheme. Take

 FN\shortminus\coloneqq{\vectorsymξ∈Zd:−N/2<\vectorsymξm

and Willot’s scheme states as algorithm 2.

It has been numerically demonstrated that convergences of Moulinec-Suquent’s scheme will deteriorate as the contrast ratio of grows, and an extreme example shows that the basic scheme indeed fails to converge for a porous material Schneider et al. (2016), while Willot’s scheme stays a stable convergence history in high contrast settings Willot (2015). Meanwhile, computational overheads of the basic and Willot’s schemes are almost same. Those advantages make Willot’s scheme become another standard method in discretizing the Lippmann-Schwinger equation eq. 3.

Interestingly, it is pointed in Schneider et al. (2017) that Willot’s scheme is related to a reduced integration variational problem: find with such that for any ,

 N−d∑\vectorsymI∈IN∇s\vectorsymvN(\vectorsymx\vectorsymI):\tensorsymC⋆N[\vectorsymI]:∇s\vectorsymuN(\vectorsymx\vectorsymI)=−N−d∑\vectorsymI∈IN∇s\vectorsymvN(\vectorsymx\vectorsymI):\tensorsymC⋆N[\vectorsymI]:\matrixsymE,

where is the trilinear Lagrange finite element space of . The reduced integration technique is significantly efficient in constructing stiffness matrices comparing to full integration, while the latter needs evaluations on all the Gaussian quadrature points ( for trilinear elements) in each element. However, due to that the stiffness matrix derived by reduced integration may be singular, there exist so-called “hourglassing” or “zero-energy” modes Koh and Kikuchi (1987); Pugh et al. (1978). For example when is an even integer, it is easy to find nonzero such that for all , and this is closely related with the frequency , which also implies there does not exist a positive constant independent of satisfying

 c∥∥\vectorsymvN∥∥2H1(Y)≤N−d∑\vectorsymI∈IN∇s\vectorsymvN(\vectorsymx\vectorsymI):\tensorsymC⋆N[\vectorsymI]:∇s\vectorsymvN(\vectorsymx\vectorsymI).

In Schneider et al. (2017), Schneider et al. proposed a novel scheme by replacing reduced integration in Willot’s scheme with full integration, i.e., find with such that for any ,

 (2N)−d∑\vectorsymI∈IN∑\vectorsymb∈{−1,1}d∇s\vectorsymvN(\vectorsymx\vectorsymI\vectorsymb):\tensorsymC⋆N[\vectorsymI]:∇s\vectorsymuN(\vectorsymx\vectorsymI\vectorsymb)=⟨∇s\vectorsymvN:\tensorsymCN:∇s\vectorsymuN⟩ = −N−d∑\vectorsymI∈IN∇% s\vectorsymvN(\vectorsymx\vectorsymI):\tensorsymC⋆N[\vectorsymI]:\matrixsymE=−⟨∇s\vectorsymvN:\tensorsymCN⟩:\matrixsymE,

where is the Gaussian quadrature point on the element determined by symbols of . The FFT technique is utilized in solving reference problem: find with such that for any ,

 (2N)−d∑\vectorsymI∈IN∑\vectorsymb∈{−1,1}d∇s\vectorsymvN(\vectorsymx\vectorsymI\vectorsymb):\tensorsymCref:∇s\vectorsymuN(\vectorsymx\vectorsymI\vectorsymb)=⟨∇s\vectorsymvN:\tensorsymCref:∇s\vectorsymwN⟩ = (2N)−d∑\vectorsymI∈IN∑\vectorsymb∈{−1,1}d∇s\vectorsymvN(\vectorsymx\vectorsymI\vectorsymb):\matrixsymF\vectorsymb,⋆N[\vectorsymI],

where . The motivation is converting into

via the Discrete Fourier Transform (DFT), where

is the symmetric Kronecker product and . Combining the symmetric relations of , we can hence derive an explicit formula of in the Fourier space as for , where and . The scheme is summarized in algorithm 3, we called it “the FEM scheme” for it is essentially a FEM with the FFT acting as a preconditioner Saad (2003).

Currently, most researches on FFT-based homogenization are focused on developing fast algorithms to accelerate convergences of iterating processes and designing schemes to stabilize performances on high/infinite contrast materials Eyre and Milton (1999); Michel et al. (2001); Zeman et al. (2010); Eloh et al. (2019); Schneider (2020), while few of them consider the convergences of those methods with respect to the spatial resolution , which could be viewed as an analogy of -estimate theories in FEMs Ciarlet (1991); Brenner and Scott (2008). To our knowledge, the work by Schneider Schneider (2015) is only published one which seriously discusses such a question. In his article, a priori error estimate of Moulinec-Suquet’s scheme is formulated with the trigonometric interpolation operator Zygmund (1968). However, because of involving pointwise evaluations, the trigonometric interpolation operator is to some certain “incompatible” with Lebesgue integrable functions.

Since the main objective in this article is studying convergences with respect to spatial resolutions, it is necessary to clarify the relation between the provided information and the ground truth . As mentioned previously, the ground truth is hidden from the scheme part, we hence cannot hold an optimistic anticipation on convergence rates like immersed FEMs Li et al. (2003); Chen et al. (2009) or unfitted FEMs Burman et al. (2014); Huang et al. (2017); Chen et al. (2021). The assumptions for and is elucidated as follows.

### Assumption A

The RVE domain consists of Lipschitz subdomains that are labeled by , and only takes a constant tensor in each subdomain .

### Assumption B

The constant tensors satisfy properties: (symmetricity) for any ,

 [\tensorsymΛl]mnpq=[\tensorsymΛl]nmpq=[\tensorsymΛl]mnqp=[\tensorsymΛl]pqmn,∀1≤m,n,p,q≤d;

(coercivity) there exist positive constants and such that for any and ,

 Λ′\matrixsymF:\matrixsymF≤\matrixsymF:\tensorsymΛl:\matrixsymF≤Λ′′\matrixsymF:\matrixsymF.

### Assumption C

For any , if then , else will be arbitrarily chosen from .

In some cases, we will use the notation to represent a positive constant which depends on .

1. In section 2, we rebuild the theories in Schneider (2015) on the convergence of Moulinec-Suquet’s scheme by introducing a new operator which is will-defined for Lebesgue integrable functions.

2. In section 3, we prove the convergence of the effective coefficients obtained by Willot’s scheme.

3. In section 4, by assuming some suitable regularities and combining several priori estimates in FEM theories, we present convergence rates of the solution and the effective coefficients derived by the FEM scheme.

## 2 Convergence of Moulinec-Suquent’s scheme

Let be a complex trigonometric polynomial space

 ⎧⎨⎩f(\vectorsymx)=∑\vectorsymξ∈FNc[\vectorsymξ]exp(2πi\vectorsymξ⋅\vectorsymx):c[\vectorsymξ]∈C⎫⎬⎭,

and be the orthogonal projection onto Conway (1990). Take where the set is the cube with as its center and as its edge length. In the following analysis, for a function , the convolution should be understood as applying to the periodic extension of , which induces that

 GN∗f(\vectorsymx)=∑\vectorsymξ∈ZdˆGN[\vectorsymξ]ˆf[\vectorsymξ]exp(2πi\vectorsymξ⋅\vectorsymx),

where

 ˆGN[\vectorsymξ]=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩d∏m=1Nsin(π\vectorsymξmN)π\vectorsymξm,if d∏m=1\vectorsymξm≠0,1,otherwise. (5)

It is easy to see that for , and we can hence define an operator that plays an essential role in our analysis:

###### Definition 1.

For any ,

 QNf(\vectorsymx)=∑\vectorsymξ∈FNˆGN−1[\vectorsymξ]exp(−πi\vectorsymξ⋅\vectorsym1N)ˆf⋆[\vectorsymξ]exp(2πi\vectorsymξ⋅\vectorsymx),

where and

 f⋆[\vectorsymI]=GN∗f(\vectorsymx\vectorsymI)=\fintY\vectorsymIf(\vectorsymx)d\vectorsymx.

We can show the following facts of :

###### Proposition 1.

The following statements hold true for

• Let , then , and if then ;

• Let , then and

 ∥f−QNf∥L2(Y)≤{1+(π/2)d}∥f−PNf∥L2(Y);
• for a series of with , then .

###### Proof.

From the definition of , it follows that . Moreover, if , we have

 f⋆[\vectorsymI] =∑\vectorsymξ∈FNˆGN[\vectorsymξ]ˆf[\vectorsymξ]exp(2πi\vectorsymξ⋅\vectorsymx\vectorsymI) =∑\vectorsymξ∈FNˆf⋆[\vectorsymξ]exp(2πi\vectorsymξ⋅\vectorsymIN).

Then holds by taking into the definition of .

From the fact for , we have

 ∥QNf∥2L2(Y) =∑\vectorsymξ∈FN∣∣∣ˆGN−1[\vectorsymξ]∣∣∣2∣∣ˆf⋆[\vectorsymξ]∣∣2 ≤(π2)2d∑\vectorsymξ∈FN∣∣ˆf⋆[\vectorsymξ]∣∣2 =(π2)2d1Nd∑\vectorsymI∈IN∣∣f⋆[\vectorsymI]∣∣2 (Parseval's theorem) =(π