Contents

Computational homogenization turns a heterogeneous microstructure β€” glass spheres embedded in epoxy, fibres in a polymer, voids in metal β€” into a single effective stiffness tensor. FFT-based methods, introduced by Moulinec and Suquet in the mid-1990s, achieve this without a mesh: the microstructure is represented on a regular voxel grid, and the governing equations are solved spectrally. The result is an algorithm that is easy to implement, naturally periodic, and fast enough to run in a browser using WebAssembly.

1. Interactive Demo

Controls (top-right panel): set volume fraction and grid resolution, choose load case and field to display, pick a slice plane, then press Run FFT. The solver runs the full 3D Moulinec-Suquet iteration in Rust/WASM. The 3D view shows the RVE cube with the cut plane and axis arrows; the lower panels compare $E^\mathrm{eff}$ and $G^\mathrm{eff}$ against all analytical bounds.

Note on shear load cases: load cases $\varepsilon_{23}, \varepsilon_{13}, \varepsilon_{12}$ return the shear modulus $G^\mathrm{eff}$, plotted separately in the second chart.

Press Run FFT to solve

Left-drag: orbit  ·  Scroll: zoom  ·  Right-drag: pan  ·  Blue = min  ·  White = mid  ·  Red = max

2. The Representative Volume Element

Consider a composite occupying a periodic domain $\Omega = [0,L]^3$ (the representative volume element, or RVE). Each point $\mathbf{x} \in \Omega$ belongs to one of two phases with local stiffness $\mathbf{C}(\mathbf{x})$.

The volume average of any field $\mathbf{f}$ is

\[\langle \mathbf{f} \rangle = \frac{1}{|\Omega|} \int_\Omega \mathbf{f}(\mathbf{x})\, d\mathbf{x}\]

The goal of homogenization is to find the effective stiffness $\mathbf{C}^\text{eff}$ such that

\[\langle \boldsymbol{\sigma} \rangle = \mathbf{C}^\text{eff} : \langle \boldsymbol{\varepsilon} \rangle\]

This is the macro-scale constitutive law that a continuum model would use. The double contraction $\mathbf{A}:\mathbf{B} = A_{ijkl}B_{kl}$ contracts over two indices.

Periodic boundary conditions

Periodicity is the key assumption. The displacement field decomposes as

\[\mathbf{u}(\mathbf{x}) = \mathbf{E} \cdot \mathbf{x} + \mathbf{u}^*(\mathbf{x})\]

where $\mathbf{E}$ is the prescribed macroscopic strain and $\mathbf{u}^*$ is a $\Omega$-periodic fluctuation with zero mean. The total strain is then

\[\boldsymbol{\varepsilon}(\mathbf{x}) = \mathbf{E} + \boldsymbol{\varepsilon}^*(\mathbf{x}), \qquad \langle \boldsymbol{\varepsilon} \rangle = \mathbf{E}\]

Computing $\mathbf{C}^\text{eff}$

Apply six unit load cases $\mathbf{E}^{(k)}$ (one for each independent Voigt component). For each, solve the equilibrium problem and compute $\langle \boldsymbol{\sigma} \rangle$. The $k$-th column of $\mathbf{C}^\text{eff}$ is $\langle \boldsymbol{\sigma} \rangle^{(k)}$.

3. The Reference Medium and Polarization

Instead of solving the heterogeneous problem directly, introduce a homogeneous reference medium with stiffness $\mathbf{C}^0$ (chosen as the arithmetic mean of the two phase stiffnesses). Define the polarization stress

\[\boldsymbol{\tau}(\mathbf{x}) = \bigl(\mathbf{C}(\mathbf{x}) - \mathbf{C}^0\bigr) : \boldsymbol{\varepsilon}(\mathbf{x})\]

The local stress can now be written as

\[\boldsymbol{\sigma}(\mathbf{x}) = \mathbf{C}^0 : \boldsymbol{\varepsilon}(\mathbf{x}) + \boldsymbol{\tau}(\mathbf{x})\]

