# FFT-based surrogate modeling of auxetic metamaterials with real-time prediction of effective elastic properties and swift inverse design

Hooman Danesh, Daniele Di Lorenzo, Francisco Chinesta, Stefanie Reese, Tim Brepols

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

To cite this version:

Hooman Danesh, Daniele Di Lorenzo, Francisco Chinesta, Stefanie Reese, Tim Brepols. FFT-based surrogate modeling of auxetic metamaterials with real-time prediction of effective elastic properties and swift inverse design. Materials & Design, 2024, 248, pp.113491. 10.1016/j.matdes.2024.113491. hal-05133379

HAL Id: hal-05133379

https://cnam.hal.science/hal-05133379v1

Submitted on 27 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
<figure><img src="image_3.png" /></figure>

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

# FFT-based surrogate modeling of auxetic metamaterials with real-time prediction of effective elastic properties and swift inverse design

Hooman Danesh<sup>a,*</sup>, Daniele Di Lorenzo<sup>b,c</sup>, Francisco Chinesta<sup>b,c,d</sup>, Stefanie Reese<sup>a,e</sup>, Tim Brepols<sup>a</sup>

## ARTICLE INFO

### Keywords:

Auxetic structures

Surrogate models

Effective properties

Inverse analysis

FFT-based homogenization

## ABSTRACT

Auxetic structures, known for their negative Poisson's ratio, exhibit effective elastic properties heavily influenced by their underlying geometry and base material properties. While periodic homogenization of auxetic unit cells can be used to investigate these properties, it is computationally expensive and limits design space exploration and inverse analysis. In this paper, the fast Fourier transform (FFT)-based homogenization approach is adopted to efficiently generate data for developing surrogate models, bypassing concerns about periodic mesh generation and boundary conditions typically associated with the finite element method (FEM). Surrogate models are developed for the real-time prediction of the effective elastic properties of auxetic unit cells with orthogonal voids of different shapes. The generated surrogate models accept geometric parameters and base material properties as inputs to predict the effective elastic constants in real-time. This rapid evaluation enables a practical inverse analysis framework for obtaining the optimal design parameters that yield the desired effective response. The performance of the generated surrogate models is rigorously examined through a train/test split methodology, a parametric study, and an inverse problem. Finally, a graphical user interface (GUI) is developed, offering real-time prediction of the effective tangent stiffness and performing inverse analysis to determine optimal geometric parameters.

## 1. Introduction

### 1.1. Context and rationale

Auxetic materials represent a fascinating frontier in materials science and engineering, distinguished by their unique and counter-intuitive property of exhibiting negative Poisson's ratios. Unlike conventional materials that contract laterally when stretched, auxetic materials exhibit the remarkable property of expanding in the direction transverse to the applied tension. This counter-intuitive feature has a wide range of innovative applications, including energy absorption [1–3], flexible electronics [4–6], protective sports equipment [7–9], and medical implants [10–12]. Such an extraordinary behavior, i.e., negative Poisson's ratio, stems from the unique structural geometry of auxetic materials. The effective response of the designed auxetic structures is substantially impacted by changing the base material properties and the unit

cell geometry. Each engineering application demands a unique set of auxetic structures that possess certain effective features. For this reason, it would be desirable to establish a framework that can extract the effective properties of auxetics with varying geometries and base constituents.

Extensive research has been conducted on the homogenization and extraction of the effective properties of auxetic metamaterials using various analytical, experimental, and numerical approaches. An analytical model to predict the effective elastic constants and Poisson's ratio of hexagonal and re-entrant honeycomb structures was developed in [13]. Using a discrete asymptotic homogenization approach, Dos Reis and Ganghoffer [14] provided closed form expressions for the effective elastic moduli and Poisson's ratio of re-entrant and chiral structures formed from a network of beams. A theoretical model for beam networks of 3D re-entrant honeycomb auxetic structures was developed in

<footnote><sup>a</sup> Institute of Applied Mechanics, RWTH Aachen University, Mies-van-der-Rohe-Str. 1, 52074 Aachen, Germany</footnote>

<footnote><sup>b</sup> ESI Group, Symbiose 2, 10 Av. Aristide Briand, 92220 Bagneux, France</footnote>

<footnote><sup>c</sup> PIMM Lab, ENSAM Institute of Technology, 151 Boulevard de l'Hôpital, 75013 Paris, France</footnote>

<footnote><sup>d</sup> CNRS @CREATE, 1 CREATE Way, 04-05 CREATE Tower, 138602, Singapore</footnote>

<footnote><sup>e</sup> University of Siegen, Adolf-Reichwein-Str. 2a, 57076 Siegen, Germany</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: hooman.danesh@rwth-aachen.de (H. Danesh).</footnote>
[15], providing analytical solutions for the elastic modulus, Poisson's ratio, and yield strength. This model was validated through comparison with experimental and finite element (FE) simulation results. Chan and Evans [16,17] carried out experimental studies to investigate the elastic response of auxetic polymeric foams under tension, compression, and shear. Combined FE simulations and experimental testing investigations of the mechanical properties of re-entrant cell honeycombs under uniaxial loading were performed in [18]. Dirrenberger et al. [19] employed the finite element method (FEM) combined with periodic homogenization to extract the effective elastic moduli of chiral and honeycomb lattices and to examine their anisotropy. The homogenized mechanical properties of hyperelastic inclusion-matrix composites were obtained via numerical homogenization by Kochmann and Venturini [20]. Their study revealed that auxeticity in such finitely strained composites emerged for certain geometries and base material properties. Slann et al. [21] developed a numerical asymptotic homogenization FE approach to conduct a parametric study on the mechanical properties of rectangular void auxetic unit cells. The model's accuracy was verified by comparing its predictions with experimental results obtained from uniaxial tensile tests on perforated sheets. Numerical and experimental investigations of auxetic structures with randomly oriented perforations were performed by Grima et al. [22]. They employed the FE approach and periodic boundary conditions to extract the average response of these auxetics.

With regard to the previously mentioned analytical, numerical, and experimental methods used to study the effective response of auxetic metamaterials, it is crucial to mention their limitations in specific situations. Analytical approaches are mainly limited to the extraction of the effective properties for simple discrete structures, such as those formed from a network of beams. Furthermore, when it comes to design space exploration, real-time prediction of effective properties, or performing inverse design, numerical simulations or experimental studies present significant challenges. Numerical simulations, based on either FE or fast Fourier transform (FFT) methods, are computationally expensive, particularly for inverse problems. Additionally, experimental methods are impractical and often too costly for iterative design processes. Consequently, there is a critical need for surrogate models that can provide real-time evaluations of effective properties for auxetic structures. These models, which are trained on high-fidelity datasets, enable rapid design space exploration, sensitivity analyses, and inverse design of auxetic structures.

Random forests regression was used by Tajalasir et al. [23] to develop machine learning models for predicting the dynamic impact stress of hierarchical auxetic structures. In a related study, Ben-Yelun et al. [24] employed the random forests method to create surrogate models for estimating the effective elastic modulus and Poisson's ratio of a class of tunable 3D auxetic metamaterials, using geometric parameters as inputs. Neural networks were used in [25] to develop machine learning models, trained on large datasets produced by FE simulations, to predict Poisson's ratio of cementitious oval void axuetics. In another study [26], machine learning models were developed based on experimental data to predict the strength, stiffness, and energy absorption of re-entrant auxetic structures manufactured by fused deposition modeling. A machine learning model to predict the effective Poisson's ratio of oval void auxetic structures was developed in [27]. The accuracy of the model was examined against experiments and FE simulations. Du et al. [28] employed recurrent neural networks to predict the effective Poisson's ratio of kirigami auxetics based on four geometric input parameters. A neural-network-based surrogate model, trained on datasets from isogeometric analysis homogenization, was developed in [29] to establish an efficient inverse design framework for designing tetra-chiral auxetics with target Poisson's ratios. Chang et al. [30] employed back-propagation neural networks and the genetic algorithm to develop a model for the inverse design of re-entrant auxetics with zero Poisson's ratio. Similarly, the same machine learning algorithms were used in [31] to set up an efficient data-driven framework for the optimization design of auxetic

structures with peanut-shaped perforations. Bastek et al. [32] utilized a deep-learning framework combining neural networks with physical constraints for the inverse design of truss lattices, achieving a fully tailorable 3D anisotropic stiffness for non-auxetic structures. A neural-network-based inverse design approach was developed in [33] to achieve specified stress-strain and/or Poisson's ratio-strain curves for a novel semi-auxetic mechanical metamaterial, validated through FE simulations and experimental tests on 3D-printed specimens. In another study, Peng and Xu [34] developed neural network surrogate models for both forward prediction and inverse design of bi-material composite triangular lattice structures, encompassing a range of tunable effective elastic properties, including non-auxetic and auxetic cases. Machine-learning-based surrogate models for metamaterials have been extensively developed in recent literature, which are not within the scope of this paper. For a comprehensive overview of these advancements, interested readers are referred to in-depth reviews conducted in this field. For instance, Zheng et al. [35] reviewed the use of deep learning for property prediction, geometry generation, and inverse design, emphasizing its potential to replace traditional design methods with data-driven approaches. Similarly, Lee et al. [36] provided a thorough overview of data-driven design methodologies, focusing on unit cell design and multiscale optimization across various physical domains. Lastly, a more focused review on machine-learning-based surrogate models specifically for auxetic structures can be found in [37].

## 1.2. *In silico* data generation

Efficient data generation is crucial in the development of accurate and robust surrogate models, including those for auxetic metamaterials. While FEM is commonly employed for generating synthetic data in the development of surrogate models, an alternative approach exists that has been largely overlooked in this context. The FFT-based method offers a highly efficient and rapid means of data generation, presenting significant potential for enhancing the process of building surrogate models for auxetic metamaterials. Owing to their matrix-free formulation, FFT-based methods can be pretty memory-efficient, allowing for rather large systems to be homogenized in comparison to FEM. In FFT-based homogenization, periodicity follows naturally from the definition of the Fourier transform and does not need extra care, as in FEM. This intrinsic periodicity eliminates the need to check for symmetric node placement on boundaries, which is crucial in FEM for applying periodic boundary conditions.

