Deep learning for the rare-event rational design of 3D printed multi-material mechanical metamaterials

Architected and biomimetic materials can achieve exceptional combinations of mechanical properties through rational micro-architecture design, such as high stiffness with low weight or joint increases in stiffness and toughness. In metamaterials, targeting unusual property combinations like auxetic responses in multiple directions with high stiffness requires identifying specific micro-architectures, which is challenging due to the vast design space enabled by multi-material 3D printing. Multi-material printing allows complex geometries, hierarchical features, gradients, and soft-hard composites, dramatically expanding design degrees of freedom and complicating optimization. The design space becomes extremely large (e.g., 7.7 × 10^43 configurations in this study’s context), making brute-force simulations intractable and motivating fast, parallelizable surrogate models. Double-auxeticity in 2D lattices is inherently rare (<1.6% of designs), and combining it with high stiffness increases rarity. The study leverages deep learning to map from design parameters (random soft/hard assignments in lattices with varying unit-cell angles and hard-phase volume fractions) to anisotropic elastic properties, drastically accelerating evaluations and enabling rare-event discovery. The objectives are to: (i) develop accurate deep learning models trained on finite element (FE) data to predict elastic moduli and Poisson’s ratios of multi-material lattices; (ii) use these models to explore billions of designs to find rare, desirable property combinations; (iii) validate predictions experimentally via multi-material 3D printing and mechanical testing, including DIC; and (iv) extend to tiled combinations (four- and nine-tile assemblies) to further expand achievable anisotropic properties.

Prior work has demonstrated that multi-material 3D printing and rational micro-architecture design enable tailored elastic responses, including auxetic behavior and independent tuning via randomness. Foundational theories (e.g., Hashin–Shtrikman bounds) describe stiffness–Poisson’s ratio trade-offs in composites. Machine learning has achieved substantial advances in materials design and property prediction, including deep-learning-enabled metamaterial design, acceleration of auxetic metamaterial discovery, Bayesian optimization for unusual behaviors, predicting stress–strain fields in composites, and learning hidden elasticity to bypass inverse problems. However, the application of AI to design multi-material mechanical metamaterials targeting extremely rare properties (e.g., double-auxeticity with high stiffness) had not been demonstrated. This work addresses that gap by training deep neural networks on large FE-generated datasets to predict anisotropic elastic properties from microstructural design parameters, enabling rapid, large-scale exploration.

Design space and lattices: Planar lattices composed of 5×5 unit cells were considered from three unit cell angle groups: re-entrant (θ=60°, typically negative Poisson’s ratio), orthogonal (θ=90°, near-zero Poisson’s ratio), and honeycomb (θ=120°, positive Poisson’s ratio). Overall lattice dimensions (W, C), unit cell dimensions (w, c), and in-plane (t) and out-of-plane (T) thicknesses were fixed (see Supplementary Table 5). Each lattice’s struts were randomly assigned hard or soft phases to achieve hard-phase volume fractions p_h (%) ∈ {5,10,20,30,40,50,60,70,80,90,95}. Combinations of unit cells were also studied by arranging multiple tiles into four-tile and nine-tile structures. Adjacent tiles were connected by a row of hard-material struts. Finite element modeling: Custom MATLAB codes generated geometries, assigned material phases, and performed FE simulations to estimate elastic properties (E11, E22, ν12, ν21). Three-node quadratic Timoshenko beam elements with rectangular cross-sections were used, with DOFs ux, uy (translations) and uz (rotation) per node. Materials were linear elastic with ν=0.48 for both phases and Young’s moduli Es=0.60 MPa (soft) and Eh=60 MPa (hard), i.e., Eh/Es=100. Boundary conditions estimated properties in each orthogonal direction by applying 3% engineering strain and appropriate constraints: (i) for ε22 loading, top nodes displaced to achieve 3% strain while bottom nodes fully constrained; (ii) for ε11 loading, right nodes displaced to 3% strain while left nodes fully constrained. Global stiffness and load were assembled (K, F), and displacements computed (F=Kd). Stresses and moduli were calculated from reaction forces and applied strains: E11=σ11/ε11 and E22=σ22/ε22, with σ based on mean reaction force over loaded-edge nodes divided by cross-sectional area; Poisson’s ratios ν12=ε22/ε11 and ν21=ε11/ε22 from measured transverse strains. Stress analyses included von Mises stress distributions and axial strut stresses to assess stress uniformity and concentration. Deep learning models: Two neural networks were implemented in TensorFlow/Keras: (1) Single unit cell model mapping design parameters to lattice elastic properties; (2) Four-tile model mapping properties of four constituent tiles to overall four-tile properties. For the single unit cell model, inputs were 151 features: 150 binary indicators for strut phase (1 hard, 0 soft) plus the unit cell angle (encoded as θ=60°, 90°, 120°). Outputs were E11, E22, ν12, ν21 (4 total). The FE dataset included 18,150,000 lattices (training: 16,500,000; testing: 1,650,000). For the four-tile model, 90 single tiles with properties uniformly covering the attainable range were selected (their properties computed by FE). All 4-combinations with permutation reduced by symmetry yielded n=6×C(90,4)=15,331,140 samples, split 90%/10% for training/testing. Inputs were 16 features (E11, E22, ν12, ν21 of each of the four tiles), outputs were 4 overall properties (E11, E22, ν12, ν21) of the assembled structure. Training setup and hyperparameters: Both models used fully connected feedforward architectures with ReLU activations, RMSprop optimizer, learning rate 1e-4, 200 training epochs, min–max feature scaling for outputs (single unit cell) and for both inputs/outputs (four-tile). Validation sets were 20% of training data for early detection of overfitting and hyperparameter tuning. Optimized architectures: Single unit cell model hidden layers 256-128-128-64-32; Four-tile model hidden layers 256-128-64-32-16-8. Metrics included MSE and MAE. Experimental validation: Multi-material specimens were fabricated using a PolyJet printer (Stratasys Objet500 Connex3) with VeroCyan (hard, RGD841) and Agilus30 white (soft, FLX985), targeting Eh≈60 MPa and Es≈0.6 MPa (ratio≈100). A total of 10 specimens per design were tested under uniaxial tension (LLOYD LRSK, 100 N load cell; stroke rate 1 mm/min). Stress–strain curves were obtained from force/area and displacement/length; stiffness from slope. Transverse strains were measured via imaging and custom MATLAB analysis; Poisson’s ratio ν=−εtrans/εaxial. Digital image correlation (LIMESS Q400-3D-12MP with Istra4D) was used to capture full-field strain for selected single-unit-cell designs after applying a speckle pattern. Tested designs included three single-unit-cell lattices (θ=60°, 90°, 120°), four four-tile structures, and three nine-tile structures.