Equilibrium $\nabla \cdot \boldsymbol{\sigma} = \mathbf{0}$ becomes

\[\nabla \cdot \bigl(\mathbf{C}^0 : \boldsymbol{\varepsilon}\bigr) = -\nabla \cdot \boldsymbol{\tau}\]

The left-hand side is the equilibrium operator for the homogeneous reference medium. Its Green’s function in Fourier space is the Green operator $\hat{\boldsymbol{\Gamma}}^0(\boldsymbol{\xi})$.

4. The Lippmann-Schwinger Equation

Inverting the equilibrium operator gives the Lippmann-Schwinger equation:

\[\boldsymbol{\varepsilon}(\mathbf{x}) + \bigl(\boldsymbol{\Gamma}^0 * \boldsymbol{\tau}\bigr)(\mathbf{x}) = \mathbf{E}\]

where $*$ denotes periodic convolution. This is the central equation of FFT homogenization: the actual strain equals the macroscopic strain minus the strain induced by the heterogeneity polarization, transmitted through the reference-medium Green operator.

Substituting the polarization definition gives a fixed-point equation for $\boldsymbol{\varepsilon}$:

\[\boldsymbol{\varepsilon} + \boldsymbol{\Gamma}^0 * \bigl[(\mathbf{C} - \mathbf{C}^0):\boldsymbol{\varepsilon}\bigr] = \mathbf{E}\]

5. The Green Operator in Fourier Space

The convolution becomes a pointwise multiplication in Fourier space:

\[\hat{\boldsymbol{\varepsilon}}(\boldsymbol{\xi}) + \hat{\boldsymbol{\Gamma}}^0(\boldsymbol{\xi}) : \hat{\boldsymbol{\tau}}(\boldsymbol{\xi}) = \hat{\mathbf{E}}(\boldsymbol{\xi})\]

For an isotropic reference medium with Lame parameters $\lambda^0, \mu^0$, the Green operator in the Fourier domain is

\[\hat{\Gamma}^0_{ijkl}(\boldsymbol{\xi}) = \frac{1}{4\mu^0} \bigl(\delta_{ik}\hat{n}_j\hat{n}_l + \delta_{il}\hat{n}_j\hat{n}_k + \delta_{jk}\hat{n}_i\hat{n}_l + \delta_{jl}\hat{n}_i\hat{n}_k\bigr) - \frac{\lambda^0+\mu^0}{\mu^0(\lambda^0+2\mu^0)}\,\hat{n}_i\hat{n}_j\hat{n}_k\hat{n}_l\]

where $\hat{\mathbf{n}} = \boldsymbol{\xi}/\lvert\boldsymbol{\xi}\rvert$ is the unit wave vector. This is singular at $\boldsymbol{\xi} = \mathbf{0}$ (the DC component), which is simply set to zero β€” this enforces $\langle\boldsymbol{\varepsilon}^*\rangle = \mathbf{0}$, i.e., the fluctuation has zero mean.

In Voigt notation ${11,22,33,23,13,12}$ the Green operator becomes a $6\times6$ matrix. Off-diagonal shear entries acquire factors of 2 (normal-shear) or 4 (shear-shear) from the symmetry of the strain tensor.

6. The Moulinec-Suquet Algorithm

The fixed-point iteration leads directly to the basic scheme of Moulinec and Suquet (1994):

Physical interpretation:

  • Polarization: compute how much the local stiffness deviates from the reference medium, weighted by the current strain guess.
  • Forward FFT: transform to Fourier space where the convolution with $\boldsymbol{\Gamma}^0$ is cheap.
  • Green operator: each frequency mode of the polarization induces a strain correction via the reference-medium response.
  • Inverse FFT: bring the strain correction back to real space.
  • Update: the new total strain is macroscopic strain plus fluctuation.