Regarding the accuracy and computational efficiency of FFT-based methods compared to FEM, Lucarini et al. [38] employed the FFT-based homogenization for metamaterial unit cells and conducted a comprehensive comparison with FEM simulations. They demonstrated that FFT becomes competitive with FEM for relative densities greater than 7% when using the same discretization within the unit cell interior. For instance, for relative densities around 30% (70% porosity), FFT simulations were found to be 4 to 8 times faster than FEM. Moreover, when evaluating the effective elastic response, the FFT results converged to values similar to those obtained by using FEM, with the maximum difference remaining below 0.6%.

The FFT-based homogenization dates back to the primary works of Moulinec and Suquet [39,40] on the effective mechanical response of composites. The FFT-based scheme is formulated in a convolution integral form, which is known as the Lippmann-Schwinger equation, and is then transformed to a more manageable algebraic product form in Fourier space (see, e.g., [41,42] for review). In spite of its well-known advantages, the basic FFT scheme of Moulinec and Suquet [39,40] reaches a converged solution notably dependent on the difference in stiffness between the phases found in the microstructure. This means that the method converges faster when the stiffness contrast between the phases is smaller. In the context of metamaterial unit cells, including the auxetic structures in the present study, one of these phases symbolizes an empty space possessing zero stiffness. Consequently, an infinite
stiffness contrast arises, impeding the convergence of the basic scheme when applied to such materials [38,43]. To address this limitation, Lucarini et al. [38] proposed a solution involving the utilization of the Galerkin FFT homogenization method [44,45], along with the rotated finite difference grid [46] and the minimal residual (MINRES) solver. Such an approach is employed in the present work to generate large high-fidelity datasets by solving the homogenization problem for different auxetic structures.

## 1.3. Objectives and contributions

Although FFT-based homogenization approaches are relatively straightforward to implement and use, the computational cost becomes prohibitive when it comes to real-time parametric studies and inverse design. Within this context, parametric studies enable engineers to rapidly evaluate how changes in geometric and base material parameters affect the effective elastic properties of auxetic structures, facilitating quick sensitivity analyses and understanding parameter-property relationships. Equally important is inverse design, where desired material properties, such as target values for the effective Poisson's ratio or other elastic constants, are defined, and the optimal structures meeting these properties are found using optimization algorithms. All these tasks require numerous evaluations of the homogenization problem for different microstructures. Therefore, the development of fast and accurate surrogate models for real-time prediction of effective properties is crucial. These surrogates enable efficient design space exploration, rapid parametric studies, evaluation of potential unit cell configurations, and swift solution of inverse problems without resorting to time-consuming and computationally expensive simulations.

The present work contributes to the field of auxetic metamaterial design and analysis through:

* **Employing FFT-based homogenization for efficient data generation:** The efficient FFT-based homogenization method is used to generate large, high-fidelity datasets required to build the surrogate models. This approach represents a novel and advantageous method for data generation in this context. The exceptional efficiency of this method allows for the creation of substantial datasets within a relatively short period of time. Consequently, the accuracy of the surrogate models can be significantly enhanced due to the voluminous data provided, marking a notable improvement over traditional data generation techniques in this field.

* **Developing accurate surrogate models for real-time prediction:** Surrogate models are developed for four types of orthogonal void auxetic unit cells with different perforation shapes, including rectangular, diamond, oval, and peanut-shaped voids. The random forests regression method is employed to construct these surrogates, which accept geometric and base material parameters as input and predict the homogenized effective tangent stiffness in real time. This new real-time prediction capability facilitates effective probe of design spaces, allowing for rapid parametric studies and sensitivity analyses. The performance of the models is examined through a train/test split approach and a parametric study. It is demonstrated that the generated surrogate models can provide highly accurate predictions in close agreement with the actual values obtained from the FFT-based solver.

* **Establishing a swift inverse design framework:** An inverse design framework is also established using an iterative grid refinement approach. While many previous studies have primarily focused on the effective Poisson's ratio as the target value for the inverse problem, the current approach offers the novel flexibility to select any component or any combination of the components of the effective elastic stiffness as the target output. The capability of the established framework in swiftly solving the inverse problem is shown for the rectangular void unit cell.

* **Developing a GUI for practical application:** A dedicated graphical user interface (GUI) is developed to enhance the accessibility and practical application of the developed surrogate models. This interface enables the real-time prediction of effective elastic constants and facilitates the inverse design of auxetic metamaterials, providing an efficient tool for both analyses and design processes with visual feedback.

The remainder of the present paper is structured as follows: Section 2 presents the theoretical basis of the FFT-based homogenization and, specifically, the computation of the homogenized effective tangent stiffness. Section 3 is concerned with the surrogate modeling of the given auxetic structures, introducing the input parameters and output quantities of the models as well as the random forests regression method. This section also includes a parametric study to assess the accuracy of the generated models, to investigate the auxetic behavior of various unit cell geometries, and to examine the effect of dataset size. The inverse analysis framework, utilizing an iterative grid refinement approach, is explained in Section 4, followed by a practical inverse problem to achieve the geometric parameters for specific mechanical properties. Section 5 introduces the developed GUI, highlighting its user-friendly capabilities in real-time prediction of effective stiffness and inverse analyses. Finally, Section 6 summarizes concluding remarks and future research directions.

## 2. Galerkin FFT homogenization

The FFT-based homogenization was first introduced by Moulinec and Suquet for the homogenization of linear and nonlinear composites [39,40]. They reformulated the homogenization problem into a periodic Lippmann-Schwinger equation, expressed as a convolution of the local strain fluctuations arising from elastic heterogeneities with the Green's function of an elastic reference medium. With the Green's function available in closed form in Fourier space, they leveraged FFT to handle the convolution as an algebraic product in Fourier space, solving the it by a fixed-point scheme. By operating naturally on regular grids, the FFT-based homogenization enables seamless incorporation of microstructures using pixels or voxels from digital images.

Nevertheless, the convergence of the Moulinec-Suquet basic scheme is highly dependent on the stiffness contrast between phases in the microstructure, resulting in non-convergence for metamaterials with one phase representing voids. Since the development of FFT-based homogenization techniques, a number of approaches have been proposed to address high or infinite stiffness contrast [47–50]. Among these, Lucarini et al. [38] successfully employed an enhanced approach, incorporating the Galerkin FFT method [44,45] in combination with the rotated finite difference grid [46] and the MINRES solver, which proved effective for metamaterials with infinite stiffness contrast. The same approach is employed in the present work to generate high-fidelity datasets for orthogonal void auxetic metamaterials.

Following [38,44,45], the initial step in the Galerkin FFT method is to apply the principle of virtual work to formulate the weak form of the boundary value problem. Considering small strain kinematics, the weak form, in the absence of boundary traction terms due to periodic boundary conditions, is written as

$$ \int_{\Omega} \delta\epsilon(\mathbf{x}) : \sigma(\mathbf{x}) d\Omega = 0, \qquad (1) $$

where $\Omega$ represents the microscopic domain, and $\sigma(\mathbf{x})$ is the microscopic stress tensor at the microscopic position $\mathbf{x}$. The microscopic virtual strain $\delta\epsilon(\mathbf{x})$ needs to satisfy the periodicity, symmetry, and compatibility conditions. The last two conditions can be applied by employing a projection operator $\mathbb{G}$ (for more details, cf. [44]) as

$$ \delta\epsilon(\mathbf{x}) = (\mathbb{G} * \zeta)(\mathbf{x}). \qquad (2) $$
In the above, $\zeta(\mathbf{x})$ is an arbitrary second-order tensor field, whose symmetry and compatibility are imposed by the fourth-order projection tensor $\mathbb{G}$ through the convolution $*$ in real space. By inserting Eq. (2) into Eq. (1) and considering the major symmetry of $\mathbb{G}$ (i.e., $G_{ijkl} = G_{klij}$), the following integral equation emerges:

$$ \int_{\Omega} \zeta(\mathbf{x}) : [(\mathbb{G} * \sigma)(\mathbf{x})] d\Omega = 0. \quad (3) $$

As Eq. (3) needs to be fulfilled for any arbitrary test function $\zeta(\mathbf{x})$, the equilibrium in its weak form reduces to the following equation, which is valid at every point in domain $\Omega$:

$$ \mathbb{G} * \sigma = \mathbf{0}. \quad (4) $$

Performing the forward ($\mathcal{F}$) and backward ($\mathcal{F}^{-1}$) Fourier transforms and taking the convolution property of Fourier transform [42] into account, Eq. (4) can be written as

$$ \mathbb{G} * \sigma = \mathcal{F}^{-1} [\hat{\mathbb{G}} : \mathcal{F}(\sigma)] = \mathbf{0}, \quad (5) $$

where $\hat{\mathbb{G}}$ is the projection tensor in the Fourier space, whose components are expressed as follows [38,44,45]:

$$ \hat{\mathbb{G}}_{ijkl} = \begin{cases} 0, & \text{for zero and Nyquist frequencies} \\ \left[ I_{pkq}^s \xi_p \xi_q \right]^{-1} \xi_j \xi_l, & \text{for other frequencies} \end{cases} \quad (6) $$

In Eq. (6), $I^s$ represents the fourth-order symmetric identity tensor, and $\xi$ is the frequency vector, which can be written as

$$ \begin{aligned} \xi_m &= iq_m N_m / L_m, & \text{with} \quad q_m &= 2\pi (n_m - N_m/2) / N_m, \\ n_m &= 0, 1, \dots, N_m - 1. & \end{aligned} \quad (7) $$

Here, $m = 1, 2$ and $m = 1, 2, 3$ for 2D and 3D domains, respectively, $i$ is reserved for the unit imaginary number, $L_m$ represents the length of the unit cell edge in the $m$ direction, and $N_m$ expresses the number of voxels/pixels in this direction.