Computation speed and accuracy: On a workstation (Intel Core i9-8950HK, 32 GB RAM), individual FE simulations took 6.2×10^-2 ± 2.7×10^-3 to 6.5×10^-3 ± 8.2×10^-4 s depending on parallelization; deep learning (DL) predictions took 8.3×10^-2 ± 2.9×10^-3 to 1.2×10^-5 ± 1.2×10^-6 s. With 10^5 FE and 10^6 DL scripts run in parallel, per-design times were 5.0×10^-3 ± 4.9×10^-4 to 2.5×10^-3 ± 1.3×10^-4 s (FE) and 1.3×10^-5 ± 1.8×10^-7 to 2.4×10^-6 ± 1.2×10^-7 s (DL). Thus, DL reduced per-design evaluation to ~2.4 μs and enabled massive parallelization. Training performance: Single unit cell model errors improved over 200 epochs from MAE 6.6×10^-4 and MSE 1.38×10^-2 to MAE 1.05×10^-4 and MSE 6.0×10^-3; validation from MAE 4.25×10^-4 and MSE 1.19×10^-2 to MAE 1.14×10^-4 and MSE 6.38×10^-3. Four-tile model errors improved from training MSE 3.1×10^-4 and MAE 1.24×10^-2 to MSE 3.27×10^-5 and MAE 4.3×10^-3; validation to MSE 3.68×10^-5 and MAE 4.6×10^-3. Both models achieved R^2≈0.998. Predicted property ranges (1.5×10^9 random single tiles): E11 ∈ [0, 10.94] MPa; E22 ∈ [0, 0.55] MPa; ν12 ∈ [−1.24, 1.16]; ν21 ∈ [−0.53, 0.51]. Orthogonal unit cells yielded the highest E11 (stretch-dominated). Re-entrant lattices produced highly negative ν; honeycomb produced highly positive ν. For p_h up to ~80%, E11 and |ν12| were inversely correlated, consistent with composite bounds; above ~80%, a direct correlation emerged due to non-affinity effects from random hard-phase distributions. Non-affinity peaked for p_h≈75–90%. Effect of multi-material ratio and tiling: In single-unit-cell lattices, both E11 and E22 increased monotonically with p_h, while ν showed non-monotonic behavior with extrema for p_h≈60–80% (where both phases comparably influence response). For fixed ν and p_h, achievable E spanned widely; e.g., orthogonal lattices at p_h=80% exhibited up to 10.7-fold E variation with minimal ν change. The four-tile model greatly expanded attainable anisotropy, especially in direction 2 where the property region approached a near-square, allowing high E to co-occur with highly positive or negative ν. Double-auxetic stiffness rarity and enhancement: Using α=Ē11×Ē22×ν̄12×ν̄21 (normalized properties) as a combined stiffness–auxeticity metric, 0.08% (single tiles) and 0.58% (four-tile assemblies) of double-auxetic designs achieved α at least 5 standard deviations above the mean, showing four-tile combinations substantially increase access to rare double-auxetic yet stiff designs. Design combinations and rotation: Combining extreme-property unit cells demonstrated synergistic gains. Mixing re-entrant with orthogonal increased E11 of the re-entrant constituent by 75.6% (four-tile structure 1) and 91.4% (nine-tile structure 5) while preserving highly negative ν12 (reductions of only 6.6% and 9.4%). Mixing honeycomb with orthogonal increased E11 by 40% (structure 2) and 40.4% (structure 6) with minimal ν12 change (reductions 1.4% and 2.8%). Rotating tiles by 90° increased E22 and improved isotropy: in structure 3, E22 exceeded constituent E22 by 56.2% (re-entrant) and 31% (orthogonal); in structure 4, E22 exceeded constituent E22 by 4.5% (honeycomb) and 75.3% (orthogonal). FE and experiments agreed closely (e.g., four-tile structure 1 FE: E11=1.63 MPa, E22=0.22 MPa, ν12=−1.10, ν21=−0.21; experiment: E11=1.64 MPa, E22=0.38 MPa, ν12=−0.96, ν21=−0.32; similar agreement across seven multi-tile structures). Stress uniformity as an additional criterion: Among 207 designs (θ=60°) within a narrow target property window (0.15<E11<0.25 MPa, −1.1<ν12<−1, 0.01<E22<0.08 MPa, −0.45<ν21<−0.3), maximum von Mises stresses varied up to 2.5× and 6.5× under ε11=3% and ε22=3%, respectively, despite similar global properties. Two p_h=70% designs showed stark contrast: one exhibited stress concentration (potential crack initiation), while the other had uniform stresses. Comparing axial stress distributions to an equivalent homogeneous material indicated the lower-stress-riser design’s stresses lay within the homogeneous-material confidence ellipse, suggesting stress-uniformity criteria can filter designs for improved integrity.