Convergence is linear with rate $\rho = (C_\mathrm{max}/C_\mathrm{min} - 1)/(C_\mathrm{max}/C_\mathrm{min} + 1)$, where $C_\mathrm{max}/C_\mathrm{min}$ is the phase contrast. For glass/epoxy (${\approx}20\times$) this gives $\rho \approx 0.9$ β€” slow but guaranteed convergent.

7. Worked Example: 2Γ—2Γ—2 Toy RVE

To show exactly what is computed at each step, consider a $2\times2\times2$ voxel RVE (8 voxels) with one inclusion voxel at position $(0,0,0)$ and seven matrix voxels.

Material parameters (GPa):

  • Matrix: $E_m = 3.5$, $\nu_m = 0.35$ gives $\lambda_m = 3.03$, $\mu_m = 1.30$
  • Fiber: $E_f = 72$, $\nu_f = 0.22$ gives $\lambda_f = 23.2$, $\mu_f = 29.5$

Reference medium (arithmetic mean): $\lambda^0 = 13.1$ GPa, $\mu^0 = 15.4$ GPa

Load case: unit macroscopic strain $\mathbf{E} = \mathbf{e}_1 \otimes \mathbf{e}_1$, Voigt: $[1,0,0,0,0,0]$.

Iteration 0 to 1

Step 1: Polarization $\boldsymbol{\tau}^{(0)}$

All voxels start with $\boldsymbol{\varepsilon}^{(0)} = [1,0,0,0,0,0]$.

For a matrix voxel, with $C^m_{11} = \lambda_m + 2\mu_m = 5.63$ GPa and $C^0_{11} = \lambda^0 + 2\mu^0 = 43.9$ GPa:

\[\tau^\mathrm{matrix}_{11} = (C^m_{11} - C^0_{11}) \cdot 1 = 5.63 - 43.9 = -38.3 \text{ GPa}\] \[\tau^\mathrm{matrix}_{22} = \tau^\mathrm{matrix}_{33} = (C^m_{12} - C^0_{12}) \cdot 1 = 3.03 - 13.1 = -10.1 \text{ GPa}\]

For the fiber voxel at $(0,0,0)$, with $C^f_{11} = 23.2 + 2(29.5) = 82.2$ GPa:

\[\tau^\mathrm{fiber}_{11} = 82.2 - 43.9 = +38.3 \text{ GPa}, \qquad \tau^\mathrm{fiber}_{22} = 23.2 - 13.1 = +10.1 \text{ GPa}\]

The polarization is equal and opposite in the two phases β€” a direct consequence of the arithmetic-mean reference: $C^0 = (C^m + C^f)/2$.

Step 2: Forward FFT of $\boldsymbol{\tau}$

A $2\times2\times2$ DFT has frequencies $\boldsymbol{\xi} \in {0,\pi}^3$. The DC component of $\tau_{11}$ is:

\[\hat{\tau}_{11}(\mathbf{0}) = 1 \cdot (+38.3) + 7 \cdot (-38.3) = -229.8 \text{ GPa}\]

Step 3: Apply Green operator

At $\boldsymbol{\xi} = \mathbf{0}$: $\hat{\boldsymbol{\Gamma}}^0 = \mathbf{0}$ (DC set to zero), so no strain correction at DC.

At $\boldsymbol{\xi} = (\pi, 0, 0)$, the unit normal is $\hat{\mathbf{n}} = (1,0,0)$, giving:

\[\hat{\Gamma}^0_{1111} = \frac{1}{\lambda^0 + 2\mu^0} = \frac{1}{43.9} = 0.0228 \text{ GPa}^{-1}\]

The strain correction:

\[\hat{\varepsilon}^*_{11}(\pi,0,0) = -0.0228 \times \hat{\tau}_{11}(\pi,0,0)\]

and similarly for all 8 frequency points.

Step 4: Update $\boldsymbol{\varepsilon}^{(1)} = \mathbf{E} + \boldsymbol{\varepsilon}^*$