As suggested in [38], in order to reach better convergence for materials with high stiffness contrast (e.g., the auxetic structures in the current study), Willot's rotated finite difference discretization [46] is used. In this case, the spatial frequency vector $\xi$ is modified as

$$ \xi'_m = i2N_m/L_m \tan(q_m/2) \prod_{m=1}^{d} \frac{1}{2} (1 + e^{-iq_m}). \quad (8) $$

Next, assuming linear elastic constitutive behavior and using the additive decomposition of the microscopic strain field $\varepsilon(\mathbf{x})$ into macroscopic strain $\bar{\varepsilon}$ and fluctuation strain $\bar{\varepsilon}(\mathbf{x})$, the stress-strain constitutive relation can be written as

$$ \sigma(\mathbf{x}) = \mathbb{D} : [\bar{\varepsilon} + \bar{\varepsilon}(\mathbf{x})], \quad (9) $$

where $\mathbb{D}$ is the fourth-order elastic stiffness tensor, and the macroscopic strain field $\bar{\varepsilon}$ is defined by volume averaging through

$$ \bar{\varepsilon} = \langle \varepsilon \rangle = \frac{1}{V} \int_V \varepsilon(\mathbf{x}) dV. \quad (10) $$

Here, $\langle \cdot \rangle$ denotes the average over the volume. Substituting Eq. (9) into Eq. (5), the equilibrium equation can be rewritten in the form

$$ \mathcal{F}^{-1}[\hat{\mathbb{G}} : \mathcal{F}(\mathbb{D} : \bar{\varepsilon})] = -\mathcal{F}^{-1}[\hat{\mathbb{G}} : \mathcal{F}(\mathbb{D} : \bar{\varepsilon})]. \quad (11) $$

Finally, Eq. (11) can be solved using any Newton-Krylov solver; however, the MINRES solver is employed here for auxetic structures with infinite stiffness contrast [38]. After computing the microscopic stress field $\sigma(\mathbf{x})$, the macroscopic average stress $\bar{\sigma}$ is obtained by

$$ \bar{\sigma} = \langle \sigma \rangle = \frac{1}{V} \int_V \sigma(\mathbf{x}) dV. \quad (12) $$

## 2.1. Homogenized effective tangent stiffness

The formulations presented in Section 2 can be used to compute the full-field solutions for the microscopic strain field $\varepsilon(\mathbf{x})$ and the macroscopic stress field $\sigma(\mathbf{x})$ as well as their corresponding macroscopic average values $\bar{\varepsilon}$ and $\bar{\sigma}$. However, in the current study, the homogenized effective tangent stiffness $\mathbb{C}$ is the main quantity of interest to compute. This homogenized effective tangent stiffness $\mathbb{C}$ is defined by

$$ \mathbb{C} = \frac{\partial \bar{\sigma}}{\partial \bar{\varepsilon}}. \quad (13) $$

When linear elasticity is concerned, a straightforward way to determine the homogenized effective stiffness tensor $\mathbb{C}$ is through a simple perturbation method [51,52]. This involves calculating the macroscopic average stress $\bar{\sigma}_{ij}$ induced by six (in 3D) or three (in 2D) macroscopic strain loadings $\bar{\varepsilon}_{kl} = \beta$. Here, $\beta$ is a real constant, which can be simply set to $\beta = 1$ since its value is irrelevant within linear elasticity. The unperturbed state is considered to be the state of zero strain. Consequently, the components of the homogenized effective tangent stiffness $C_{ijkl}$ can be computed by

$$ C_{ijkl} = \frac{\bar{\sigma}_{ij(kl)}}{\beta}, \quad (14) $$

where $\bar{\sigma}_{ij(kl)}$ refers to the macroscopic stress component $ij$ resulting from the application of a perturbation in the macroscopic strain component $kl$ while all other strain components remain zero.

An alternative approach to computing the homogenized effective stiffness tensor $\mathbb{C}$ is to directly solve for it using the projection operator $\mathbb{G}$ [53–55]. To this end, keeping Eqs. (10) and (12) in mind and taking advantage of the additive split of strain, i.e., $\varepsilon(\mathbf{x}) = \bar{\varepsilon} + \bar{\varepsilon}(\mathbf{x})$, Eq. (13) can be rewritten as

$$ \mathbb{C} = \left\langle \frac{\partial \sigma}{\partial \varepsilon} : (\mathbb{I} + \frac{\partial \bar{\varepsilon}}{\partial \bar{\varepsilon}}) \right\rangle = \frac{1}{V} \int_V \left[ \frac{\partial \sigma}{\partial \varepsilon} : (\mathbb{I} + \frac{\partial \bar{\varepsilon}}{\partial \bar{\varepsilon}}) \right] dV. \quad (15) $$

Given that $\partial\sigma/\partial\varepsilon = \mathbb{D}$, the only term in Eq. (15) that needs to be computed to achieve the homogenized effective stiffness $\mathbb{C}$ is the sensitivity of the strain fluctuation field with respect to the strain average field, i.e., $\partial\bar{\varepsilon}/\partial\bar{\varepsilon}$. This term will appear by substituting Eq. (9) into Eq. (4) and, subsequently, differentiating with respect to the macroscopic strain field $\bar{\varepsilon}$:

$$ \mathbb{G} * [\mathbb{D} : (\mathbb{I} + \frac{\partial \bar{\varepsilon}}{\partial \bar{\varepsilon}})] = \mathbf{0}. \quad (16) $$

Similar to the equilibrium equation, i.e., Eqs. (4), (5) and (11), using the forward and inverse Fourier transforms, Eq. (16) can also be rewritten as

$$ \mathcal{F}^{-1}[\hat{\mathbb{G}} : \mathcal{F}(\mathbb{D} : \frac{\partial \bar{\varepsilon}}{\partial \bar{\varepsilon}})] = -\mathcal{F}^{-1}[\hat{\mathbb{G}} : \mathcal{F}(\mathbb{D})]. \quad (17) $$

Eventually, analogous to the solution method employed for Eq. (11), Eq. (17) is also solved using the MINRES Newton-Krylov solver, which is particularly effective for auxetic structures with infinite stiffness contrast. In the present study, this second approach, i.e., adopting the projection operator $\mathbb{G}$ and Eqs. (15) and (17), is used to compute the homogenized effective tangent stiffness $\mathbb{C}$.

## 3. Surrogate modeling

This section is devoted to the surrogate modeling of four different auxetic structures shown in Fig. 1. These auxetic structures are formed from orthogonal voids of different shapes, including rectangular, diamond, oval, and peanut-shaped voids. The solid regions are assumed to be composed of isotropic material characterized by Young's modulus $E$.
<figure><img src="image_5.png" /><figcaption>Fig. 1. Auxetic unit cells with orthogonal voids of different shapes, characterized by the unit cell length L and void diameters d and D.</figcaption></figure>

and Poisson's ratio $\nu$. The input parameters and the output quantities for the desired surrogate models are discussed in the following sections.

### 3.1. Input parameters

As shown in Fig. 1, the geometric input parameters for each of the four structures are the relative void diameters $d/L$ and $D/L$. The geometric descriptions of the rectangular, diamond, and oval void unit cells are relatively straightforward. However, the peanut-shaped void requires a more complex mathematical representation. For this particular void geometry, the following mathematical expression in the polar coordinate system is used [56]:

$$ \rho(\theta) = \sqrt{\cos(2\theta) + \sqrt{1 - \sin^2(2\theta)}}, \qquad (18) $$

which can be transferred to the Cartesian coordinate system and scaled with the relative void diameters $d/L$ and $D/L$ as

$$ \begin{aligned} x_1 &= \frac{d}{L} \frac{1}{\sqrt{1 + \sqrt{1.1}}} \rho(\theta) \cos \theta, \\ x_2 &= \frac{D}{L} \frac{1}{\sqrt{1 + \sqrt{1.1}}} \rho(\theta) \sin \theta, \quad 0 \le \theta \le 2\pi. \end{aligned} \qquad (19) $$

In the above, $\rho(\theta)$ represents the radial distance from the origin to a point on the peanut-shaped curve as a function of the polar angle $\theta$. The variables $x_1$ and $x_2$ denote the Cartesian coordinates of points on the scaled peanut-shaped void. Moreover, the factor $\frac{1}{\sqrt{1+\sqrt{1.1}}} $ is a normalization constant to ensure proper scaling.

Regarding the material input parameters, by examining Hooke's law, it is evident that the microscopic elastic tensor $D$ is linear in terms of Young's modulus $E$ and nonlinear in terms of Poisson's ratio $\nu$. Consequently, Young's modulus functions as a scaling parameter. This property allows for the simplification of the surrogate models by considering only the effect of Poisson's ratio $\nu$ as the material input parameter, as the effect of Young's modulus can be accounted for through simple scaling of the results subsequent to the prediction.

### 3.2. Output quantities

Considering linear elastic behavior, the constitutive relationship between the macroscopic average stress $\bar{\sigma}$ and the macroscopic average strain $\bar{\epsilon}$, i.e.,

$$ \bar{\sigma} = C : \bar{\epsilon} \qquad (20) $$

can be written in its most general anisotropic form by using Voigt's notation for the 2D case as

$$ \begin{pmatrix} \bar{\sigma}_{11} \\ \bar{\sigma}_{22} \\ \bar{\sigma}_{12} \end{pmatrix} = \begin{pmatrix} C_{11} & C_{12} & C_{13} \\ C_{21} & C_{22} & C_{23} \\ C_{31} & C_{32} & C_{33} \end{pmatrix} \begin{pmatrix} \bar{\epsilon}_{11} \\ \bar{\epsilon}_{22} \\ 2\bar{\epsilon}_{12} \end{pmatrix} \qquad (21) $$

Here, $C_{ij}$ denotes the components of the homogenized effective tangent stiffness $C$. Analyzing the auxetic structures considered in this study, it is realized that for all four structures, the homogenized effective tangent stiffness exhibits the following form within the stress-strain relationship:

