# FFT-accelerated modeling of mechanical metamaterials: A solver for torsion and bending

Indre Jödicke, Martin Ladecký, Lars Pastewka

<figure><img src="image_1.png" /></figure>

## To cite this version:

Indre Jödicke, Martin Ladecký, Lars Pastewka. FFT-accelerated modeling of mechanical metamaterials: A solver for torsion and bending. 2025. hal-05387359

HAL Id: hal-05387359

https://hal.science/hal-05387359v1

Preprint submitted on 28 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# FFT-accelerated modeling of mechanical metamaterials:
A solver for torsion and bending

Indre Jödicke,¹ Martin Ladecký,¹,² and Lars Pastewka¹,²,∗

¹Department of Microsystems Engineering, University of Freiburg,
Georges-Köhler-Allee 103, 79110 Freiburg, Germany

²Cluster of Excellence livMatS, Freiburg Center for Interactive
Materials and Bioinspired Technologies, University of Freiburg,
Georges-Köhler-Allee 105, 79110 Freiburg, Germany

(Dated: November 2025)
Abstract

Fast Fourier transform (FFT)-accelerated solvers provide a computationally efficient framework for mechanical modeling of periodic representative volumes of material microstructures. These solvers model the deformation of a heterogeneous unit cell under uniform macroscopic strain. In this study, we use an FFT-accelerated finite element discretization capable of simulating voids with zero stiffness. These void regions decouple periodic images, enabling the modeling of free (non-periodic) boundaries and thus the application of non-periodic and non-uniform strain fields required for torsion and bending. We validate our small-strain formulation through benchmark tests on beam torsion and bending. We demonstrate the capabilities of the solver on a known chiral metamaterial and compare the results with conventional finite element simulations. Our solver enables efficient simulation of complex metamaterial architectures under realistic loading conditions, bridging the gap between homogenization theory and practical engineering applications.

## I. INTRODUCTION

Mechanical metamaterials are “materials” whose behavior is largely determined by their structure as opposed to by the mechanical properties of the constituent material [1–3]. They can exhibit properties that are impossible to realize in bulk materials, such as extremely high stiffness-density ratios [4], negative Poissons ratios [5–8] and negative compressibility [9]. In the last decade, chiral mechanical metamaterials have been proposed which allow for the direct coupling of a torque and a force along the same axis [10–13].

The overwhelming majority of metamaterials are composed of periodically repeated unit cells. To characterize the metamaterial, the effective properties of one unit cell are often studied with some homogenization method. One class of efficient solvers for homogenization problems employs the fast Fourier transform (FFT) for solution [14]. However, traditional FFT solvers are limited to periodic boundary conditions, precluding application to torsion and bending—two fundamental loading modes essential for characterizing metamaterial performance. This work overcomes this limitation.

The FFT-based solvers for mechanical equilibrium originate in the works of Moulinec and Suquet in the 1990s [15, 16]. Since then, FFT-based solvers have greatly evolved

<footnote>* Corresponding author: lars.pastewka@imtek.uni-freiburg.de</footnote>
26 (see Ref. [17–19] for current reviews) and have been successfully used for complex homog-
27 enization problems, e.g. Refs. [20, 21]. The main advantages of FFT-based solvers include:
28 i) excellent convergence rates that are independent of system size; ii) quasi-linear compu-
29 tational complexity enabled by the FFT algorithm and its efficient implementations, such
30 as FFTW [22]; iii) greater ease of implementation relative to conventional finite element
31 methods.

During the last decade, a variant of the method using *compatibility projection* was de-
32 veloped by Vondrejc, Zeman, de Geus et al. [23–25]. These solvers are often classified as
33 strain-based, i.e. they work with some strain measure (e.g. the deformation gradient) rather
34 than the displacements as the unknown variable. To ensure that the strain field is indeed a
35 (symmetrized) gradient of displacement, the strain needs to fulfill an additional condition:
36 the compatibility condition. In strain-based FFT solvers, compatibility can be enforced by
37 a projection operator. The operator is straightforward to evaluate in a Fourier space, and
38 no additional compatibility equation must be solved. It has been shown that the projection
39 serves also as a preconditioner, ensuring the good convergence properties of the solver [26].
40

Since all Fourier-based solvers rely on FFT, they have implicit periodic boundary con-
ditions. While this is often an advantage in homogenization problems that are usually
41 formulated in terms of periodically repeated representative volume elements, it also limits
42 the number of applicable problems to those that can be represented with Bravais lattices.
43
44 A problem outside this category is the application of torsion or bending to mechanical
45 structures.
46

Recently, methods have been proposed to extend the possible boundary conditions of
47 FFT-based solvers to include Dirichlet boundary conditions (e.g. [27, 28]), Neumann bound-
48 ary conditions (e.g [29]), mixed boundary conditions (e.g. [30]) and imposed velocities
49 (e.g [31]). Here, we will propose an approach, which relies on modeling free surfaces, to
50 model torsion with FFT-based solvers. This approach has the advantage that it is straight-
51 forward and easy to implement.
52

Our method is based on an FFT-accelerated formulation that uses a strain-based, finite-
53 element projection at its core recently proposed by us [32]. The discretization is completely
54 equivalent to canonical finite-elements, and the FFT-acceleration can be interpreted as pre-
55 conditioning [26]. Our formulation allows modeling voids of zero stiffness [32], which in
56 more traditional formulations would lead to artificial oscillations [33] that persist in many
57
58 FFT-based solvers (see e.g. Refs.[34, 35]). Zero-stiffness regions can be used to introduce
59 free surfaces (see e.g. Refs. [36, 37]). In this paper, we show how to exploit the possibility
60 to model free surfaces to allow torsional and bending deformation of a microstructure. The
61 approach was inspired by efforts to break the symmetry imposed by Bravais lattices in clas-
62 sical [38] and quantum [39–41] atomic modeling. We demonstrate our solver on the torsion
63 of chiral metamaterials and bending of beams.

