FFT 方法求解单胞问题
参考文献：MOULINEC H, SUQUET P. A numerical method for computing the overall response of nonlinear composites with complex microstructure[J/OL]. Computer Methods in Applied Mechanics and Engineering, 1998, 157(1-2): 69-94. DOI:10.1016/S0045-7825(97)00218-1.
力学列式
对于非均质线弹性材料，其模量为 \(\mathbb{L}:Y \mapsto \mathbb{S}^4\)，单胞上的应变与应力为 \(\varepsilon, \sigma:Y \mapsto \mathbb{S}^2\)，在单胞上满足无体力的平衡方程： \[ \nabla \cdot \sigma = 0 \tag{1} \] 以及线弹性本构方程 \[ \sigma = \mathbb{L} : \varepsilon \tag{2} \] 引入一个参考均质模量 \(\mathbb{L}^0 \in \mathbb{S}^2\)，并定义 \[ \delta \mathbb{L} \triangleq \mathbb{L}- \mathbb{L}^0 \] 以及极化应力 \(p:Y \mapsto \mathbb{S}^2\) \[ p \triangleq \delta \mathbb{L}:\varepsilon \] 代入本构方程（2）中得到 \[ \sigma = p + \mathbb{L}^{0}:\varepsilon \] 为了定解上述方程，我们将应变 \(\varepsilon(x)\) 分解为远场处的常值应变 \(\varepsilon_0 \in \mathbb{S}^2\) 和满足周期性边界条件的应变扰动 \(\tilde{\varepsilon}:\Omega \mapsto \mathbb{S}^2\) ： \[ \varepsilon \triangleq \tilde{\varepsilon} + \varepsilon^0 \] 由此得到 \[ \sigma = p + \mathbb{L}^{0}:\tilde{\varepsilon} + \mathbb{L}^{0}:\varepsilon^{0} \] 再代入到平衡方程（1）中得到 \[ \nabla \cdot \left( \mathbb{L}^{0}:\tilde{\varepsilon} \right) + \nabla \cdot p = 0 \] 写成关于位移 \(u:Y \mapsto \mathbb{R}^3\) 的分量形式 \[ L_{ijkl}^0 \tilde{u}_{k,lj} + p_{ij,j} = 0 \tag{3} \] 若求得偏微分方程的基本解 \[ L_{ijkl}^0 G_{nk,lj}(x-x^{\prime}) + \delta_{in} \delta(x-x^{\prime}) = 0 \] 那么位移场可以表示为基本解和极化应力的卷积形式。对式（3）进行 Fourier 变换得到 \[ \xi_j\xi_l L_{ijkl}^{0}~ \hat{\tilde u}_{k}(\xi) = i\xi_j~ \hat{p}_{ij}(\xi) \tag{4} \] 这里用到的 Fourier 变换定义为 \[ \mathscr{F}[f(x)](\xi) = \hat{f}(\xi) \triangleq \int_{\mathbb{R}^3} f(x) e^{-i \xi \cdot x} \mathrm{d}x \] 以及相应的 Fourier 逆变换 \[ \mathscr{F}^{-1}[f(\xi)](x) = \check{f}(\xi) \triangleq \frac{1}{(2\pi)^3}\int_{\mathbb{R}^3} f(\xi) e^{i \xi \cdot x} \mathrm{d}\xi \] 定义关于代数方程（4）的系数矩阵以及其它相关量 \[ K_{ik} \triangleq \xi_j\xi_l L_{ijkl}, \quad N_{ik} \triangleq \mathrm{adj}(K_{ik}), \quad D \triangleq \det(K_{ik}) \] 由此方程（4）的解可以表示为 \[ \hat{\tilde u}_{k}= i\xi_j~ D^{-1} N_{ki} \hat{p}_{ij} \]
为了接下来可以得到具有对称性的算子，我们利用 \(\hat{p}_{ij}\) 的对称性将上式写为 \[ \hat{\tilde u}_{k} = \frac{i}{2}~ D^{-1} \left( \xi_j N_{ki} + \xi_i N_{kj} \right) \hat{p}_{ij} \] 所以应变扰动 \(\tilde{\varepsilon}\) 的 Fourier 变换为 \[ \hat{\tilde\varepsilon}_{ij} = \frac{i}{2} (\xi_i \hat{\tilde u}_{j} + \xi_j \hat{\tilde u}_{i}) = - \frac{1}{4} D^{-1} \left( \xi_l N_{ik} + \xi_k N_{il} + \xi_l N_{jk} + \xi_k N_{jl} \right) \hat{p}_{kl} \tag{5} \] 因此我们通过上式将极化应力与应变扰动相关联。为了写成更简洁的形式，定义 \(\hat{\Gamma}_{ijkl}\) 为 \[ \hat{\Gamma}_{ijkl} \triangleq \frac{1}{4} D^{-1} \left( \xi_l N_{ik} + \xi_k N_{il} + \xi_l N_{jk} + \xi_k N_{jl} \right) \] 所以式（5）可以写成如下更加紧凑的形式： \[ \hat{\tilde\varepsilon}_{ij} = - \hat{\Gamma}_{ijkl} \hat{p}_{kl} \] 或者写成卷积的形式： \[ \tilde\varepsilon_{ij} = - \Gamma_{ijkl} * p_{kl} \] 最终我们得到关于应变 \(\varepsilon\) 的积分方程： \[ \varepsilon_{ij} = \varepsilon_{ij}^{0} - \Gamma_{ijkl} * (\delta L_{klmn} \varepsilon_{mn}) \]
若参考材料为各项同性，其 Lame 系数和剪切模量分别为 \(\lambda^0\) 和 \(\mu^0\)，那么可以显式地写出相空间中 \(\hat{\Gamma}_{ijkl}^0\)： \[ \hat{\Gamma}_{ijkl}^{0} =\frac{1}{4 \mu^0|\xi|^2} \left(\delta_{ik} \xi_j \xi_l + \delta_{il} \xi_j \xi_k +\delta_{jk} \xi_i \xi_l + \delta_{jl} \xi_i \xi_k \right) -\frac{\lambda^0+\mu^0}{\mu^0\left(\lambda^0+2 \mu^0\right)} \frac{\xi_i \xi_j \xi_k \xi_l}{|\xi|^4} . \]
算法
| 初始化 | \(\varepsilon^{(0)}(x) = \varepsilon^0, \quad \sigma^{(0)}(x) = \mathbb{L}(x):\varepsilon^{(0)}(x), \quad \forall x \in Y\) |
| 第 \(i+1\) 次迭代 | 已知 \(\varepsilon^{(i)}\) 与 \(\sigma^{(i)}\) |
| (a) | \(p^{(i)}(x) = \sigma^{(i)}(x) - \mathbb{L}_0 : \varepsilon^{(i)}(x)\) |
| (b) | \(\hat{p}^{(i)} = \mathscr{F}[p^{(i)}]\) |
| (c) | 收敛性测试，\(e^{(i)} < \mathrm{Tol}\) |
| (d) | \(\hat{\varepsilon}^{(i+1)}(\xi) = -\hat{\Gamma}^{0}(\xi):\hat{p}^{(i)}\) |
| (e) | \(\varepsilon^{(i+1)}(\xi) = \mathscr{F}^{-1}[\hat{\varepsilon}^{(i+1)}(\xi)]\) |
| (f) | \(\sigma^{(i+1)} = \mathbb{L}(x):\varepsilon^{(i+1)}(x)\) |
当得到的应力场 \(\sigma^{(i+1)}\) 满足平衡方程时算法收敛，因此检验算法收敛性的误差公式为 \[ e^{(i)} = \frac{\left( \langle \| \nabla \cdot \sigma^{(i)} \|^2 \rangle \right)^{\frac{1}{2}}} {\langle \| \sigma^{(i)} \| \rangle} \]