$$ \begin{pmatrix} \bar{\sigma}_{11} \\ \bar{\sigma}_{22} \\ \bar{\sigma}_{12} \end{pmatrix} = \begin{pmatrix} C_{11} & C_{12} & 0 \\ C_{12} & C_{11} & 0 \\ 0 & 0 & C_{33} \end{pmatrix} \begin{pmatrix} \bar{\epsilon}_{11} \\ \bar{\epsilon}_{22} \\ 2\bar{\epsilon}_{12} \end{pmatrix} \qquad (22) $$

Such an orthotropic behavior is due to the presence of the inherent symmetry axes within the considered auxetic unit cells, as previously demonstrated for chiral auxetic structures in [19]. It is also to be noted that in contrast to isotropic materials, in which $C_{33} = (C_{11} - C_{12})/2$, such an equality does not hold for the auxetic structures at hand. Therefore, as can be realized from Eq. (22), the three distinct constants $C_{11}$, $C_{12}$, and $C_{33}$ are required to fully define the homogenized effective tangent stiffness for the given auxetic structures. However, as mentioned in Section 3.1, due to the linear effect of Young's modulus $E$ on the effective elastic stiffness, the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ are considered as the outputs of the surrogate model. This normalization allows for a more generalized representation of the unit cell's behavior, independent of the specific value of Young's modulus. Besides, the effect of Young's modulus can be incorporated by applying it as a simple scaling factor to the outputs of the generated surrogate models after prediction.

### 3.3. Random forests regression

Among others, the random forests regression method [57] is used in the current work to generate surrogate models for the four auxetic unit cells presented in Fig. 1. This choice was made due to the method's suitability for handling nonlinear relationships with a relatively small number of input parameters (three in the present work), its ease of implementation with minimal hyperparameter tuning, and its ability to provide accurate predictions with minimal overfitting. For each unit cell, three separate surrogate models are generated, one for each of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$, resulting in a total of twelve surrogate models. Each model predicts its respective effective elastic constant based on given values of relative void diameters $d/L$ and $D/L$, as well as Poisson's ratio of the base material $\nu$. Random forest regression is an advanced machine learning method that combines multiple decision trees to generate a powerful predictive model. At the training stage, the method creates multiple decision trees, each on a different subset of data, whose predictions are then averaged as the output of the model at the prediction stage. Thus, this makes the random forests approach quite robust in nature by mitigating overfitting of individual decision trees, thereby improving generalization to unseen data.

The FFT-based homogenization scheme is used to generate datasets comprising 48k (i.e., 48000) samples for each of the four unit cells. FFT-based simulations are performed on a 256×256 grid discretization of the geometry to compute the normalized effective elastic constants (through Eqs. (15) and (17)) using the MINRES solver with a tolerance of $10^{-6}$. The 48k data points are generated by randomly sampling from the input
<figure><img src="image_6.png" /><figcaption>Fig. 2. Plots of the predicted values (obtained from the surrogate models) versus the actual values (obtained from the FFT-based solver) of the normalized effective elastic constants C₁₁/E, C₁₂/E, and C₃₃/E for the rectangular void unit cell. The diagonal dashed line depicts the perfect prediction.</figcaption></figure>

parameter space. This space is defined by three key parameters: Poisson's ratio of the base material $v$, ranging from 0.2 to 0.4, and the relative void parameters $d/L$ and $D/L$, both spanning from 0.0 to 1.0. An additional geometric constraint is imposed such that $d/L + D/L < 1$. This constraint ensures that the voids within the unit cell do not interconnect, maintaining the structural integrity of the unit cell. This comprehensive sampling strategy, incorporating both parameter ranges and geometric constraints, ensures that the surrogate models are trained on a diverse yet physically meaningful set of unit cell configurations.

It is worth highlighting that the utilization of such voluminous datasets is justified by the remarkable efficiency of the FFT-based solver, which enables the production of such large volumes of data in a relatively short period of time. To be more specific, it just took a maximum of around 8 hours of wall-clock time on a computing system equipped with two Intel Xeon Platinum 8160 processors, totaling 48 physical cores (i.e., $2 \times 24$) to generate such a large dataset for each of the four unit cells. This time includes both the geometry generation and the FFT-based homogenization procedures. Notably, the 8-hour duration represents the maximum time required, observed for the peanut-shaped void unit cell due to its complex geometry that demands more time for generation. In comparison, the rectangular void unit cell, with its simpler geometry, required significantly less time, taking just over 4 hours for both geometry generation and FFT-based homogenization. Additionally, if only the full-field solution for a given macroscopic strain field is required, and not the homogenized effective tangent stiffness, the solver time (excluding geometry generation) is reduced to approximately one-third of the time needed for tangent stiffness evaluation. Moreover, employing a substantial dataset size is crucial for achieving highly accurate surrogate models that yield smooth solution fields (see

Section 3.4). Nonetheless, studies on the influence of dataset size on the accuracy of the surrogate models indicate that even significantly smaller datasets can produce models with reasonable accuracy. The effect of dataset size on the accuracy of the trained surrogate models is investigated in more detail in Section 3.4.1.

For the training of the surrogate models, the following mean squared error $\mathcal{L}_{\text{model}}$ is employed as the loss function:

$$ \mathcal{L}_{\text{model}} = \frac{1}{n_s} \sum_{i=1}^{n_s} (z_i - \hat{z}_i)^2, \qquad (23) $$

where $n_s$ denotes the number of samples, $z_i$ is the $i$-th actual value obtained from the FFT-based solver, and $\hat{z}_i$ is the $i$-th predicted value obtained from the trained model. In this context, the values $z_i$ and $\hat{z}_i$ specifically correspond to the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, or $C_{33}/E$, depending on the particular model being trained. Moreover, the performance of the models is evaluated using the coefficient of determination $R^2$, which is defined as

$$ R^2 = 1 - \frac{\sum_{i=1}^{n_s} (z_i - \bar{z}_i)^2}{\sum_{i=1}^{n_s} (z_i - \bar{z}_i)^2}. \qquad (24) $$

Here, $\bar{z}$ is the mean of the actual values, which can be expressed by

$$ \bar{z} = \frac{1}{n_s} \sum_{i=1}^{n_s} z_i. \qquad (25) $$

The key hyperparameters used for the random forests regressor from the scikit-learn library are listed in Table 1. To evaluate the predictive capability of the generated surrogate models, a train/test split methodology is employed. The dataset is randomly partitioned using
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Plots of the predicted values (obtained from the surrogate models) versus the actual values (obtained from the FFT-based solver) of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ for the diamond void unit cell. The diagonal dashed line depicts the perfect prediction.</figcaption></figure>

a 90%/10% distribution, with 90% of the data allocated for training and 10% reserved for testing. Figs. 2–5 illustrate the performance of the trained surrogate models in predicting the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ for unit cells with rectangular, diamond, oval, and peanut-shaped voids, respectively. In each figure, subfigure (a) presents the results for the training set, while subfigure (b) shows the results for the test set. In these figures, the actual values and the corresponding predictions are plotted on the x and y axes, respectively, with the diagonal line (i.e., $x = y$) depicting the perfect prediction. Furthermore, to quantify the accuracy of the models, the coefficient of determination $R^2$ is employed as the primary criterion and prominently displayed up to four decimal places within the figures. As evident from these plots, the data points show minimal deviation from the diagonal, indicating remarkably high accuracy across all surrogate models, regardless of the unit cell geometry or the specific effective elastic constant being predicted. As anticipated, the models demonstrate superior performance on the training sets, given their prior exposure to this data. Nevertheless, the test sets also display promising accuracy despite being unseen, with $R^2$ values no lower than 0.9993 for all cases. This consistently high level of performance across both seen and unseen data, and across all unit cell types and elastic constants, underscores the robustness and generalizability of the generated surrogate models.

Among different unit cell types, the diamond void unit cell exhibits the best performance, as evidenced by both the minimal deviation of data points from the diagonal line and the highest coefficient of determination $R^2$. In contrast, while still maintaining solid predictive capabilities, the oval void unit cell shows slightly more variability in its predictions. Nevertheless, it is noteworthy that even for the oval void models, the $R^2$ values do not fall under 0.9993, highlighting the con-

sistently high precision of the generated surrogate models across all four unit cell geometries. When comparing the performance of the surrogate models for different effective elastic constants, the predictions for $C_{33}/E$ consistently achieve the highest accuracy across all unit cell types. Conversely, while $C_{12}/E$ predictions show the largest deviations from the ideal line, the generated models still achieve reliable accuracy with $R^2$ values no lower than 0.9993. This demonstrates the overall reliability of the generated surrogate models for all three effective elastic constants. It is worth noting that the high precision of the developed surrogate models can be attributed to the limited number of input parameters, consisting of two geometric parameters and a single material input parameter. Increasing the number of input parameters, as necessary, e.g., when considering more complex geometries, may reduce model accuracy. For instance, surrogate models developed for the complex geometry of tetra-chiral auxetics [29], parameterized with eight design parameters [58], demonstrated slightly less accuracy.

It is also interesting to note that all four unit cells exhibit very similar tunable ranges for their effective elastic constants. The maximum values across all unit cells are approximately $C_{11}/E \approx 2.14$, $C_{12}/E \approx 1.43$, and $C_{33}/E \approx 0.42$. While the minimum values for $C_{11}/E$ and $C_{33}/E$ approach the theoretical limit of zero, $C_{12}/E$ has a minimum value around -0.15. These results indicate that each void shape provides nearly identical flexibility in tuning the mechanical properties.

### 3.4. A parametric study

In this section, a parametric study is performed to both evaluate the accuracy of the developed surrogate models and examine the auxetic behavior of the unit cells for different geometric parameters. To this
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Plots of the predicted values (obtained from the surrogate models) versus the actual values (obtained from the FFT-based solver) of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ for the oval void unit cell. The diagonal dashed line depicts the perfect prediction.</figcaption></figure>