64 The ability to model torsion and bending within FFT frameworks is significant because:
65 (i) FFT solvers offer 10–100 times speedup over conventional FEM for metamaterial homog-
66 enization; (ii) extending them to non-periodic boundary conditions unlocks analysis of chiral
67 and other architecturally complex metamaterials; and (iii) the computational efficiency en-
68 ables rapid design iteration for metamaterial optimization.

## 69 II. METHOD

### 70 A. Mechanical equilibrium

71 Our goal is to simulate periodic material micro-structures in a three-dimensional compu-
72 tational domain $\Omega_0 = [0, l_x] \times [0, l_y] \times [0, l_z]$, where $l_x, l_y$, and $l_z$ are lengths of the domain
73 in $X, Y$, and $Z$ directions, respectively. To do so, we solve for static mechanical equilibrium
74 under the assumption of small strains, which means solving for nonzero strain field $\varepsilon(\vec{X})$
75 that satisfy the equation

$$ \nabla \cdot \boldsymbol{\sigma}(\varepsilon(\vec{X})) = \vec{0}, \quad \text{for all } \vec{X} \in \Omega_0, \qquad (1) $$

76 where $\nabla \cdot$ is the gradient with respect to the undeformed positions $\vec{X}$ and the operator $\cdot$ is
77 the dot product: $[\vec{a} \cdot \mathbf{B}]_j = a_i B_{ij}$. Here, and in what follows, we use the Einstein summation
78 convention. We use bold symbols to represent second-order tensors and arrows to represent
79 vectors. $\boldsymbol{\sigma}(\vec{X})$ is the Cauchy stress tensor, which depends on the local small-strain $\varepsilon(\vec{X})$
80 through an appropriate constitutive law. Note that $\boldsymbol{\sigma}$ and $\varepsilon$ are symmetric in the small
81 strain formulation.

82 We assume that the stress field $\boldsymbol{\sigma}(\vec{X})$ is a $\Omega_0$-periodic function: $\boldsymbol{\sigma}(\vec{X}) = \boldsymbol{\sigma}(\vec{X} + n_1\vec{a}_1 + n_2\vec{a}_2 + n_3\vec{a}_3)$ for integer $n_i \in \mathbb{Z}$. Here, three vectors $\vec{a}_1, \vec{a}_2$, and $\vec{a}_3$ define the unit cell $\Omega_0$.
83 The Fourier basis, which automatically fulfills the periodicity condition, is a natural basis
84 function for the expansion of $\boldsymbol{\sigma}(\vec{X})$ in analytic and numerical techniques.
85
86 **Weak equilibrium.**

87 We use the finite element discretization, therefore we proceed to the weak form by intro-
88 ducing the test vectors $\vec{v}$. Using the divergence theorem, we obtain the weak form of the
89 problem as:

$$
0 = \int_{\Omega_0} d^3 X \, \vec{v} \cdot \nabla \cdot \boldsymbol{\sigma} \left( \boldsymbol{\varepsilon}(\vec{X}) \right) = \underbrace{\int_{\partial\Omega_0} d^3 X \, \vec{v} \cdot \boldsymbol{\sigma} \left( \boldsymbol{\varepsilon}(\vec{X}) \right) \cdot \hat{n}}_{=0} - \int_{\Omega_0} d^3 X \, (\nabla \otimes \vec{v}) \cdot \boldsymbol{\sigma} \left( \boldsymbol{\varepsilon}(\vec{X}) \right), \quad (2)
$$

90 where $\partial\Omega_0$ is the surface of $\Omega_0$, $\hat{n}$ is the normal vector to that surface, and the operator : designates the double contraction $\mathbf{A} : \mathbf{B} = A_{ij}B_{ji}$. The operator $\otimes$ denotes the outer product, 91 i.e. $[\vec{a} \otimes \vec{b}]_{ij} = a_i b_j$. Due to the periodicity of the stress $\boldsymbol{\sigma}$, the integral over the boundary 92 $\partial\Omega_0$ is zero. We emphasize that this remains valid even when the stress at the boundary 93 is zero, which can occur because of the zero stiffness of the material in the vicinity of the 94 boundary. 95

96 **Compatibility projection.** The term $\nabla \otimes \vec{v}$ represents a compatible test strain field.
97 To use an arbitrary test strain field $\delta\boldsymbol{\varepsilon}(\vec{X})$, we must enforce the compatibility. We achieve
98 this by using a self-adjoint projection operator $\mathcal{G}$, which maps arbitrary strain fields into
99 the compatible one. Accordingly, the equation (2) takes the form

$$
0 = \int_{\Omega_0} d^3 X \, \delta\boldsymbol{\varepsilon}(\vec{X}) : \mathcal{G} * \boldsymbol{\sigma}(\boldsymbol{\varepsilon}(\vec{X})), \quad (3)
$$

100 where * denotes convolution. For a detailed derivation, we refer interested reader to papers
101 by Leute et al. [32] and de Geus et al. [25]. The projection operator inherently depends on the
102 chosen discretization scheme. We use the finite element discretization with five tetrahedral
103 elements per voxel and linear shape functions [26, 32]. Finite element discretization allows
104 us to model voids, i.e., material with zero stiffness, which we will use to model internal voids
105 in metamaterials and to break the periodicity.