The study addresses the challenge of efficiently exploring enormous multi-material metamaterial design spaces to find rare, high-value property combinations (e.g., double-auxeticity with high stiffness). By training highly accurate deep neural networks on extensive FE-generated datasets, the authors map complex design parameters (random phase assignments and unit-cell geometry) to anisotropic elastic properties with microseconds per-design evaluation time and excellent predictive fidelity (R^2≈0.998). This acceleration enables screening of billions of designs and supports rare-event discovery. The exploration reveals how p_h governs stiffness monotonically but Poisson’s ratio non-monotonically, with non-affinity effects driving atypical stiffness–Poisson correlations at high p_h. Tiling strategies (four- and nine-tile assemblies) significantly broaden the property space, especially in achieving high stiffness concurrently with extreme positive or negative Poisson’s ratios and in enabling double-auxetic-yet-stiff designs at higher frequency than single tiles. Rotational transformations can improve isotropy by elevating the weaker directional modulus. Importantly, the giga-scale database allows imposing new constraints like stress uniformity, revealing that designs with similar bulk properties can differ by multiples in peak local stress, which has major implications for durability and failure risk. Experimental validation via multi-material printing, mechanical testing, and DIC confirms model predictions for both single and tiled designs, reinforcing the approach’s reliability for rational metamaterial design.

Deep learning models accurately predict the anisotropic elastic properties of multi-material mechanical metamaterials while reducing per-design evaluation time to the microsecond scale and enabling highly parallel exploration of vast design spaces (10^10–10^20 candidates within reach). This capability facilitates rational discovery of rare-event designs, such as double-auxetic lattices that are simultaneously stiff, and supports additional selection criteria like stress-uniformity for improved structural integrity. The work demonstrates that combining and rotating tiles with extreme properties can significantly enlarge the attainable property space and tune anisotropy. Future extensions include expanding input parameters (e.g., allowing variable hard/soft modulus ratios beyond Eh/Es=100), exploring broader geometric families beyond θ=60°, 90°, 120° and fixed 5×5 unit cells, integrating fracture and fatigue performance metrics into the design objective, and leveraging the ML surrogate to guide inverse design and multi-objective optimization.

The simulations assumed linear elastic behavior with identical Poisson’s ratio (ν=0.48) for both phases and a fixed stiffness ratio Eh/Es=100; the authors note that varying this ratio can affect non-affinity and overall properties and should be studied further. Geometrical parameters (overall dimensions, unit-cell dimensions, thicknesses) and lattice size (5×5 unit cells) were fixed, and only three unit-cell angles (60°, 90°, 120°) were considered, which may limit generalizability to other architectures. Properties were evaluated under 3% strain with beam-element FE idealizations, which may not capture nonlinearities or failure mechanisms. Stress-uniformity analysis, while informative, was demonstrated on selected subsets and boundary conditions. Data and code are available upon reasonable request rather than fully open, potentially limiting immediate reproducibility.