<table><caption>**Table 1**
Key hyperparameters for the random forests regressor from the scikit-learn library.</caption><thead><tr><td>Hyperparameter</td><td>Value</td><td>Description</td></tr></thead><tbody><tr><td>n_estimators</td><td>100</td><td>Number of trees in the forest.</td></tr><tr><td>max_depth</td><td>None</td><td>Maximum depth of each tree (None indicates fully grown trees).</td></tr><tr><td>min Samples Split</td><td>2</td><td>Minimum number of samples required to split a node.</td></tr><tr><td>max_features</td><td>sqrt</td><td>Number of features to consider when looking for the best split (sqrt refers to the square root of the total features).</td></tr><tr><td>bootstrap</td><td>True</td><td>Specifies whether bootstrap samples are used when building trees.</td></tr></tbody></table>

end, by fixing the Poisson's ratio of the base material as $v = 0.3$ and one of the relative void diameters as $d/L = 0.05$, the other relative void diameter $D/L$ is changed within the interval $[0.05, 0.9]$ with a step size of 0.01.

The normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ are plotted versus the relative void diameter $D/L$ in Figs. 6–9 for all four unit cells. It is observed that for all three effective elastic constants and all four structures, the surrogate models provide results in close agreement with those obtained from the FFT-based solver. The curves exhibit a consistent overall trend, indicating a general similarity in the effect of the void size (i.e., $D/L$) on the effective elastic constants, whereas subtle differences are apparent in the slopes of these curves. Interestingly, despite the slight differences in the slopes, the maximum and minimum values of each effective elastic constant remain comparable across different auxetic unit cells.

The auxeticity of the structures is identified from the negative values of $C_{12}$ in Figs. 6–9 for a specific range of the void diameter $D/L$. In order to provide a clearer perspective on the auxeticity of the structures for different void sizes, the effective Poisson's ratio is computed

as $v_{eff} = C_{12}/C_{11}$ and plotted versus the relative void diameter $D/L$ in Fig. 10 for all four unit cells. It is observed that for the unit cells with rectangular, diamond, oval, and peanut-shaped voids, the structures enter the auxetic region for the relative void diameters $D/L$ equal and greater than 0.42, 0.46, 0.45, and 0.44, respectively. It is worth highlighting that both the FFT-based solver and surrogate models yield identical values of the aforementioned critical void diameters, thereby demonstrating the high accuracy of the generated surrogate models. However, it can be seen that for the rectangular void and peanut-shaped void unit cells, the accuracy of the surrogate models slightly decreases when the relative void diameter $D/L$ approaches 0.9. This can be well justified by the fact that by increasing the void diameter, the effective elastic constants approach zero as the amount of material decreases in the structure (see Figs. 6–9). For such large values of the relative void diameter $D/L$, the effective Poisson's ratio $v_{eff} = C_{12}/C_{11}$ is obtained by dividing two small numbers close to zero. Therefore, even a small error in each of the elastic constants $C_{12}$ and $C_{11}$ may lead to a higher error in their fraction. In case there is a need for highly accurate solutions for the effective
<figure><img src="image_9.png" /><figcaption>Fig. 5. Plots of the predicted values (obtained from the surrogate models) versus the actual values (obtained from the FFT-based solver) of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ for the peanut-shaped void unit cell. The diagonal dashed line depicts the perfect prediction.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Plots of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ versus the relative void diameter $D/L$ for the rectangular void unit cell. Poisson's ratio of the base material and the other relative void diameter are fixed as $v = 0.3$ and $d/L = 0.05$, respectively.</figcaption></figure>

Poisson's ratio $v_{\text{eff}}$, a separate surrogate model to predict this value can be generated as well.

To illustrate the effect of geometric parameters on the anisotropy of the effective Poisson's ratio $v_{\text{eff}}$, Fig. 11 depicts polar plots of the effective Poisson's ratio $v_{\text{eff}}$ with different void diameter ratios of $D/d = 1, 6, 12$, and $18$ (computed by keeping $d/L = 0.05$ constant and increasing $D/L$) for all four considered unit cells. To avoid presenting an excessive number of figures, only the results of the FFT-based solver are shown due to their close similarity with the results obtained using the surrogate models. As evident from the plots, all unit cells exhibit orthotropic behavior, with the effective Poisson's ratio $v_{\text{eff}}$ in the x and y directions ($\theta = 0^\circ$ and $90^\circ$) being identical. For all unit cells, increasing the void diameter ratio $D/d$ leads to a greater degree of anisotropy

and enhanced auxetic behavior in the x- and y-directions. For smaller $D/d$ ratios, the effective Poisson's ratio $v_{\text{eff}}$ can become negative along the diagonal direction ($\theta = 45^\circ$), while remaining positive in the x and y directions. However, as $D/d$ increases, the auxetic behavior in the x- and y-directions becomes more pronounced, while it decreases along the diagonal direction.

In the previously presented figures in this section (i.e., Figs. 6–10), the effective elastic constants and Poisson's ratio were shown with a fixed relative void diameter $d/L$ and a variable $D/L$. To examine the effect of varying both void diameters on the effective elastic constants and Poisson's ratio, the distribution of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$, along with the effective Poisson's ratio $v_{\text{eff}} = C_{12}/C_{11}$, is plotted as a function of the relative void diam-
<figure><img src="image_11.png" /><figcaption>Fig. 7. Plots of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ versus the relative void diameter $D/L$ for the diamond void unit cell. Poisson's ratio of the base material and the other relative void diameter are fixed as $v = 0.3$ and $d/L = 0.05$, respectively.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 8. Plots of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ versus the relative void diameter $D/L$ for the oval void unit cell. Poisson's ratio of the base material and the other relative void diameter are fixed as $v = 0.3$ and $d/L = 0.05$, respectively.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 9. Plots of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$ versus the relative void diameter $D/L$ for the peanut-shaped void unit cell. Poisson's ratio of the base material and the other relative void diameter are fixed as $v = 0.3$ and $d/L = 0.05$, respectively.</figcaption></figure>

eters $d/L$ and $D/L$ for both the FFT-based solver (Fig. 12a) and the surrogate model (Fig. 12b). As the results are closely similar for all four unit cells, only the plots for the rectangular void unit cell are shown to avoid an unnecessary number of figures. A close agreement can be observed between the results obtained from the surrogate models and those from the FFT-based solver. Additionally, the combinations of the relative void diameters $d/L$ and $D/L$ resulting in specific effective elastic constants or Poisson's ratio can be identified from the presented plots.

### 3.4.1. Dataset size

As mentioned in Section 3.3, the surrogate models are generated from datasets of size 48k. The reasons to consider such large datasets are twofold: Firstly, the utilized FFT-based solver is so efficient that generating such a large volume of data is possible within a few hours (i.e., around 8 hours of wall-clock time with 48 core processors; see Section 3.3). Secondly, utilizing a large volume of data ensures obtaining

highly accurate surrogate models that can provide smooth solution fields, such as those presented in Figs. 6–9.

Nevertheless, investigating the effect of dataset size on the accuracy of the surrogate models reveals that even significantly smaller datasets can still generate adequately accurate surrogate models. To explain more thoroughly, surrogate models for the $C_{11}$ component of the effective elastic stiffness are generated using dataset sizes of 3k, 6k, 12k, and 24k in addition to the already available models with 48k data points. Then, these new surrogate models are utilized for the same problem of Section 3.4 to compute the normalized elastic constant $C_{11}/E$ for all four unit cells. Completely similar to the previous parametric study (Section 3.4), the models are evaluated for fixed values of $v = 0.3$ and $d/L = 0.05$ while $D/L$ is varied in the range [0.05, 0.9] with a step size of 0.01.

The relative errors between the values obtained from the surrogate models trained on varying dataset sizes and the FFT-based solver, de-
<figure><img src="image_14.png" /><figcaption>Fig. 10. Plots of the effective Poisson's ratio $v_{eff}$ versus the relative void diameter $D/L$ for the auxetic unit cells with rectangular, diamond, oval, and peanut-shaped voids. The Poisson's ratio of the base material and the other relative void diameter are fixed as $v = 0.3$ and $d/L = 0.05$, respectively.</figcaption></figure>

fined as $|\frac{z_i - \hat{z}_i}{z_i}| \times 100$, are computed for all the values of relative void diameter $D/L$ and presented as box plots in Fig. 13 for all four unit cells. A noticeable trend is the reduction in the median relative error as the dataset size increases, indicating an overall improvement in model accuracy. Additionally, the interquartile range (IQR) becomes narrower, highlighting a decrease in the variability of the relative error. The length of the whiskers also reduces, suggesting a tightening in the overall spread of the error distribution. It is worth noting that, although utilizing a large dataset consisting of 48k data points can result in higher precision with a relative error distribution below 2%, even significantly reduced datasets, such as a dataset consisting of 3k data points, are capable of producing reasonably accurate surrogate models with a relative error distribution around or below 6%.

Regarding the performance of surrogate models trained on the smallest dataset size (i.e., 3k data points) for different unit cell types, a comparison of the box plots reveals that the oval void unit cell achieves the highest accuracy with smaller datasets, exhibiting a relative error distribution below 4%. In contrast, the peanut-shaped void unit cell shows the highest error, slightly exceeding 6%. The rectangular and diamond void unit cells demonstrate similar performance, with errors just below 6%. However, for all four unit cells, using the dataset with 48k data points results in highly accurate models with relative errors below 2%.

## 4. Inverse analysis

The possibility of evaluating the generated surrogate models in real time enables exploring the design space and finding the specific geometric and material input parameters that result in the desired effective stiffness components. For instance, see Section 5, where even an inexperienced user can probe the design space in real time using a newly developed GUI.

While the real-time evaluation of the surrogate models is highly advantageous, establishing an inverse analysis framework utilizing these