106 **B. Strain**

107 In our simulations, we impose a reference strain field $\boldsymbol{\varepsilon}_0(\vec{X})$, such that

$$
\boldsymbol{\varepsilon}(\vec{X}) = \boldsymbol{\varepsilon}_0(\vec{X}) + \Delta\boldsymbol{\varepsilon}(\vec{X}), \quad (4)
$$

108 where $\Delta\boldsymbol{\varepsilon}(\vec{X})$ is a field with zero mean describing fluctuations around the reference strain
109 $\boldsymbol{\varepsilon}_0(\vec{X})$.
110 Equivalently, we can formulate the mechanical equilibrium in terms of the displacement field $\vec{u}(\vec{X})$, because the strain $\varepsilon$ is the symmetrized gradient of the displacement $\vec{u}$,

$$ \varepsilon = \frac{1}{2} \left[ \nabla \otimes \vec{u} + (\nabla \otimes \vec{u})^T \right]. \qquad (5) $$

112 In analogy to the strain decomposition in Eq. (4), the displacement field can be written as

$$ \vec{u}(\vec{X}) = \vec{u}_0(\vec{X}) + \Delta\vec{u}(\vec{X}), \qquad (6) $$

113 where $\vec{u}_0(\vec{X})$ is a reference displacement field that yields the reference strain $\varepsilon_0(\vec{X})$. The unknown displacement fluctuation $\Delta\vec{u}(\vec{X})$ is related to the strain fluctuation $\Delta\varepsilon(\vec{X})$ by Eq. (5); 114 both are $\Omega_0$-periodic. In contrast, the reference displacement field $\vec{u}_0(\vec{X})$ is not periodic. 115

116 **Linear displacement.** Traditionally in micromechanics (homogenization), we prescribe 117 a linear reference displacement

$$ \vec{u}_0(\vec{X}) = \bar{\varepsilon} \cdot \vec{X}, $$

118 where $\varepsilon_0(\vec{X}) = \bar{\varepsilon}$ is a spatially uniform (constant) reference strain field.

119 **Imposed displacement.** In our simulations, we aim to explore more complex deformation scenarios. Therefore, we write the reference displacements as

$$ \vec{u}_0(\vec{X}) = \bar{\varepsilon} \cdot \vec{X} + \vec{u}_{\text{imposed}}(\vec{X}), $$

121 where $\vec{u}_{\text{imposed}}(\vec{X})$ is an imposed displacement field that may, for instance, characterize 122 bending or torsional deformation of the domain. The corresponding reference strain field is 123 given by:

$$ \varepsilon_0(\vec{X}) = \bar{\varepsilon} + \varepsilon_{\text{imposed}}(\vec{X}). $$

124 **Linearization.** Based on the prescribed decomposition of the strain field, the weak form 125 of micromechanical equilibrium Eq. (3) yields

$$ 0 = \int_{\Omega_0} d^3 X \delta\varepsilon(\vec{X}) : \mathbb{G} * \sigma(\underbrace{\bar{\varepsilon} + \varepsilon_{\text{imposed}}}_{=\varepsilon_0}) + \Delta\varepsilon. \qquad (7) $$

126 This is similar to the weak formulation in other papers like Refs. [25, 32]. Newton-Raphson 127 iteration iteratively solves a linearized problem of the form

$$ \mathbb{G} * \frac{\partial \sigma}{\partial \varepsilon}(\varepsilon_0 + \Delta\varepsilon_i) : \delta\Delta\varepsilon_{i+1} = -\mathbb{G} * \sigma(\varepsilon_0 + \Delta\varepsilon_i), \qquad (8) $$
<figure><img src="image_3.png" /><figcaption>FIG. 1. Beam torsion schematic. (a) A 3D illustration shows the twisting deformation of a beam with a square cross-section. The beam rotates around the axis $\hat{n}$, centered at point $\vec{X}_0$ within the bar's section. The green arrows represent the displacements of the corners. Column (x.1) shows top (b.1) and side (c.1) view of the undeformed beam. Voids (elements with zero stiffness) are positioned around the beam's geometry to break the periodic boundary conditions. Panels (b.2) and (c.2) show the deformed beam.</figcaption></figure>

<footnote n="128">128 for an unknown strain fluctuation increment $\delta\Delta\varepsilon_{i+1}$. Here $\frac{\partial\sigma}{\partial\varepsilon}(\varepsilon_0 + \Delta\varepsilon_i)$ is the stiffness
129 tangent at the current strain. The strain fluctuation is then obtained iteratively, such that
130 $\Delta\varepsilon_{i+1} = \Delta\varepsilon_i + \delta\Delta\varepsilon_{i+1}$ where $i$ counts the iteration. We solve the linearized Eq. (8) using a
131 conjugate-gradient solver.</footnote>

<footnote n="132">132 **Breaking periodicity.** To accommodate non-periodic imposed strains, we introduce
133 void layers (zero-stiffness regions) at the domain boundaries. This decouples the periodic
134 images: since stress vanishes in void regions, $\sigma = 0$, the periodicity assumption is satisfied
135 even though the imposed strain is non-periodic elsewhere. The key enabling feature is our
136 solver's capability to handle zero-stiffness elements exactly [32].</footnote>
### 137 C. Torsion

One possible application of this approach is the modeling of torsional deformation. For torsion centered at a position $\vec{X}_0$ around an axis with normal $\hat{n}$ as represented in Fig. 1 (a), we can write the reference displacements as