After the inverse FFT, the strain in the stiff fiber voxel drops below 1 and in the soft matrix rises above 1 β€” physically correct: the stiff inclusion attracts stress and sheds strain. After ${\approx}40$ iterations the residual falls below $10^{-4}$ and the average stress gives

\[C^\mathrm{eff}_{11} \approx \langle\sigma_{11}\rangle.\]

8. Analytical Bounds

Voigt and Reuss

The Voigt bound (iso-strain, upper) and Reuss bound (iso-stress, lower), where $V_f$ is the fiber (inclusion) volume fraction and $V_m = 1 - V_f$:

\[E^\mathrm{Voigt} = V_m E_m + V_f E_f, \qquad \frac{1}{E^\mathrm{Reuss}} = \frac{V_m}{E_m} + \frac{V_f}{E_f}\]

The true modulus lies strictly between these for a two-phase composite.

Hashin-Shtrikman Bounds

The tightest possible bounds for a statistically isotropic composite are the Hashin-Shtrikman bounds (1963). In terms of bulk modulus $K$ and shear modulus $G$:

Upper bounds (stiff phase as reference):

\[K^+ = K_f + \frac{V_m}{(K_m - K_f)^{-1} + \dfrac{3V_f}{3K_f + 4G_f}}\] \[G^+ = G_f + \frac{V_m}{(G_m - G_f)^{-1} + \alpha_f}, \qquad \alpha_f = \frac{6V_f(K_f+2G_f)}{5G_f(3K_f+4G_f)}\]

Lower bounds (soft phase as reference):

\[K^- = K_m + \frac{V_f}{(K_f - K_m)^{-1} + \dfrac{3V_m}{3K_m + 4G_m}}\] \[G^- = G_m + \frac{V_f}{(G_f - G_m)^{-1} + \alpha_m}, \qquad \alpha_m = \frac{6V_m(K_m+2G_m)}{5G_m(3K_m+4G_m)}\]

Convert back to Young’s modulus: $E = 9KG/(3K+G)$.

Mori-Tanaka

The Mori-Tanaka model gives a single point estimate. For spherical inclusions in a softer matrix it coincides with the HS lower bound (Benveniste 1987):

\[K^\mathrm{MT} = K_m + \frac{V_f(K_f-K_m)}{1 + V_m(K_f-K_m)/(K_m + 4G_m/3)}\] \[G^\mathrm{MT} = G_m + \frac{V_f(G_f-G_m)}{1 + V_m(G_f-G_m)/(G_m + f_m)}, \qquad f_m = \frac{G_m(9K_m+8G_m)}{6(K_m+2G_m)}\]

Physical picture: each inclusion is embedded in a matrix subjected to the average matrix strain, accounting for inclusion-inclusion shielding.

See also: Analytical bounds reference β€” all models (Voigt, Reuss, HS, MT, Halpin-Tsai, dilute, self-consistent), derivations, interactive comparison charts, and DOI-linked references.

9. Material System and Results

The demo uses glass inclusions ($E_f = 72$ GPa, $\nu_f = 0.22$) in an epoxy matrix ($E_m = 3.5$ GPa, $\nu_m = 0.35$), a contrast ratio of about $20\times$.

For normal load cases ($\varepsilon_{11}$, $\varepsilon_{22}$, $\varepsilon_{33}$), the FFT result for a random isotropic arrangement of spheres should land between the HS bounds and close to the Mori-Tanaka / HS lower bound. Deviation from the lower bound grows at higher volume fractions as inclusion-inclusion interactions become significant.

For shear load cases ($\varepsilon_{23}$, $\varepsilon_{13}$, $\varepsilon_{12}$), the returned quantity is the effective shear modulus $G^\mathrm{eff} = \langle\sigma_{ij}\rangle / \gamma_{ij}$ where $\gamma_{ij} = 1$ is the applied unit engineering shear. These should also fall within the HS bounds for $G$, shown in the second chart.