models remains highly beneficial. In such inverse approaches, input parameters that yield desired outputs are determined through iterative schemes, such as gradient-based or heuristic optimization methods. In the present paper, the efficiency of real-time model evaluation provides an opportunity to implement an iterative grid refinement strategy, solving the inverse problem with high precision and minimal computational cost.

This iterative grid refinement approach begins by evaluating the model for a coarse grid that covers the entire parameter space, offering comprehensive initial exploration. Within this grid, regions with high potential for optimal results are identified. At each refinement stage, the grid is adaptively redefined to focus on these promising areas, progressively concentrating on regions likely to meet the target criteria. By iterating in this targeted manner, the approach efficiently minimizes evaluations in less relevant parts of the design space. The refinement continues until a specified accuracy is achieved or after completing a set number of stages, resulting in a flexible and effective search method. This adaptive grid refinement approach fully leverages the speed of the surrogate models, producing optimal solutions with focused exploration. A numerical example in Section 4.1 demonstrates the effectiveness of this method.

In the current work, considering that the base material for the fabrication of the structures, and thereby its Poisson's ratio $v$, is known, the inverse problems involves finding the optimal relative void diameters $d/L$ and $D/L$ that result in the desired values for the components of the normalized effective tangent stiffness $\hat{y}_i$, defined as

$$ \hat{y}_1 = \frac{\hat{C}_{11}}{E}, \quad \hat{y}_2 = \frac{\hat{C}_{12}}{E}, \quad \hat{y}_3 = \frac{\hat{C}_{33}}{E}. \qquad (26) $$

Here, $\hat{\cdot}$ is used to denote the desired output values. The predicted values $y_i$ are then obtained by evaluating the corresponding surrogate models $S_i$ for each component of the effective stiffness matrix as
<figure><img src="image_15.png" /><figcaption>Fig. 11. Polar plots of the effective Poisson's ratio $v_{eff}$ with different void diameter ratios of $D/d = 1, 6, 12$, and 18 for the auxetic unit cells with rectangular, diamond, oval, and peanut-shaped voids. The Poisson's ratio of the base material is set to $v = 0.3$, the relative void diameter $d/L$ is kept constant at 0.05, and the other relative void diameter $D/L$ is adjusted to achieve different $D/d$ ratios.</figcaption></figure>

$$ 
\begin{aligned}
y_1 &= \frac{C_{11}}{E} = S_1 \left( \frac{d}{L}, \frac{D}{L}, v \right), & y_2 &= \frac{C_{12}}{E} = S_2 \left( \frac{d}{L}, \frac{D}{L}, v \right), \\
y_3 &= \frac{C_{33}}{E} = S_3 \left( \frac{d}{L}, \frac{D}{L}, v \right).
\end{aligned} 
\quad (27) $$

Then, by computing the mean squared error $L_{inv}$ between the predicted values $y_i$ and the desired output values $\hat{y}_i$ using

$$ L_{inv} \left( \frac{d}{L}, \frac{D}{L}, v \right) = \frac{1}{n_c} \sum_{i=1}^{3} w_i (y_i - \hat{y}_i)^2, \quad (28) $$

the optimal set of input parameters $d/L$ and $D/L$ resulting in an error smaller than a chosen threshold can be extracted. In the above, $n_c = \sum_{i=1}^{3} w_i$, where $w_i$ is a boolean variable indicating whether a particular component or any combination of the components of the effective stiffness matrix (7 possible combinations) is specified as the desired output, and it is defined as follows:

$$ w_i = \begin{cases} 1 & \text{if } \hat{y}_i \text{ is specified} \\ 0 & \text{if } \hat{y}_i \text{ is not specified} \end{cases} \quad (29) $$

When specifying the desired output as a combination of the components of the effective stiffness matrix, it is important to note that the solution to the inverse problem for such a combination may or may not exist within the defined error threshold. Conversely, if only a particular

component is specified as the desired output, the likelihood of the inverse problem having a solution is higher, provided that the desired value is both physical and possible. A practical example with specifying two components of the effective tangent stiffness as the desired output can be followed in Section 4.1 for the rectangular void unit cell.

### 4.1. An inverse problem

In this section, an inverse problem is solved to examine the performance of the generated surrogate models for inverse design. It is considered that the base material for fabricating the structure is the aditively manufactured poly-lactic acid (PLA) with the elastic properties $E = 3500$ MPa and $v = 0.36$ chosen from [59]. The goal here is to obtain the relative void diameters $d/L$ and $D/L$ for the rectangular void unit cell that yield the desired effective stiffness components $C_{11} = 400$ MPa and $C_{12} = -200$ MPa. Such values of the effective stiffness components result in a negative effective Poisson's ratio $v_{eff} = C_{12}/C_{11} = -0.5$.

The inverse problem is solved using the iterative grid refinement approach described in Section 4. As outlined, the accuracy of the solutions obtained depends on the density of the initial coarse grid and the parameters used for refinement. In this case study, the initial coarse grid consists of 400 evaluation points, providing comprehensive coverage of the parameter space in the first round. At each refinement stage,
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Distribution of the normalized effective elastic constants $C_{11}/E$, $C_{12}/E$, and $C_{33}/E$, along with the effective Poisson's ratio $v_{eff} = C_{12}/C_{11}$ as a function of the relative void diameters $d/L$ and $D/L$ for the rectangular void unit cell. The Poisson's ratio of the base material is set to $v = 0.3$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Box plots for the distribution of the relative error in the normalized effective elastic constant $C_{11}/E$ between the actual values obtained from the FFT-based solver and the predicted values obtained from the surrogate models trained on varying dataset sizes. Results are presented for the auxetic unit cells with rectangular, diamond, oval, and peanut-shaped voids.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 14. Geometry and boundary conditions for the specimens made form a 10 × 10 distribution of the rectangular void unit cell with geometric parameters obtained form the inverse problem. In both cases, the left and bottom boundaries are fixed in the x and y directions, respectively, and a displacement of $u_x = 1$ mm is applied to the right boundary. This is while, the top boundary is free for the first case (subfigure a) and fixed in the y direction for the second case (subfigure b).</figcaption></figure>

the grid is adaptively focused around the top-performing 5% of points with the lowest loss values, reducing the search area by a factor of 0.2 around these promising regions. The refinement process continues iteratively until a predefined loss threshold of $10^{-5}$ is met. Although the 48k dataset is available, the 3k dataset is used here to solve the inverse problem, showcasing the effectiveness of surrogate models even with a smaller dataset. The inverse problem is solved efficiently with 489 evaluations and a CPU time of 0.03 seconds, yielding relative void diameters of $d/L = 0.34$ and $D/L = 0.50$. To validate, the inverse problem is also solved with the 48k dataset, producing the same values.

After obtaining the values $d/L = 0.34$ and $D/L = 0.50$ from the inverse problem, a 10 × 10 distribution of the rectangular void unit cell with a side length of 10 mm is simulated with FEM to examine the accuracy of the obtained design parameters. As shown in Fig. 14, two different sets of boundary conditions are considered to obtain $v$, $C_{11}$ and $C_{12}$ from FE simulations.

In the first case (Fig. 14a), the left and bottom boundaries of the specimen are fixed in the x and y directions, respectively, and a displacement of $u_x = 1$ mm is applied to the right boundary. Then, by computing the displacement of the transverse boundary, i.e., the top boundary, the effective Poisson's ratio can be computed by

$$v_{\text{eff}} = -\frac{\bar{\epsilon}_{yy}}{\bar{\epsilon}_{xx}} \quad \text{with} \quad \bar{\epsilon}_{xx} = \frac{u_x}{L_x}, \quad \bar{\epsilon}_{yy} = \frac{u_y}{L_y}, \qquad (30)$$

where the effective strains of the entire specimen in the x and y directions are represented by $\bar{\epsilon}_{xx}$ and $\bar{\epsilon}_{yy}$, respectively. These strains are related to the displacements of the specimen's boundaries: $u_x$ refers to the displacement of the right boundary in the x direction, while $u_y$ corresponds to the displacement of the top boundary in the y direction. Moreover, $L_x$ and $L_y$ denote the widths of the specimen in the x and y directions, respectively.

In the second case (Fig. 14b), the boundary conditions are similar to those in the first case but the top boundary is additionally fixed in the y direction. Here, the effective stiffness components $C_{11}$ and $C_{12}$ can be computed by

$$C_{11} = \frac{\bar{\sigma}_{xx}}{\bar{\epsilon}_{xx}}, \quad C_{12} = \frac{\bar{\sigma}_{yy}}{\bar{\epsilon}_{xx}} \quad \text{with} \quad \bar{\epsilon}_{xx} = \frac{f_x}{tL_y}, \quad \bar{\epsilon}_{yy} = \frac{f_y}{tL_x}, \qquad (31)$$

where $t$ is the thickness of the specimen. Furthermore, the sum of reaction forces in the x direction on the right boundary is denoted by $f_x$, while $f_y$ represents the sum of reaction forces in the y direction on the top boundary.

The obtained values from FE simulations and the corresponding error with respect to the desired values are summarized in Table 2. It

<caption>Table 2</caption>

Comparison of the obtained values from FE simulations using the geometric parameters achieved from the inverse problem for the rectangular void unit cell against the desired values and their corresponding relative error.

<table><thead><tr><th>Quantity</th><th>Desired value</th><th>Obtained value (FEM)</th><th>Error</th></tr></thead><tbody><tr><td>$v_{\text{eff}}$ [–]</td><td>–0.500</td><td>–0.500</td><td>0.00%</td></tr><tr><td>$C_{11}$ [MPa]</td><td>400.00</td><td>402.69</td><td>0.67%</td></tr><tr><td>$C_{12}$ [MPa]</td><td>–200.00</td><td>–199.89</td><td>0.06%</td></tr></tbody></table>

can be seen that the effective Poisson's ratio $v_{\text{eff}}$ shows an exact match between the desired and obtained values up to three decimal places, resulting in a zero error. For the effective elastic constants $C_{11}$ and $C_{12}$, the obtained values from FEM deviate slightly from the desired values with percentage errors of 0.67% and 0.06%, respectively, indicating a high degree of accuracy in the obtained design parameters from the inverse problem.