$$ \vec{u}_0(\vec{X}) = \bar{\epsilon} \cdot \vec{X} + \vec{u}_{\text{twist}}(\vec{X}), \qquad (9) $$

with

$$ \vec{u}_{\text{twist}}(\vec{X}) = w \left[ \hat{n} \otimes \hat{n} \cdot (\vec{X} - \vec{X}_0) \right] \times \left[ (\mathbf{1} - \hat{n} \otimes \hat{n}) \cdot (\vec{X} - \vec{X}_0) \right], \qquad (10) $$

where $\times$ is the cross product, $\mathbf{1}$ the second order identity tensor and $w$ the twist angle per length. Eq. (10) follows the right hand rule, i.e. positive $w$ leads to counter-clockwise twist.

The corresponding strain field is given by:

$$ \epsilon_0 = \bar{\epsilon} + \epsilon_{\text{twist}} = \bar{\epsilon} + \frac{w}{2} \left( \hat{n} \otimes \left( \hat{n} \times (\vec{X} - \vec{X}_0) \right) + \left( \hat{n} \times (\vec{X} - \vec{X}_0) \right) \otimes \hat{n} \right). \qquad (11) $$

We will twist around the z-axis of our simulation cell, such that

$$ \vec{u}_{\text{twist}}(\vec{X}) = \begin{pmatrix} -w(Y - Y_0)Z \\ w(X - X_0)Z \\ 0 \end{pmatrix}, \qquad (12) $$

where $X, Y, Z$ are the components of $\vec{X}$ and we set $Z_0 = 0$. The strain field we impose during simulation is therefore:

$$ \epsilon_{\text{twist}}(\vec{X}) = \frac{w}{2} \begin{pmatrix} 0 & 0 & -(Y-Y_0) \\ 0 & 0 & X-X_0 \\ -(Y-Y_0) & X-X_0 & 0 \end{pmatrix}. \qquad (13) $$

As discussed in Sec. II B, the reference strain $\epsilon_0(\vec{X})$ is not $\Omega_0$-periodic in $X$ and $Y$ direction while considering twist. In consequence, it is necessary to break the periodic boundary conditions in $X$ and $Y$ direction direction. This, however, does not effect the stress field $\sigma$ at these parts of the boundary, because we place a layer of void pixels at these boundaries where $\sigma \equiv 0$.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Beam bending schematic. (a) A 3D illustration shows the bending deformation of a beam with a rectangular cross-section. The prescribed bending radius is denoted by $R$, and $\bar{z}_0$ represents the neutral axis of bending. Column (x.1) shows the undeformed state of the beam: (b.1) displays the side view, while (c.1) shows the cross-sectional view. The right column shows the deformed beam: (b.2) displays the side view, and (c.2) shows the corresponding cross-section. Voids (elements with zero stiffness) are positioned around the beam's geometry to break the periodic boundary conditions.</figcaption></figure>

## D. Bending

The presented method can be easily extended to problems other than torsion. The only condition is that the applied deformation must be described using a strain field which is periodic in at least one single direction. As an example, we consider the plane stress bending of an Euler-Bernoulli beam with a constant bending moment around the y-axis as shown in Fig. 2 (a). The strain representing this load case is given by:

$$ \epsilon_{\text{bending}} = \begin{pmatrix} \kappa(Z - Z_0) & 0 & 0 \\ 0 & -\nu\kappa(Z - Z_0) & 0 \\ 0 & 0 & -\nu\kappa(Z - Z_0) \end{pmatrix}, \qquad (14) $$

with the Poisson's ratio of the material $\nu$ and the imposed curvature at the origin $\kappa = 1/R$.

Here, $Z_0$ denotes the prescribed neutral axis of the imposed bending deformation: the axis
162 with zero imposed strain. Note that the actual neutral axis of the beam, with vanishing longitudinal stress may differ depending on its geometry. In the bending simulations described 163 below, we break the periodicity in z-direction and y-direction. Therefore the stiffnesses in 164 these directions are zero, which means that of the three non-zero entries of $\varepsilon_{\text{bending}}$ only 165 $\varepsilon_{\text{bending},xx}$ is relevant. 166

### III. RESULTS

168 In this section we study some cases to validate our method. All examples were calculated 169 using the FFT-based solver µSpectre [42]. We have used the five tetrahedral finite element 170 per voxel with linear shape functions [26, 32] for discretization. In all examples, the imposed 171 uniform strain was zero: $\bar{\varepsilon} = 0$. In each case, the material was assumed to be linear elastic 172 and isotropic.

173 We considered a three-dimensional computational domain $\Omega_0 = [0, l_x] \times [0, l_y] \times [0, l_z]$, 174 where $l_x, l_y$, and $l_z$ are lengths of the domain in X, Y, and Z directions, respectively. Note 175 that, if we do not break the periodic boundary conditions in direction $i$, $l_i = L_i$, where $L_i$ 176 is the length of the geometry in direction $i$. In contrast, if we break the periodic boundary 177 conditions in direction $i$, $l_i = L_i + \Delta L_i$, where $\Delta L_i$ is the length of void in direction $i$.

#### A. Validation with analytical solutions

179 First we studied three cases for which analytical solutions exist: (1) the torsion of a beam 180 with a circular cross-section, (2) the torsion of a beam with a square cross-section, and (3) 181 the bending of a beam with a rectangular cross-section.

182 **Torsion.** As parameters for the torsion problems we have used a beam of a homogeneous 183 material with a Poisson's ratio of $\nu = 0$. The radius of the circular cross-section is $r = 0.04L_z$ 184 and the length of the square cross-section is $L_x = L_y = 0.06L_z$, where the length of the beam 185 in z-direction $L_z$ serves as the length unit. The beam was loaded with a twist angle per 186 length of $w = 1.5/L_z$. The number of voxels in z-direction is $n_z = 30$.

