A fast and robust discrete FFT-based solver for computational homogenization

TL;DR

This work introduces a fast, robust FFT-based solver for computational homogenization using a tetrahedral stencil on regular grids, ensuring mechanical equilibrium is free of nonphysical artefacts across finite and infinite elastic contrasts. By defining all tensor components on a common grid and collapsing two dual discretizations into a single Lippmann-Schwinger equation on full grids, the method enables straightforward FFT implementation and rapid convergence with the basic scheme. Numerical experiments across cubic inclusions, voids, and spherical inclusions demonstrate accurate stress and strain fields without checkerboarding or ringing, outperforming the Moulinec-Suquet and rotated schemes, especially at large contrasts. The approach promises broad applicability to micromechanics and can be extended to related fields requiring robust discrete Green operators and efficient FFT-based solvers.

Abstract

We propose a new discrete FFT-based method for computational homogenization of micromechanics on a regular grid that is simple, fast and robust. The discretization scheme is based on a tetrahedral stencil that displays three crucial properties. First, and most importantly, the Fourier representation of the associated Green operator is defined for any finite q-vector generated by the periodic boundary conditions and that does not belong to the Reciprocal Lattice of the discrete grids. As shown in the paper, this property guaranties that, for any elastic contrats, even infinite, mechanical equilibrium is always mathematically stable, i.e. free of any unphysical patterns, such as oscillations, ringing or checkerboarding, a property which is not shared by the original Moulinec-Suquet method \cite{moulinec1994fast,moulinec1998numerical} nor by the rotated scheme proposed by Willot \cite{willot2015fourier}. Second, the components of tensorial quantities are all defined on the same location, which permits the use of any elastic anisotropy and any spatial variation of the material fields. Third, convergence to equilibrium using the simplest iterative scheme, the "basic scheme", is fast and the number of iterates stabilizes at high contrasts, so that infinite contrast is obtained without additional computational cost.

Figures (10)

• Figure 1: Tetrahedral stencil for a simple cubic discretization. Only two voxels along axis (100) are shown. The simple cubic displacement grid, whose sites are located at the corners of the voxels, is split into two FCC subgrids whose sites are labelled by the indexes 1 and 2, respectively. Similarly, the simple cubic strain grid, whose sites are located at the centers of the voxels, is split into two FCC subgrids whose sites are labelled by the indexes A and B, respectively. The four translations that relate a strain site of type A to the four sites of a tetrahedron of type 1 (blue lines that emerge from a site A) are related through a central symmetry operation to the four translations that relate a strain site of type B to the four sites of a tetrahedron of the same type 1 (red lines that emerge from a site B). This central symmetry relation is reflected in the corresponding finite difference operators $\bm{\underline{D}}^{(A1)}$ and $\bm{\underline{D}}^{(B1)}$ that link the strains on subgrids A and B to the displacements on subgrid 1, as shown in Eq. (\ref{['eq:DB1']}). Similarly, the translation sets that relate strain sites of type A and B to tetrahedra of type 2 (red and blue lines that emerge from sites of type A and B, respectively) are related by the central symmetry operation mentioned above. This is reflected into the relation between the finite difference operators $\bm{\underline{D}}^{(A2)}$ and $\bm{\underline{D}}^{(B2)}$ given in Eq. (\ref{['eq:DA2']}). Finally, we observe that the four translations that relate a strain site A to displacement sites of type 1 are identical to the ones that relate a strain site B to displacement sites of type 2, which implies the equality of the operators $\bm{\underline{D}}^{(A1)}$ and $\bm{\underline{D}}^{(B2)}$, as stated in Eq. (\ref{['eq:DB2']}).
• Figure 2: A cubic inclusion with eigenstrain $\epsilon^0_{33}$ in an isotropic medium with a zero applied stress: internal stresses $\sigma_{11}$, $\sigma_{22}$ and $\sigma_{33}$, in units of $\mu\epsilon^0_{33}$, along the line parallel to direction $(100)$ and that goes through the center of the cell. The numerical results obtained with the tetrahedral stencil (red symbols) are compared to the analytic solution (blue lines); the right panel shows a zoom on the small zone marked in the left panel.
• Figure 3: A cubic void within an isotropic medium subjected to an hydrostatic pressure $P=300$ MPa: error to equilibrium, defined in (\ref{['eq:tolerance']}), as a function of iterates; comparison between the tetrahedral stencil (red), the rotated stencil (green) and Moulinec-Suquet method (blue). The dotted line indicates the precision required to reach, with the tetrahedral stencil, an accuracy of 1 MPa on the stress field.
• Figure 4: A cubic void within an isotropic medium subjected to an hydrostatic pressure $P=300$ MPa : internal stresses $\sigma_{11}$, $\sigma_{22}$ and von Mises stress (MPa) along the line parallel to direction $(100)$ and passing through the center of the void. The numerical results obtained with the tetrahedral stencil (red lines) are compared to the results obtained with the rotated scheme (green lines) and with the Moulinec-Suquet scheme (blue lines); the cubic cavity consists in $31^3$ voxels; full lines represent results obtained with a unit cell containing $64^3$ voxels and doted lines those obtained with a unit cell containing $65^3$ voxels; the middle and bottom panels show zooms on the zones marked 1 and 2, respectively, in the top panel.
• Figure 5: A cubic void within an isotropic medium subjected to an hydrostatic pressure $P=300$ MPa : von Mises stress maps along the (001) plane passing through the center of computational box; comparaison between the results obtained with the Moulinec-Suquet method (top left panel), the rotated scheme (top right panel) and the tetrahedral solver (bottom panel). The computational box is discretized into $64^3$ cells and the cubic void consists in $31^3$ voxels .