## 5. Graphical user interface (GUI)

To facilitate the practical usage of the generated surrogate models for less experienced users, a GUI has been developed using Python's Tkinter library. The developed GUI consists of two tabs, including the **Surrogate Model** tab (Fig. 15a) and the **Inverse Analysis** tab (Fig. 15b).

In the Surrogate Model tab, users can extract the effective tangent stiffness $C$ and the resulting effective Poisson's ratio $v_{\text{eff}}$ in real time for the desired auxetic structure. Users can change the geometric input parameters $d/L$ and $D/L$ as well as the material input parameters $E$ and $v$. Moreover, the unit cell geometry is plotted in real time for users.

In the Inverse Analysis tab, users can extract the relative void diameters $d/L$ and $D/L$ for a selected unit cell, given the desired components of the effective elastic stiffness and the base material properties $E$ and $v$. Additionally, the corresponding unit cell geometry with the obtained geometric parameters is plotted. It is worth mentioning that, as detailed in Section 4, users are given the flexibility to select any component or any combination of the components of the effective elastic stiffness as the target output. However, it is important to recognize the inherent constraints of the inverse problem. Users should be reminded that not all values of desired effective elastic properties are achievable within the physical limitations of the unit cell design. The space of all possible solutions is bounded by the base material properties and geometric constraints of the unit cell structure. Therefore, if the desired values lie outside this feasible region, the inverse analysis fails to provide a valid solution. Users are advised to consider these physical limitations when
<figure><img src="image_19.png" /><figcaption>(a) Surrogate Model tab for the real-time prediction of the effective tangent stiffness.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Inverse Analysis tab for finding the relative void diameters *d/L* and *D/L*.</figcaption></figure>

<caption>Fig. 15. GUI for the real-time prediction of the effective tangent stiffness and the inverse analysis of the auxetic unit cells.</caption>

setting the desired outputs to ensure that the inverse problem remains
solvable.

**6. Conclusions**

The present study focused on developing surrogate models for real-
time prediction of the effective elastic properties of auxetic unit cells
with orthogonal voids of different shapes, including rectangular, dia-
mond, oval, and peanut-shaped voids. These surrogate models accept

the relative void diameters *d/L* and *D/L* along with Poisson's ratio of the base constituent material *v* as input parameters and deliver the normalized effective elastic constants *C*₁₁/*E*, *C*₁₂/*E*, and *C*₃₃/*E* as their outputs. These three constants represent the only distinct values necessary to fully define the homogenized effective elastic stiffness tensor.

The FFT-based homogenization scheme was employed to efficiently
generate the datasets required for training the surrogate models. This
particular solver facilitated the rapid generation of extensive datasets
within a relatively short time frame, circumventing the challenges typ-
ically associated with FEM. More specifically, it eliminated the need for careful consideration of symmetrical meshes at the unit cell boundaries and the application of periodic boundary conditions, which are often cumbersome in FE approaches. The efficiency of this approach was demonstrated by the generation of a comprehensive dataset comprising 48k data points, accomplished within a maximum of 8 hours of wall-clock time using a computer with 48 CPU cores.

The random forests regression method was used to develop surrogate models based on the datasets generated by the FFT-based solver. The performance of the generated surrogate models was examined through a train/test split methodology. By comparing the predicted values obtained from the models with the actual values obtained from the FFT-based solver, it was demonstrated that the generated surrogate models are capable of achieving a high level of accuracy in the prediction of the effective elastic constants even for the test set, which contained unseen data. In addition, a parametric study was also conducted to evaluate the accuracy of the developed surrogate models and investigate the auxetic behavior of the unit cells across diverse geometric configurations. This analysis not only showcased the high precision of the surrogate models but also provided valuable insights into the relationship between structural geometry and auxetic properties. The study investigated the impact of dataset size on the accuracy of the surrogate models, revealing that even substantially reduced datasets can produce satisfactory results.

An inverse design framework has been established to determine the geometric parameters that yield desired effective stiffness constants. This inverse problem was solved using an iterative grid refinement algorithm, which efficiently explores the design space by progressively focusing on promising regions, minimizing evaluations in less relevant areas and leveraging the speed of surrogate models to extract the optimal set satisfying the desired requirements. To assess the capability of the established surrogate models in inverse design, a practical inverse problem was solved for the rectangular void unit cell. Despite using the small 3k-sized dataset, the inverse problem was solved efficiently with a short CPU time of 0.03 seconds, achieving the same results as the 48k-sized dataset, which demonstrates the effectiveness of surrogate models with minimal data. FE simulations were performed using the obtained geometric parameters from the inverse problem to examine the performance. It was shown that the obtained values from FE simulations, i.e., the effective Poisson's ratio and the effective elastic constants, were in close agreement with the target values for the inverse problem with a maximum error of 0.67%. This underscores the benefit of the generated surrogate models in designing auxetic structures with tailored mechanical properties.

A user-friendly GUI has also been developed to enhance the accessibility and practical application of surrogate models for auxetic metamaterials. The GUI features two main components: a tab for real-time prediction of effective elastic constants based on user-defined geometric and material inputs, and a second tab for determining the optimal unit cell geometry to achieve desired effective elastic properties. This tool not only facilitates rapid effective properties prediction and inverse design but also provides visual feedback through real-time unit cell geometry plotting.

Moving forward, this study can be extended in several directions for subsequent research. Developing surrogate models capable of predicting full-field strain/stress distributions would be valuable for applications where such detailed information is crucial. These models would incorporate macroscopic strain field values in addition to geometric and material input parameters. Another possible extension is to further increase the tunability of the considered auxetic unit cells by replacing the void areas with soft inclusions [20,60]. For such composite auxetic unit cells, the ratio of the Young's moduli between the two phases is also considered an input parameter, thereby expanding the design space. Moreover, creating surrogate models for the finite deformation behavior of both auxetic and non-auxetic lattice structures could significantly reduce computational effort in two-scale homogenization approaches. Such models would take the macroscopic deformation

gradient as input and yield the corresponding homogenized tangent stiffness at each macroscopic integration point during incremental loading. Future work could extend this framework to nonlinear mechanical response, as demonstrated by studies that explore inverse design approaches to tailor nonlinear stress-strain responses for applications such as soft robotics and impact mitigation [61]. Incorporating auxetic geometries in energy harvesting systems has proven effective in enhancing energy output and efficiency [62], highlighting the potential benefits of expanding the current framework to include time-dependent and dynamic behaviors. Eventually, another avenue for future research is to establish an active learning framework with on-the-fly data generation, in which the FFT-based solver is used to generate data points dynamically during the training process. This approach would enhance the accuracy of the generated surrogate models and the quality of fit, while avoiding the need to generate an extensive dataset prior to the training stage, potentially leading to more efficient and adaptive model development.

### CRediT authorship contribution statement

**Hooman Danesh:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**Daniele Di Lorenzo:** Visualization, Software, Methodology, Investigation, Conceptualization.

**Francisco Chinesta:** Supervision, Resources, Project administration, Funding acquisition.
**Stefanie Reese:** Supervision, Resources, Project administration, Funding acquisition.
**Tim Brepols:** Writing – review & editing, Supervision, Resources, Funding acquisition, Project administration.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

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

The authors gratefully acknowledge the funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 956401 (XS-Meta).

### Data availability

Data will be made available on request.

### References

[1] J. Zhang, G. Lu, Z. You, Large deformation and energy absorption of additively manufactured auxetic materials and structures: a review, Composites, Part B, Eng. 201 (2020) 108340.

[2] M. Najafi, H. Ahmadi, G. Liaghat, Experimental investigation on energy absorption of auxetic structures, Mater. Today Proc. 34 (2021) 350–355.

[3] N.K. Choudhry, B. Panda, S. Kumar, In-plane energy absorption characteristics of a modified re-entrant auxetic structure fabricated via 3D printing, Composites, Part B, Eng. 228 (2022) 109437.

[4] Z. Yan, B. Wang, K. Wang, C. Zhang, A novel cellular substrate for flexible electronics with negative Poisson ratios under large stretching, Int. J. Mech. Sci. 151 (2019) 314–321.

[5] S. Jiang, X. Liu, J. Liu, D. Ye, Y. Duan, K. Li, Z. Yin, Y. Huang, Flexible metamaterial electronics, Adv. Mater. 34 (52) (2022) 2200070.

[6] B. Jang, S. Won, J. Kim, J. Kim, M. Oh, H.-J. Lee, J.-H. Kim, Auxetic meta-display: stretchable display without image distortion, Adv. Funct. Mater. 32 (22) (2022) 2113299.

[7] T. Allen, N. Martinello, D. Zampieri, T. Hewage, T. Senior, L. Foster, A. Alderson, Auxetic foams for sport safety applications, Proc. Eng. 112 (2015) 104–109.
[8] O. Duncan, T. Shepherd, C. Moroney, L. Foster, P.D. Venkatraman, K. Winwood, T. Allen, A. Alderson, Review of auxetic materials for sports applications: expanding options in comfort and protection, Appl. Sci. 8 (6) (2018) 941.

[9] D. Tahir, M. Zhang, H. Hu, Auxetic materials for personal protection: a review, Phys. Status Solidi (b) 259 (12) (2022) 2200324.

[10] H.M. Kolken, S. Janbaz, S.M. Leeflang, K. Lietaert, H.H. Weinans, A.A. Zadpoor, Rationally designed meta-implants: a combination of auxetic and conventional meta-biomaterials, Mater. Horiz. 5 (1) (2018) 28-35.

[11] H. Kolken, K. Lietaert, T. Van der Sloten, B. Pouran, A. Meynen, G. Van Loock, H. Weinans, L. Scheys, A. Zadpoor, Mechanical performance of auxetic biomaterials, J. Mech. Behav. Biomed. Mater. 104 (2020) 103658.