187 Figure 3 illustrates the discretization error by comparing the xz-component of the analytical stress with the xz-component of the numerical stress for a discretization of 110x110x30 188 voxels. It shows that the numerical stress agrees well with the analytical one. Only at some 189
<figure><img src="image_5.png" /><figcaption>FIG. 3. Stress component $\sigma_{xz}$ over the cross-section of a cylinder subjected to torsion. Since $\sigma_{xz}$ is symmetric, only one half of the cross-section is shown. (a) shows the analytical stress and (b) the numerical stress.</figcaption></figure>

191 voxels at the border of the cylinder is the numerical stress significantly higher or lower than 192 the analytical stress. This difference is due to the fact that cubic voxels are not able to 193 fully represent a circular cross-section. For the other components of stress tensor $\sigma$ we find 194 similar results.

195 The analytical torsion stiffness $S_{ana}$ is given by:

$$S_{ana} = \mu J, \qquad (15)$$

196 where $\mu$ is the shear modulus and $J$ is a torsion constant depending on the cross-section. 197 We chose to rotate around the geometric center of the cross-section. In this case, the torsion 198 constant of a circular cross-section is

$$J_{\text{circular}} = \frac{1}{2}\pi r^4 \qquad (16)$$

199 with $r$ the radius. For a square cross-section, no exact analytical solution exists but the 200 approximation is

$$J_{\text{square}} \approx 0.141a^4, \qquad (17)$$

201 where $a$ denotes the length of the square (see e.g. Ref. [43, Chapter 2.5]).

202 The numerical torsion stiffness was calculated from the numerical stresses via the average 203 numerical moment:

$$S_{\text{num}} = \frac{M_{\text{num}}}{w} = \frac{1}{w} \frac{1}{L_z} \int_{0}^{L_z} dZ \int_{0}^{L_x} dX \int_{0}^{L_y} dY \sigma_{yz}(X - X_0) - \sigma_{xz}(Y - Y_0) \qquad (18)$$
<figure><img src="image_6.png" /><figcaption>FIG. 4. Relative error of the numerical stiffness with respect to the analytical stiffness for torsion and bending. (a) Error for torsion of two beams with different cross-sections. The number of voxels $n$ corresponds to the number of voxels in $x$-direction and in $y$-direction $n = n_x = n_y$. (b) Error for bending of one beam with rectangular cross-section. The number of voxels $n$ corresponds to the number of voxels in each cartesian direction $n = n_x = n_y = n_z$.</figcaption></figure>

Figure 4 (a) shows the difference between the analytical stiffnesses and the numerical stiffnesses for different discretizations. The differences of the torsion stiffnesses do not monotonically decrease. However, if the discretization is fine enough (50x50x30 voxels or finer), the difference is less than 1.5%. For a square cross-section, which can be discretized exactly by our voxel grid, the difference seems to converge towards around 0.3%. Note that we get similar results for different material parameters or different geometric parameters.

That the differences for the torsional stiffness do not disappear completely is expected: First, for a circular cross-section, cubic voxels are not able to fully represent a circular cross-section. The resulting discretization error at the borders prevents the numerical stiffness from quickly converging towards the analytical solution. Second, for a square cross-section the analytical solution is only an approximation. Therefore the numerical stiffness does not exactly converge towards the analytical solution.

**Bending.** For the bending problem we consider a Bernoulli-beam of a homogeneous material with a Poisson's ratio of $\nu = 0.3$. The size of the rectangular cross section is given by $L_y = 0.05L_x \times L_z = 0.08L_x$. As curvature we imposed $\kappa = 10/(9L_x)$.
<figure><img src="image_7.png" /><figcaption>FIG. 5. Chiral metamaterial. (a) Single unit cell of the metamaterial. The geometric parameters are indicated in the illustration. (b) Voxelized unit cell with 100x100x100 voxels. (c) Metamaterial consisting of 3x3x2 unit cells.</figcaption></figure>

The analytical bending stiffness for a rectangular cross-section is

$$ S_{\text{ana}} = EI_y \\ \text{with } I_y = \int_{0}^{L_y} dy \int_{0}^{L_z} dz z^2 = \frac{1}{3} L_z^3 L_y, \qquad (19) $$

where E is the Young's modulus and $L_y$ and $L_z$ are the length of the cross-section in y-direction and z-direction respectively. The numerical bending stiffness was calculated from the numerical stresses via the average numerical moment:

$$ S_{\text{num}} = \frac{M_{\text{num}}}{\kappa} = \frac{1}{\kappa} \frac{1}{L_x} \int_{0}^{L_x} dX \int_{0}^{L_y} dY \int_{0}^{L_z} dZ \sigma_{xx}(Z - Z_0) \qquad (20) $$

Figure 4 (b) shows the difference between the analytical stiffnesses and the numerical stiffnesses for different discretizations. We observe that the error of the bending stiffness continuously converges towards zero. Even for relatively coarse discretizations (20x20x20 voxels), the difference is already less than 0.1%. The result for the bending stiffness directly confirms the validity of our method.

## B. Chiral metamaterial

To demonstrate the capabilities of our approach to simulating complex, non-periodic structures, we modeled a chiral metamaterial. As an example, we chose the chiral metamaterial introduced by Frenzel et al. in [12]. Figure 5 shows the geometry of this metamaterial.
<figure><img src="image_8.png" /><figcaption>FIG. 6. Results of mesh refinement for 1x1x1 unit cells. (a) Average force in z-direction resulting from COMSOL calculations. (b) Average force in z-direction resulting from µSpectre calculations.</figcaption></figure>

In the simulations, the length of one unit cell in each direction is $a = 0.5\text{mm}/N_{xy}$ where $N_x = N_y = N_{xy}$ is the number of unit cells forming the metamaterial in the $x$ and $y$-directions. This means we decrease the size of the unit cell as we increase the number of cells, keeping the overall size of the structure constant. Accordingly, the lengths of the simulated metamaterial are $L_x = L_y = 0.5\text{ mm}$ and $L_z = aN_z$. The beam thickness $t$ is $0.06a$, the outer radius for all circles $r$ is $0.4a$ and the angle $\alpha$ describing the inclination of the beams is $80^\circ$. The base material is homogeneous, linear elastic and isotropic and has a Young's modulus of $2.6 \text{ GPa}$ and a Poisson's ratio of $0.4$. To study the behavior of this metamaterial we considered a torsion around the central axis $\vec{z}_0$, so that $X_0 = L_x/2$ and $Y_0 = L_y/2$, similar to the torsion of the beam shown in Fig. 1.

Since no analytical solution exists for such complex geometries, we have validated the results of our method by comparing them to finite element simulations with COMSOL (version 5). For all COMSOL simulations we have used an unstructured mesh of first order tetrahedral elements. Like those of µSpectre, the COMSOL simulations are made with the assumption of small-strains and for the quasi-static case.

In COMSOL model, we have used Dirichlet boundary conditions on lower and upper surfaces. We fix the lower surface be setting a zero displacement, i.e. $\vec{u} = 0$ for $Z = 0$, and we rotated the upper surface be prescribing displacement of $\vec{u} = (-w(Y - L_y/2)L_z, w(X - L_x/2)L_z, 0)^T$ for $Z = L_z$. All other surfaces are modeled as free surfaces.

First, we have performed mesh refinement studies for both µSpectre and COMSOL. Fig-
<figure><img src="image_9.png" /><figcaption>FIG. 7. Distribution of the stress $\sigma_{zz}$ in one side of the chiral metamaterial for a twist angle per length of 0.05/mm. The deformation of the geometry was plotted with a scaling factor of 5. The undeformed geometry is shown with thin black lines.</figcaption></figure>

ure 6 represents the result for one unit cell. It shows that the COMSOL simulations converge relatively smoothly with the number of degrees of freedom. In contrast, the µSpectre simulations converge, but not monotonously. Similar to the torsion of the cylinder, this is likely due to the difficulties of representing the smooth curves in the geometry of the metamaterial with voxels.

All following COMSOL simulations were made with $8 \times 10^5$ degrees of freedom per unit cell, while all following µSpectre simulations were made with a discretization of 100x100x100 voxels per unit cell, corresponding to $3 \times 10^7$ degrees of freedom.

Figure 7 shows the distribution of the stress component $\sigma_{zz}$ in one side of the deformed metamaterial, both for µSpectre and COMSOL. Overall, the results of µSpectre and COMSOL look similar: The deformation seems practically identical and both have values of $\sigma_{zz}$ in the same order of magnitude. In addition, the maxima and minima of $\sigma_{zz}$ are mostly in the same locations. Only at the corners does the COMSOL simulation has some additional extrema we do not observe in the µSpectre simulation. Also, the absolute values of $\sigma_{zz}$ in the µSpectre simulation seem altogether smaller than those in the COMSOL simulation. Considering that on one hand, the COMSOL simulation differs in the boundary conditions from the µSpectre simulation, and on the other hand, the voxel discretization can only approximate the complex geometry of the metamaterial, these differences are well within the acceptable range.

The COMSOL simulations differ in their boundary conditions from the µSpectre simulation.
<figure><img src="image_10.png" /><figcaption>FIG. 8. Average force in z-direction of the twisted chiral metamaterial. (a) Force as a function of the number of unit cells in z-direction for a twist angle per length of $w = 0.05$/mm. (b) Force as a function of the twist angle per length $w$. The number of unit cells in z-direction is $N_z = 1$ for µSpectre and $N_z = 4$ for COMSOL.</figcaption></figure>

273 lations, where we have periodic boundary conditions in z-direction. To evaluate the importance of this difference, we have studied the influence of the number of unit-cells in z-direction $N_z$ on the simulation results. Figure 8a shows the average force in z-direction for an applied twist angle per length $w = 0.05$/mm for different $N_z$. It illustrates that the results of µSpectre are independent of $N_z$, which corresponds to the periodic boundary conditions in Z direction used in µSpectre. In contrast, the COMSOL results are highly dependent on $N_z$. However, for larger $N_z$, the result seem to converge towards the solutions of µSpectre and the difference is less than 5% between $N_z = 3$ and $N_z = 4$. We conclude that the COMSOL simulations for $N_z = 4$ are close enough to the periodic boundary condition to allow a comparison with the µSpectre simulations.