[12] M. Shirzad, A. Zolfagharian, M. Bodaghi, S.Y. Nam, Auxetic metamaterials for bone-implanted medical devices: recent advances and new perspectives, Eur. J. Mech. A, Solids 98 (2023) 104905.

[13] I. Masters, K. Evans, Models for the elastic deformation of honeycombs, Compos. Struct. 35 (4) (1996) 403-422.

[14] F. Dos Reis, J.-F. Ganghoffer, Equivalent mechanical properties of auxetic lattices from discrete homogenization, Comput. Mater. Sci. 51 (1) (2012) 314-321.

[15] L. Yang, O. Harrysson, H. West, D. Cormier, Mechanical properties of 3D re-entrant honeycomb auxetic structures realized via additive manufacturing, Int. J. Solids Struct. 69 (2015) 475-490.

[16] N. Chan, K. Evans, The mechanical properties of conventional and auxetic foams. Part I: compression and tension, J. Cell. Plast. 35 (2) (1999) 130-165.

[17] N. Chan, K. Evans, The mechanical properties of conventional and auxetic foams. Part II: shear, J. Cell. Plast. 35 (2) (1999) 166-183.

[18] F. Scarpa, P. Panayiotou, G. Tomlinson, Numerical and experimental uniaxial loading on in-plane auxetic honeycombs, J. Strain Anal. Eng. Des. 35 (5) (2000) 383-388.

[19] J. Dirrenberger, S. Forest, D. Jeulin, Effective elastic properties of auxetic microstructures: anisotropy and structural applications, Int. J. Mech. Mater. Des. 9 (2013) 21-33.

[20] D.M. Kochmann, G.N. Venturini, Homogenized mechanical properties of auxetic composite materials in finite-strain elasticity, Smart Mater. Struct. 22 (8) (2013) 084004.

[21] A. Slann, W. White, F. Scarpa, K. Boba, I. Farrow, Cellular plates with auxetic rectangular perforations, Phys. Status Solidi (b) 252 (7) (2015) 1533-1539.

[22] J.N. Grima, L. Mizzi, K.M. Azzopardi, R. Gatt, et al., Auxetic perforated mechanical metamaterials with randomly oriented cuts, Adv. Mater. 28 (2) (2016) 385-389.

[23] A.H. Tajalsir, K. Mustapha, T. Ibn-Mohammed, Numerical and random forest modelling of the impact response of hierarchical auxetic structures, Mater. Today Commun. 31 (2022) 103797.

[24] I. Ben-Yelun, G. Gómez-Carano, F.J. San Millán, M.Á. Sanz, F.J. Montáns, L. Saucedo-Mora, Gam: general auxetic metamaterial with tunable 3d auxetic behavior using the same unit cell boundary connectivity, Materials 16 (9) (2023) 3473.

[25] G.A. Lyngdoh, N.-K. Kelter, S. Doner, N.A. Krishnan, S. Das, Elucidating the auxetic behavior of cementitious cellular composites using finite element analysis and interpretable machine learning, Mater. Des. 213 (2022) 110341.

[26] S. Vyavahare, S. Teraiya, S. Kumar, FDM manufactured auxetic structures: an investigation of mechanical properties using machine learning techniques, Int. J. Solids Struct. 265 (2023) 112126.

[27] H. Wang, S.-H. Xiao, C. Zhang, Novel planar auxetic metamaterial perforated with orthogonally aligned oval-shaped holes and machine learning solutions, Adv. Mater. 23 (7) (2021) 2100102.

[28] C. Du, Y. Wang, Z. Kang, Auxetic Kirigami metamaterials upon large stretching, ACS Appl. Mater. Interfaces 15 (15) (2023) 19190-19198.

[29] Z. Liao, Y. Wang, L. Gao, Z.-P. Wang, Deep-learning-based isogeometric inverse design for tetra-chiral auxetics, Compos. Struct. 280 (2022) 114808.

[30] Y. Chang, H. Wang, Q. Dong, Machine learning-based inverse design of auxetic metamaterial with zero Poisson's ratio, Mater. Today Commun. 30 (2022) 103186.

[31] H. Liu, F. Hou, A. Li, Y. Lei, H. Wang, High-efficient and reversible intelligent design for perforated auxetic metamaterials with peanut-shaped pores, Int. J. Mech. Mater. Des. 19 (3) (2023) 553-566.

[32] J.-H. Bastek, S. Kumar, B. Telgen, R.N. Glaesener, D.M. Kochmann, Inverting the structure-property map of truss metamaterials by deep learning, Proc. Natl. Acad. Sci. 119 (1) (2022) e2111505119.

[33] M. Mohammadnejad, A. Montazeri, E. Bahmanpour, M. Mahnama, Artificial neural networks for inverse design of a semi-auxetic metamaterial, Thin-Walled Struct. 200 (2024) 111927.

[34] X.-L. Peng, B.-X. Xu, Data-driven inverse design of composite triangular lattice structures, Int. J. Mech. Sci. 265 (2024) 108900.

[35] X. Zheng, X. Zhang, T.-T. Chen, I. Watanabe, Deep learning in mechanical metamaterials: from prediction and generation to inverse design, Adv. Mater. 35 (45) (2023) 2302530.

[36] D. Lee, W. Chen, L. Wang, Y.-C. Chan, W. Chen, Data-driven design for metamaterials and multiscale systems: a review, Adv. Mater. 36 (8) (2024) 2305254.

[37] C. Zhang, Y.F. Zhao, A critical review on the application of machine learning in supporting auxetic metamaterial design, J. Phys. Mater. (2024).

[38] S. Lucarini, L. Cobian, A. Voitus, J. Segurado, Adaptation and validation of FFT methods for homogenization of lattice based materials, Comput. Methods Appl. Mech. Eng. 388 (2022) 114223.

[39] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and nonlinear mechanical properties of composites, C. R. Acad. Sci., Sér. 2, Méc. Phys. Chim. Astron. (1994).

[40] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Comput. Methods Appl. Mech. Eng. 157 (1-2) (1998) 69-94.

[41] M. Schneider, A review of nonlinear FFT-based computational homogenization methods, Acta Mech. 232 (6) (2021) 2051-2100.

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

[43] H. Danesh, T. Brepols, S. Reese, Challenges in two-scale computational homogenization of mechanical metamaterials, PAMM 23 (1) (2023) e202200139.

[44] J. Vondrejč, J. Zeman, I. Marek, An FFT-based Galerkin method for homogenization of periodic media, Comput. Math. Appl. 68 (3) (2014) 156-173.

[45] J. Vondrejč, J. Zeman, I. Marek, Guaranteed upper-lower bounds on homogenized properties by FFT-based Galerkin method, Comput. Methods Appl. Mech. Eng. 297 (2015) 258-291.

[46] F. Willot, Fourier-based schemes for computing the mechanical response of composites with accurate local fields, C. R., Méc. 343 (3) (2015) 232-245.

[47] J. Michel, H. Moulinec, P. Suquet, A computational scheme for linear and non-linear composites with arbitrary phase contrast, Int. J. Numer. Methods Eng. 52 (1-2) (2001) 139-160.

[48] S. Brisard, L. Dormieux, FFT-based methods for the mechanics of composites: a general variational framework, Comput. Mater. Sci. 49 (3) (2010) 663-671.

[49] V. Monchieta, G. Bonneta, A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast, Int. J. Numer. Methods Eng. 89 (11) (2012) 1419-1436.

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

[51] M.E.F. Idrissi, F. Praud, V. Champaney, F. Chinesta, F. Meraghni, Multiparametric modeling of composite materials based on non-intrusive PGD informed by multiscale analyses: application for real-time stiffness prediction of woven composites, Compos. Struct. 302 (2022) 116228.

[52] C. Gierden, J. Kochmann, J. Waimann, B. Svendsen, 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 (6) (2022) 4115-4135.

[53] F.S. Göküzüm, M.-A. Keip, An algorithmically consistent macroscopic tangent operator for FFT-based computational homogenization, Int. J. Numer. Methods Eng. 113 (4) (2018) 581-600.

[54] M. Rambausek, F.S. Göküzüm, L.T.K. Nguyen, M.-A. Keip, A two-scale FE-FFT approach to nonlinear magneto-elasticity, Int. J. Numer. Methods Eng. 117 (11) (2019) 1117-1142.

[55] V. Minh Nguyen-Thanh, L. Trong Khiem Nguyen, T. Rabczuk, X. Zhuang, A surrogate model for computational homogenization of elastostatics at finite strain using high-dimensional model representation-based neural network, Int. J. Numer. Methods Eng. 121 (2020) 4811-4842.

[56] S.Y. Reutskiy, The method of approximate fundamental solutions (MAFS) for elliptic equations of general type with variable coefficients, Eng. Anal. Bound. Elem. 36 (6) (2012) 985-992.

[57] L. Breiman, Random forests, Mach. Learn. 45 (2001) 5-32.

[58] Z.-P. Wang, Y. Wang, L.H. Poh, Z. Liu, Integrated shape and size optimization of curved tetra-chiral and anti-tetra-chiral auxetics using isogeometric analysis, Compos. Struct. 300 (2022) 116094.

[59] S. Vyavahare, S. Kumar, Numerical and experimental investigation of FDM fabricated re-entrant auxetic structures of ABS and PLA materials under compressive loading, Rapid Prototyping J. 27 (2) (2021) 223-244.

[60] X.-L. Peng, C. Soyarslan, S. Bargmann, Phase contrast mediated switch of auxetic mechanism in composites of infilled re-entrant honeycomb microstructures, Extrem. Mech. Lett. 35 (2020) 100641.

[61] J.-H. Bastek, D.M. Kochmann, Inverse design of nonlinear mechanical metamaterials via video denoising diffusion models, Nat. Mach. Intell. 5 (12) (2023) 1466-1475.

[62] S.K. Panda, J. Srinivas, Artificial neural network-based optimization studies for efficient energy harvesting from auxetic laminated composite smart beam, Mech. Adv. Mat. Struct. (2024) 1-17.