Figure 8b shows the average force in z-direction when the chiral metamaterial is deformed with different twist angles per length $w$. The first observation is that all simulations with the chiral metamaterial clearly show the coupling between the torsion around the z-axis and the force in z-direction which is the characteristic of chiral metamaterials. This is emphasized by the µSpectre simulation of a cylinder, which has no such coupling. Next, in all four cases the force is proportional to the rotation angle. This is due to the small-strain
289 assumption we have used. The stiffnesses calculated with µSpectre are ca. 0.15 N/m for
290 $N_x = N_y = 1$ and ca. 0.046 N/m for $N_x = N_y = 3$. This decrease in the stiffness with
291 increasing $N_x$ corresponds to the results from Frenzel et al. [12]. It is likely due to the
292 fact that neighboring unit cells hinder each other when rotating in the same direction. The
293 stiffnesses calculated with COMSOL are approximate 0.18 N/m for $N_x = N_y = 1$ and around
294 0.067 N/m for $N_x = N_y = 3$ and therefore slightly larger than the stiffness calculated by
295 µSpectre. However, they are in the same order of magnitude and show the same effect when
296 increasing the number of unit cells in x- and y-direction. Keeping in mind the differences
297 between the simulation set-ups, these differences in the resulting stiffnesses are acceptably
298 small.

## 299 IV. CONCLUSIONS

300 In this study, we present a simple approach for simulating torsion and bending within
301 the framework of fast Fourier transform (FFT)-based solvers incorporating compatibility
302 projection. This approach relies on breaking the periodic boundary conditions with void
303 regions and applying a non-uniform macroscopic strain to the representative volume element.
304 To validate the approach, we present three test cases with existing (semi-)analytical solution:
305 The torsion of a cylinder with circular cross-section, the torsion of a beam with square cross-
306 section, and the bending of beam with rectangular cross-section. In addition, we use this
307 method to study the behavior of a chiral metamaterial under torsion and validate the results
308 by comparing them to a classic finite element simulation.

309 Our approach unlocks efficient analysis of non-periodic metamaterial loading scenar-
310 ios, enabling rapid computational design exploration for applications ranging from impact-
311 damping structures to adaptive mechanical systems. The method’s straightforward imple-
312 mentation within existing FFT frameworks makes it accessible to the broader computational
313 mechanics community.

314 ACKNOWLEDGEMENTS

315 We are indebted to Lucia Brunold for help with Figs. 1 and 2. IJ and LP acknowledge
316 funding by the Carl Zeiss Foundation (Research cluster "Interactive and Programmable Ma-
terials - IPROM") and the Deutsche Forschungsgemeinschaft (EXC 2193/1 - 390951807).
ML acknowledges the funding by the European Commission (Marie Skłodowska-Curie Fel-
lowship 101106585 - microFFTTO). We also received infrastructure support by the state of
Baden-Württemberg through the bwHPC program and the German Research Foundation
(DFG) through grant no INST 39/963-1 FUGG (bwForCluster NEMO).

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

[1] L. Montemayor, V. Chernow, and J. R. Greer, Materials by design: Using architecture in material design to reach new property spaces, MRS Bull. **40**, 1122 (2015).

[2] M. Kadic, G. W. Milton, M. van Hecke, and M. Wegener, 3D metamaterials, Nat. Rev. Phys. **1**, 198 (2019).

[3] Z. Jia, F. Liu, X. Jiang, and L. Wang, Engineering lattice metamaterials for extreme property, programmability, and multifunctionality, J. Appl. Phys. **127** (2020).

[4] X. Zheng, H. Lee, T. H. Weisgraber, M. Shusteff, J. DeOtte, E. B. Duoss, J. D. Kuntz, M. M. Biener, Q. Ge, J. A. Jackson, et al., Ultralight, ultrastiff mechanical metamaterials, Science **344**, 1373 (2014).

[5] R. Lakes, Foam structures with a negative Poisson's ratio, Science **235**, 1038 (1987).

[6] M. J. Taylor, L. Francesconi, M. Gerendás, A. Shanian, C. Carson, and K. Bertoldi, Low porosity metallic periodic structures with negative Poisson's ratio, Adv. Mater. (2014).

[7] X. Ren, R. Das, P. Tran, T. D. Ngo, and Y. M. Xie, Auxetic metamaterials and structures: a review, Smart Mater. Struct. **27**, 023001 (2018).

[8] G. Felsch, N. Ghavidelnia, D. Schwarz, and V. Slesarenko, Controlling auxeticity in curved-beam metamaterials via a deep generative model, Comput. Methods Appl. Mech. Eng. **410**, 116032 (2023).

[9] R. Gatt and J. N. Grima, Negative compressibility, Phys. Status Solidi RRL **2**, 236 (2008).

[10] C. S. Ha, M. E. Plesha, and R. S. Lakes, Chiral three-dimensional isotropic lattices with negative Poisson's ratio, Phys. Status Solidi B **253**, 1243 (2016).

[11] I. Fernandez-Corbaton, C. Rockstuhl, P. Ziemke, P. Gumbsch, A. Albiez, R. Schwaiger, T. Frenzel, M. Kadic, and M. Wegener, New twists of 3D chiral metamaterials, Adv. Mater. **31**, 1807742 (2019).
[12] T. Frenzel, M. Kadic, and M. Wegener, Three-dimensional mechanical metamaterials with a twist, Science **358**, 1072 (2017).

[13] M. Fu, F. Liu, and L. Hu, A novel category of 3D chiral material with negative Poisson's ratio, Compos. Sci. Technol. **160**, 111 (2018).

[14] F. Roters, P. Eisenlohr, C. Kords, D. D. Tjahjanto, M. Diehl, and D. Raabe, DAMASK: the Düsseldorf Advanced MAterial Simulation Kit for studying crystal plasticity using an FE based or a spectral numerical solver, Procedia IUTAM **3**, 3 (2012).

[15] H. Moulinec and P. Suquet, A fast numerical method for computing the linear and nonlinear mechanical properties of composites, C. R. Acad. Sci. Ser. II **318**, 1417 (1994).

[16] H. Moulinec and P. Suquet, A numerical method for computing the overall response of non-linear composites with complex microstructure, Comput. Meth. Appl. Mech. Eng. **157**, 69 (1998).

[17] M. Schneider, A review of nonlinear FFT-based computational homogenization methods, Acta Mech. **323**, 2051 (2021).

[18] S. Lucarini, M. V. Upadhyay, and J. Segurado, FFT based approaches in micromechanics: Fundamentals, methods and applications, Model. Simul. Mater. Sci. Eng. **30**, 023002 (2021).

[19] C. Gierden, J. Kochmann, J. Waimann, B. Svendsen, and S. Reese, A review of FE-FFT-based two-scale methods for computational modeling of microstructure evolution and macroscopic material behavior, Arch. Comput. Methods Eng. **29**, 4115 (2022).

[20] T. W. J. de Geus, J. E. P. van Duuren, R. H. J. Peerlings, and M. G. D. Geers, Fracture initiation in multi-phase materials: A statistical characterization of microstructural damage sites, Materials Science and Engineering: A **673**, 551 (2016).

[21] S. Lucarini and J. Segurado, On the accuracy of spectral solvers for micromechanics based fatigue modeling, Comput. Mech. **63**, 365 (2019).

[22] M. Frigo and S. G. Johnson, The design and implementation of FFTW3, Proc. IEEE **93**, 216 (2005).

[23] J. Vondřejc, J. Zeman, and I. Marek, An FFT-based Galerkin method for homogenization of periodic media, Comput. Math. Appl. **68**, 156 (2014).

[24] J. Zeman, T. W. J. de Geus, J. Vondřejc, R. H. J. Peerlings, and M. G. D. Geers, A finite element perspective on nonlinear FFT-based micromechanical simulations, Int. J. Numer. Methods Eng. **111**, 903 (2017).
[25] T. W. J. de Geus, J. Vondřejc, J. Zeman, R. H. J. Peerlings, and M. G. D. Geers, Finite strain FFT-based non-linear solvers made simple, Comput. Methods Appl. Mech. Eng. 318, 412 (2017).

[26] M. Ladecký, R. J. Leute, A. Falsafi, I. Pultarová, L. Pastewka, T. Junge, and J. Zeman, An optimal preconditioned FFT-accelerated finite element solver for homogenization, Appl. Math. Comput. 446, 127835 (2023).

[27] L. Gélébart, A modified FFT-based solver for the mechanical simulation of heterogeneous materials with Dirichlet boundary conditions, Comptes Rendus. Mécanique 348, 693 (2020).

[28] L. Risthaus and M. Schneider, Imposing Dirichlet boundary conditions directly for FFT-based computational micromechanics, Comput. Mech. 74, 1089 (2024).

[29] L. Gélébart, FFT-based simulations of heterogeneous conducting materials with combined non-uniform neumann, periodic and dirichlet boundary conditions, Eur. J. Mech. A/Solids 105, 105248 (2024).

[30] J. Paux, L. Morin, and L. Gélébart, A discrete sine-cosine transforms Galerkin method for the conductivity of heterogeneous materials with mixed Dirichlet/Neumann boundary conditions, Int. J. Numer. Methods Eng. 126, e7615 (2025).

[31] M. Zecevic and R. A. Lebensohn, Extended FFT-based micromechanical formulation to consider general non-periodic boundary conditions, Int. J. Solids Struct. 311, 113225 (2025).

[32] R. J. Leute, M. Ladecký, A. Falsafi, I. Jödicke, I. Pultarová, J. Zeman, T. Junge, and L. Pastewka, Elimination of ringing artifacts by finite-element projection in FFT-based homogenization, J. Comput. Phys. 453, 110931 (2022).

[33] D. Gottlieb and C.-W. Shu, On the Gibbs phenomenon and its resolution, SIAM Rev. 39, 644 (1997).

[34] M. Schneider, F. Ospald, and M. Kabel, Computational homogenization of elasticity on a staggered grid, Int. J. Numer. Methods Eng. 105, 693 (2016).

[35] X. Ma, M. Shakoor, D. Vasiukov, S. V. Lomov, and C. H. Park, Numerical artifacts of fast fourier transform solvers for elastic problems of multi-phase materials: their causes and reduction methods, Comput. Mech. 67, 1661 (2021).

[36] T. Maiti and P. Eisenlohr, Fourier-based spectral method solution to finite strain crystal plasticity with free surfaces, Scr. Mater. 145, 37 (2018).
[37] Y. Chen, D. Vasiukov, L. Gélébart, and C. H. Park, A FFT solver for variational phase-field modeling of brittle fracture, Comput. Meth. Appl. Mech. Eng. **349**, 167 (2019).

[38] T. Dumitrică and R. James, Objective molecular dynamics, J. Mech. Phys. Solids **55**, 2206 (2007).

[39] D.-B. Zhang, M. Hua, and T. Dumitrică, Stability of polycrystalline and wurtzite Si nanowires via symmetry-adapted tight-binding objective molecular dynamics, J. Chem. Phys. **128**, 084104 (2008).

[40] P. Koskinen and O. O. Kit, Efficient approach for simulating distorted materials, Phys. Rev. Lett. **105**, 106401 (2010).

[41] O. O. Kit, L. Pastewka, and P. Koskinen, Revised periodic boundary conditions: Fundamentals, electrostatics, and the tight-binding approximation, Phys. Rev. B **84**, 155431 (2011).

[42] https://gitlab.com/muspectre/muspectre.

[43] K.-H. Grote and J. Feldhusen, *DUBBEL: Taschenbuch für den Maschinenbau* (Springer, 2005).