Robust topological designs for extreme metamaterial micro-structures

We demonstrate that the consideration of material uncertainty can dramatically impact the optimal topological micro-structural configuration of mechanical metamaterials. The robust optimization problem is formulated in such a way that it facilitates the emergence of extreme mechanical properties of metamaterials. The algorithm is based on the bi-directional evolutionary topology optimization and energy-based homogenization approach. To simulate additive manufacturing uncertainty, combinations of spatial variation of the elastic modulus and/or, parametric variation of the Poisson’s ratio at the unit cell level are considered. Computationally parallel Monte Carlo simulations are performed to quantify the effect of input material uncertainty to the mechanical properties of interest. Results are shown for four configurations of extreme mechanical properties: (1) maximum bulk modulus (2) maximum shear modulus (3) minimum negative Poisson’s ratio (auxetic metamaterial) and (4) maximum equivalent elastic modulus. The study illustrates the importance of considering uncertainty for topology optimization of metamaterials with extreme mechanical performance. The results reveal that robust design leads to improvement in terms of (1) optimal mean performance (2) least sensitive design, and (3) elastic properties of the metamaterials compared to the corresponding deterministic design. Many interesting topological patterns have been obtained for guiding the extreme material robust design.

Metastructures are metamaterial induced concepts implanted in structural design. The exploration of metamaterials was heavily steered by the path-breaking investigation of electromagnetic metamaterials, which exhibited properties like, negative permittivity and permeability 1 . This paved the way for mechanical metamaterials 2 . These metamaterials are engineered to derive their fundamental mechanical properties from the geometry of their structural building blocks (periodicity or translational symmetry), instead of the constituting materials 3,4 . Mechanical metamaterials have gained wide popularity in engineering applications due to their superior properties compared to conventional materials found in nature. It has been found that by varying one or more material constants such as, bulk modulus, shear modulus, Poisson's ratio and elastic modulus, extraordinary materials such as, ultra light weight and high strength aerospace shells, functionally graded electromagnetic sensors and energy absorbing dampers can be secured 5,6 .
For exploring such exciting possibilities of discovering new materials with extreme properties, determining their optimal configuration at the very conceptual micro-structural design stage is indispensable 7 . A recent works with multi-material lattices 8 , pre-stressed materials 9 and piezo-embedded damped lattice metamaterials 10,11 show some ways of achieving extreme homogeneous properties. In this scenario, topology optimization (TO) is well-known to be an effective computational analysis tool which renders greater design freedom and yields the best material layout given a prescribed design domain and boundary condition 12 . Depending on the algorithm's performance and versatility, various successful TO methods have emerged to the present 13 , few popular ones being the solid isotropic material with penalization (SIMP), evolutionary structural optimization (ESO), level set method (LSM), and moving morphable components (MMC). In fact, the popularity of TO has even led to its integration with relatively newer numerical discretization schemes (compared to finite element method) such as, isogeometric analysis 14 . A simple illustration of multi-scale design and topology optimization of the microstructure of metamaterial is presented in Fig. 1a and b, respectively. An illustrative visual tool guide of the key aspects of the proposed design paradigm. The focus is on optimizing the topology of the periodic unit cell as shown in (a) and determine robust micro-structural designs against uncertainties. Note that the objective of the deterministic topological design in (b) is to maximize the bulk modulus (refer Table 1 for details). The negative ordinate indicates the maximization problem. The 3D computer aided design (CAD) geometries of the robust micro-structural configurations shown in (c) are obtained using Autodesk Fusion 360 software (version 2.0.9930) (https:// www. autod esk. co. uk/ campa igns/ educa tion/ fusion-360). These geometries can be directly exported for 3D printing. The four different geometries shown here correspond to robust topologies exhibiting extreme mechanical properties (for example, maximum bulk modulus, maximum shear modulus, minimum Poisson's ratio and maximum equivalent elastic modulus).

Scientific Reports
| (2021) 11:15221 | https://doi.org/10.1038/s41598-021-94520-x www.nature.com/scientificreports/ The success of TO is clearly evident from its extensive applications for material micro-structural design in recent years. In this regard, few note-worthy literatures have been discussed as follows. Gibiansky and Sigmund 15 investigated multi-phased elastic composites and obtained topological design corresponding to extreme bulk and shear modulus. For obtaining maximum bulk and shear modulus, the optimal topology of cellular microstructures was determined using ESO 16 . A compact code in MATLAB based on the SIMP and energy-based homogenization method (EBHM) was developed for extreme material multi-scale structural design 17 . Clausen et al. 18 designed and fabricated 3D auxetic material micro-structures undergoing large deformations. Long et al. 19 performed topology optimization to maximize the effective Young's modulus, so as to obtain the optimal distribution of 3D material micro-structures whose constituent phases consist of non-identical Poisson's ratios. A novel TO method was presented based on the independent point-wise density interpolation to obtain a bimaterial chiral metamaterial 20 . Chen and Huang 21 designed 3D chiral metamaterials based on the couple-stress homogenization and maximized the chiralty (the coupling effects between tension/compression and torsion) of the cellular structure. Xu et al. 14 achieved designs of ultra-lightweight architected materials along with extremal properties such as, maximum bulk and shear modulus by using isogeometric TO. Petal-shaped auxetic metamaterials were designed where a back-propagation neural network was coupled with isogeometric TO to evaluate their macroscopic equivalent properties 22 . Zheng et al. 23 proposed two numerical schemes within the evolutionary optimization framework to obtain 2D auxetic metamaterials with favourable characteristics, and experimentally validated their special properties. Ye et al. 24 proposed a method for the simultaneous adaptive design of gradually stiffer mechanical metamaterials along with auxetic property.
Although the above published articles related to micro-structural TO explore multiple potent research prospects to achieve extreme material properties, they are restricted to the deterministic optimization framework. In this regard, the literature is found to be scarce in studies addressing uncertainties in manufacturing processes. The following few works highlight the importance of considering manufacturing variability and their effect on the structural design. Most nano, micro and macro structures fabricated using e-beam lithography, etching and milling processes, respectively are vulnerable to manufacturing uncertainties. Even the structures meticulously designed and optimized may experience inferior performance or in an extreme scenario, lose their functionality due to the wear of machining tools, under or over-etching, or malcalibrated e-beam equipment. To address few of the above aspects, robust topology optimization (RTO) for structural design was presented accounting for manufacturing errors in [25][26][27] . Specifically, the effect of over-etching (erosion) and under-etching (dilation) was modelled for structures produced by milling or etching, which may cause the structural parts to become thinner or thicker than intended. For the works 25,26 , the problem was posed as a worst case design formulation and the error effects were simulated by a Heaviside projection threshold. While they considered only uniform manufacturing errors (i.e. constant in magnitude over the entire design domain), the method was extended to account for non-uniform manufacturing errors (i.e. with spatially varying magnitude) 27 . For doing so, a probabilistic approach was followed where the projection threshold was represented by a (non-Gaussian) random field. Jansen et al. 28 developed an RTO framework which incorporated the spatially varying geometric imperfections due to misalignment and misplacement of material often encountered in slender (civil engineering) structures. These imperfections should be accounted for in accordance to the provision recommended in the Eurocode for design of steel structures and by the Joint Committee on Structural Safety. Note that the misplacement errors cause a perturbation in the spatial location of the material, whereas the errors due to etching, add or remove material at the structural surface. An ESO based RTO framework was developed for multi-material structures considering interval loading uncertainty in 29 . Efficiency was enhanced by (1) orthogonal decomposition of the loading uncertainty which decoupled the stochastic problem from the optimization loop and (2) sensitivity analysis. Material uncertainty modelled as imprecise probability was integrated with the multi-scale concurrent ESO based TO framework in 30 . The RTO methodology adopted the type I hybrid random interval model. An ESO based RTO framework was developed for multiscale 2D and 3D structures considering multiple random loading conditions in 31 . Decoupled sensitivity analysis was observed to improve the computational efficiency. An ESO based RTO methodology was proposed for actuator-coupled structures considering hybrid uncertainties in 32 . The hybrid random interval model simulated the stochastic mechanical and piezoelectric parameters and perturbation analysis was employed to evaluate the statistical quantities.
While the above works simulated manufacturing uncertainties numerically and investigated their effect on the structural topology, the following studies have carried out experimental analysis to probe manufacturing uncertainties. It has been shown that the emission rate of a 3D printer based on fused deposition modeling principle increased significantly with the increase in the extrusion temperature 33 . Melenka et al. 34 found that the per cent infill has a significant effect on the longitudinal elastic modulus and ultimate strength of the test specimens, whereas print orientation and layer thickness fail to achieve significance. Dimensional analysis of test specimens also showed that the test specimen varied significantly from the nominal print dimensions. The manufacturing tolerances of metamaterial samples produced by a selective laser sintering process were obtained to simulate the manufacturing uncertainty in mass-produced industrial applications 35 . It was found that even small levels of variability, given by less than 1 % for the mass and less than 3 % for the elastic modulus, has a significant effect on the overall vibration attenuation performance.
While the first of the above three paragraphs discusses works on deterministic TO for extreme metamaterial design, the second paragraph illustrates the numerical simulation of manufacturing uncertainty and robust topological design. Finally, the above paragraph highlights a few instances of experimental exploration of manufacturing uncertainties. Motivated by the need to consider manufacturing uncertainties for metamaterial design, this work attempts to bridge the above points. Thus, we integrate these aspects to present an RTO formulation specifically for extreme metamaterial design by minimizing the effect of potential material uncertainties. To the best of the authors' knowledge, this is one of the first applications of robust micro-structural design for obtaining extreme mechanical properties of metamaterials. Few instances of ready-to-print robust micro-structural Scientific Reports | (2021) 11:15221 | https://doi.org/10.1038/s41598-021-94520-x www.nature.com/scientificreports/ topologies of extreme metamaterials achieved as a part of this study have been presented in Fig. 1c for illustrative purpose and the interest of readers.

Methodology
The work integrates the material uncertainty model (simulated by in-house developed codes) with deterministic topology optimization 36 to develop a parallel Monte Carlo based robust topology optimization framework. A flow-diagram of the proposed RTO framework is presented in Fig. 2.The main intent is to illustrate the effect of material uncertainties on the micro-structural topologies and derive robust configurations for extreme material design. As shown in Fig. 2, the entire methodology comprises of three steps: (1) the user-defined input module where the system and the optimization parameters are initialized, (2) the analysis module where the expensive computations take place involving uncertainty quantification nested within the optimization loop. To reduce the computational cost, a parallel version of Monte Carlo simulation is employed to simulate the randomness within each optimization iteration. The details of this step regarding the optimization algorithm, finite element model, energy-based homogenization and sensitivity computations are presented subsequently. (3) The post-processing module is where the resulting robust micro-structural topologies with extreme properties are exported to a CAD software for 3D printing.
Homogenization. Considering linear elastic behaviour, the equivalent constitutive behaviour of periodic structures can be determined by homogenization technique. This work employs the energy-based homogenization approach to predict the macroscopic equivalent elastic properties of micro-structures 36 . The approach is based upon energy conservation with respect to stress and strain. The homogenized stiffness tensor E H ijkl can be expressed in terms of mutual energies as where Y represents the base cell, E pqrs is the elasticity tensor in index notation, ε A(ij) pq and ε A(kl) rs denote the superimposed strain fields. When the base cell is discretized into N finite elements, Eq. (1) can be approximated as  (3) post-processing unit for additive manufacturing, are bordered by dotted lines which constitute the approach. A square-shaped void region at the centre of the 150 × 150 unit cell is adopted as the initial design (as shown in the right side of the user-defined module block). The robust optimal topology obtained (corresponding to maximum shear modulus subjected to the spatial variability of the elastic modulus (refer Table 2)) is shown as a 20 × 30 periodic arrangement of the unit cell in the post-processing module. The 3D fabrication-ready object is obtained after extruding and material rendering in Autodesk Fusion 360 software (version 2.0.9930) (https:// www. autod esk. co. uk/ campa igns/ educa tion/ fusion-360). The abbreviations RTO and UQ represent robust topology optimization and uncertainty quantification, respectively. are the element displacement solutions corresponding to the unit test strain fields ε 0(kl) and k e denotes the element stiffness matrix. It is to be noted that for 2-D problems, 11 ≡ 1 , 22 ≡ 2 , and 12 ≡ 3 . Thus, Eq. (2) can be concisely represented as where Q ij is given by e are the element mutual energies and can be expressed as Periodic boundary conditions. The displacements on a pair of opposite boundaries of a 2-D base cell where notations k− and k+ represent the pair of two opposite parallel boundary surfaces perpendicular to the kth direction. One can remove the unknown periodic fluctuation term u * i by using the difference between the displacements as, For a given unit cell, the quantity (y k+ j − y k− j ) is constant. Thus, the right side of Eq. (7) is a constant for a specified ε 0 ij . This boundary conditions can be incorporated within the finite element model by restraining the corresponding pairs of nodal displacements. Note that this boundary configuration satisfies the periodicity and continuity conditions in terms of displacement and stress employing displacement-based finite element analysis 37 .
Optimization approach. The unit cell is discretized into N finite elements, resulting into equal number of density design variables ρ ∈ R N . The element level elastic modulus can be defined by the modified SIMP method 38 as where E min denotes the elastic modulus of the Ersatz material, an approximation for void material to prevent singularity of the stiffness matrix, 0 < ρ e < 1 , limits corresponding to Ersatz and solid materials, p is a penalization factor for driving the density towards the black and white solution and E 0 denotes the elastic modulus of solid material. The deterministic optimization problem can be expressed as where the objective function c(E H ijkl ) is a function of the homogenized stiffness tensors, K is the global stiffness matrix, u A(kl) and f kl are the global displacement and external force vectors of the test case (kl), respectively, d, v e and v denote the spatial dimension, element volume and upper limit of the volume fraction, respectively.
The motive of this work is to determine the extreme properties of material micro-structures, such as maximum bulk modulus, maximum shear modulus, minimum Poisson's ratio and maximum equivalent elastic modulus. The corresponding objective functions can be defined as • Design of maximum bulk modulus: • Design of maximum shear modulus: • Design of maximum equivalent elastic modulus: Note that the topological design of minimum and negative Poisson's ratio (E 1122 /E 1111 ) has always been a difficult issue. It is observed from the literature that auxetic properties can be obtained by imposing additional constraints on isotropy or bulk modulus 36 . An equivalent expression in the form of Eq. (12) is used for this purpose and adopted from 23 . Selection of the term δ in Eq. (12) is discussed later in the numerical illustration section.
The solution scheme for solving the equilibrium equation in Eq. (9) adopted here can be found in Xia and Breitkopf 36 . After obtaining the finite element displacement solution, the optimization problem in Eq. (9) is solved by using optimality criteria method. The heuristic updating criterion 13 is expressed as where m is the step size, η is the damping coefficient and B e can be determined by using the following optimality condition where is the Lagrange multiplier. It imposes the satisfaction of the material volume fraction constraint. The numerator ∂c ∂ρ e which is the sensitivity of the objective function can be determined by using the adjoint method as where k 0 represent the element stiffness matrix with unit elastic modulus. For a uniform mesh, the element volume v e = 1 and hence, ∂V ∂ρ e = 1 . To ensure convergence of the optimization problem Eq. (9), filtering techniques are used which are illustrated next.

Filtering.
To eliminate the formation of a checkerboard pattern and resolve the mesh-dependent issue, a filter is applied either to the sensitivities or the densities 38,39 . The sensitivity based filtering modifies the sensitivities ∂c ∂ρ e (Eq. (16) as, where γ = 10 −3 is a small positive number incorporated to avoid the denominator becoming zero. H ei is a weight term defined in Eq. (18). N e represents the set of elements i whose centre-to-centre distance to element e, �(e, i) , is less than the filter radius r min . The size of the filter radius can be used to control the minimum size of the emerging features in the design domain.
The density based filtering converts the original densities to the filtered densities ρ e as, The sensitivities with respect to ρ j can be obtained as, Note that ∂c ∂ρ e can be evaluated by Eq. (16) upon replacing ρ e with ρ e .
Robust topology optimization. The objective function considered for the robust topology optimization is (11) c = −E 1212 ρ new ∂c ∂ρ e = 1 max(γ , ρ e ) i∈N e H ei i∈N e H ei ρ i ∂c ∂ρ i

Numerical illustration: results and discussion
Implementation details. In this work, the design domain is considered to be square and discretized into square plane stress elements. The volume fraction is taken to be 0.5. The unit cell is discretized into 3 × 3 elements. The influence of an initial design can affect the final optimal topology significantly. For this work, a square shaped void region at the center of the unit cell is adopted as the initial design 40 as shown in Fig. 2. The small void is generally inserted in the initial design to stimulate the evolutionary optimization process.
It has been reported in multiple studies that topology optimization for the design of extreme material properties is prone to multiple local minima. The final optimal solution is found to be extremely sensitive to the initial topology, shape of the unit cell, filter radius (r min ) , penalization factor (p), filtering approach (sensitivity filtering or density filtering) and others. Therefore, for selecting the appropriate parameter values, few rules of thumb have been followed. The optimization problem Eq. (9) is non-convex for p ≥ 1 . Although the high values of penalization factor will result in more distinct topologies, the algorithm is likely to get trapped in a local minima 17 . A smaller filter radius value generally results in a better solution (in terms of a detailed micro-structure) as higher frequency details are visible through low-pass filter. Compared to sensitivity filtering, density filtering is observed to yield less sensitive designs and hence preferable for material micro-structure design.
Parametric study has been carried out by considering various combinations of stochastic material parameter models such as, spatial variation via random field and random parametric models. The stochastic models assumed are in line with the literature (discussed in the introduction section) where spatially varying and parametric models simulate non-homogenous and uniform manufacturing errors, respectively. For stochastic simulations, the elastic modulus and Poisson's ratio are considered to be random and follow log-normal distribution with 5 % variability. It is practical to assume that all realizations corresponding to these random parameters will be positive and therefore, they are assumed to be log-normally distributed. In doing so, 1000 MCS is performed within www.nature.com/scientificreports/ every optimization iteration to simulate the material variability. This has been carried out in a cost-effective manner using parallel computing (employing all 8 cores of an Intel Xeon processor). Four case studies have been undertaken to obtain extreme properties of material micro-structures, such as maximum bulk modulus, maximum shear modulus, minimum negative Poisson's ratio (auxetic metamaterials) and maximum equivalent elastic modulus. The initial design is the same for all the following cases studied (except type III). The nominal values of the elastic modulus and Poisson's ratio at the element level are adopted as 1 and 0.3, respectively. For the stochastic cases, their mean values are same as that of the nominal values. The length scale used for discretizing the random field model is taken as 0.2. The results are represented by reporting the optimal unit cell, 3 × 4 periodic arrangement, equivalent elastic matrices and robust optimal solutions for each of the extreme material configurations. Table 1 illustrates the optimal solutions for the case of type I. To obtain these results, density based filtering and the following parameter values are considered: p = 5 and, r min = 2. It can be observed from Table 1 that uncertainty affects both the topological configuration and the equivalent elastic matrix. The stiffness terms Q 11 and Q 22 are observed to increase for all types of uncertainty models compared to the deterministic case. Specifically, the greatest rise is observed for the random field cases. The term Q 33 decreases due to uncertainty. The off-diagonal term Q 12 reduces for all uncertainty models and the deterministic model achieves the maximum value. Interestingly, for the hybrid input uncertainty modelling case, the resulting structure shows non-zero (negative) Q 13 and Q 23 terms. This indicates coupling of axial and shear deformations/ stress states. Although there is a slight difference between the terms Q 11 and Q 22 for the deterministic and ran- Table 1. Optimal solutions for maximum bulk modulus (Type I). Results from four case studies are shown. They include, (1) deterministic design, (2) robust design considering parametric uncertainty models for elastic modulus and Poisson's ratio, (3) robust design considering random field model for elastic modulus and (4) robust design considering a combination of random field and parametric uncertainty models for the elastic modulus and Poisson's ratio, respectively. Note that the terms of the equivalent elastic matrix (given by Eq. (4)) represent their mean values. The standard deviation of the matrix terms is consistent with the level of input uncertainty and is observed to be within 5% bounds. The objective functions of the deterministic and robust optimization have been evaluated by Eqs. (10) and (21) www.nature.com/scientificreports/ dom field models, the difference is significant for the hybrid uncertainty model, exhibiting strong anisotropic behaviour. In terms of mean value of the objective function, the random field cases lead to higher value of bulk modulus compared to the deterministic and parametric uncertainty cases. Also, the random field cases result in a more robust topological design indicated by lower values of the standard deviation of objective function. Table 2 illustrates the optimal solutions for the case of type II.

Maximum shear modulus (Type II).
To obtain these results, density based filtering and the following parameter values are considered: p = 5 and, r min = 2.It can be observed from Table 2 that uncertainty affects both the topological configuration and the equivalent elastic matrix. The stiffness terms (diagonal entries) Q 11 , Q 22 and Q 33 and the off-diagonal term Q 12 are observed to increase for the random field and hybrid input uncertainty models compared to the deterministic case. These terms achieve lower values for the parametric uncertainty model compared to the deterministic case. Interestingly, for the random field and hybrid input uncertainty models, the resulting structure show (1) nonzero values of Q 13 and Q 23 terms and (2) Q 11 = Q 22 . (1) indicates coupling of axial and shear deformations/stress states and (2) illustrates anisotropic behaviour. In terms of mean value of the objective function, the random field and hybrid input uncertainty models lead to higher value of shear modulus compared to the deterministic and parametric uncertainty cases. The parametric stochastic model results in the minimum value of shear modulus. Also, the random field and hybrid input uncertainty models result in a more robust topological design indicated by lower values of the standard deviation of objective function.  (4) robust design considering a combination of random field and parametric uncertainty models for the elastic modulus and Poisson's ratio, respectively. Note that the terms of the equivalent elastic matrix (given by Eq. (4)) represent their mean values. The standard deviation of the matrix terms is consistent with the level of input uncertainty and is observed to be within 5 % bounds. The objective functions of the deterministic and robust optimization have been evaluated by Eqs. (11) and (21) Type III). Density based filtering and the following parameter values are considered: p = 5 and, r min = 5 , for the type III case. The unit cell is discretized into 200 elements along both the horizontal and vertical directions. As the term δ (refer Eq. (12)) causes significant sensitivity to the topological design of auxetic metamaterials 23 , a parametric study is performed in Table 3 to investigate the effect of δ on the resulting topologies with minimum Poisson's ratio.It can be observed from Table 3 that out of the adopted δ values, only δ = 0.5 and 1 lead to negative mean Poisson's ratio corresponding to the deterministic and each of the stochastic models. δ = 0.5 has been selected as lower Poisson's ratios are obtained compared to the case δ = 1. The corresponding row has been highlighted to indicate the best and the most consistent performance in Table 3 for clarity. Table 4 illustrates the optimal solutions for the case of type III, corresponding to δ = 0.5.It can be observed from Table 4 that uncertainty affects both the topological configuration and the equivalent elastic matrix. Although the overall resulting topologies from the deterministic and parametric uncertainty models are very similar with fine differences at the boundary. The stiffness terms Q 11 and Q 22 are observed to increase for the random field and hybrid uncertainty models and decrease for parametric uncertainty compared to the deterministic case. The stiffness term Q 33 is observed to increase for all stochastic models compared to the deterministic case. The off-diagonal term Q 12 is negative for all the cases reported and observed to increase for the parametric stochastic model and decrease for the random field and hybrid uncertainty cases than that of the deterministic one. Interestingly, for the random field and hybrid input uncertainty models: (1) the terms Q 13 and Q 23 are non-zero (negative), which illustrates coupling of axial and shear deformations/stress states and (2) Q 11 = Q 22 , which indicates anisotropic behaviour. The random field and hybrid input uncertainty models lead to lower values of the mean objective function compared to the deterministic and parametric uncertainty cases. Specifically, the parametric and hybrid stochastic models result in the maximum and minimum value of Poisson's ratio, respectively. The deterministic, random field and hybrid input uncertainty models result in the same level of robustness of the topological design indicated by the standard deviation of the objective function. The parametric uncertainty model achieves the most robust design by a small margin. The results in Table 4 illustrate that negative Poisson's ratio (given by Q 12 /Q 11 ) has been achieved for all the case studies undertaken. This is worth mentioning as the design of auxetic metamaterials by topology optimization is a challenging task on its own and even more difficult with uncertainties.
Maximum equivalent elastic modulus (Type IV). Table 5 illustrates the optimal solutions for the case of type IV. To obtain these results, density based filtering and the following parameter values are considered: p = 5 and, r min = 2.It can be observed from Table 5 that uncertainty affects both the topological configuration and the equivalent elastic matrix. The resulting topologies obtained from the random field and hybrid random input models are the same. The stiffness terms (diagonal entries) Q 11 , Q 22 and Q 33 are observed to increase for the random field compared to the deterministic case. The stiffness terms Q 11 and Q 33 are observed to increase for the hybrid random input model compared to the deterministic case. These terms achieve lower values for the parametric uncertainty model compared to the deterministic case. The off-diagonal term Q 12 is observed to be higher for all the stochastic models than that of the deterministic one. Interestingly, for the hybrid input www.nature.com/scientificreports/ uncertainty models, the resulting structure shows anisotropic property indicated by Q 11 = Q 22 and illustrates coupling of axial and shear deformations/stress states indicated by non-zero Q 13 term. In terms of mean value of the objective function, the random field and hybrid input uncertainty models lead to higher value of equivalent elastic modulus compared to the deterministic and parametric uncertainty cases. The parametric stochastic model results in the minimum value of equivalent elastic modulus. Also, the random field and hybrid input uncertainty models result in a more robust topological design indicated by lower values of the standard deviation of objective function.
Conclusions. This section summarizes the results obtained from the numerical investigation. Few critical points have been enumerated below: • It has been observed that the material uncertainties affect the resulting micro-structural topologies and mechanical performance significantly. • Out of the input stochastic material models, the random field model followed by the hybrid one lead to desirable performance in terms of (1) optimal mean performance, (2) most robust performance and (3) elastic property, indicated by the mean value of objective function, standard deviation of objective function and stiffness terms of the equivalent elastic matrix, respectively, for types I, II and IV. • For type III, the hybrid uncertainty model followed by the random field leads to desirable output performance (based upon combination of the above points (1)- (3)). This point along with the above are indicative  (4) robust design considering a combination of random field and parametric uncertainty models for the elastic modulus and Poisson's ratio. Note that the terms of the equivalent elastic matrix (given by Eq. (4)) represent their mean values. The standard deviation of the matrix terms is consistent with the level of input uncertainty and is observed to be within 5 % bounds. The objective functions of the deterministic and robust optimization have been evaluated by Eqs. (12) and (21) www.nature.com/scientificreports/ of the degree of sensitivity of input stochastic models on the optimal micro-structural topology. Thus, these observations have a practical relevance as one needs to be more cautious in fabricating the metastructures which follow a more sensitive input stochastic model. • The robust design leads to improved material properties compared to the deterministic design.
• Multiple interesting topological patterns have been found for guiding the robust metamaterial design. Importantly, all of the resulting topologies have tessellating property, which make them feasible for 3D printing. Thus, there can be a possible extension of the present computational work as the resulting robust topological configuration of materials with extreme properties are practically realizable. • In addition to the few resulting micro-structural configurations showing isotropic properties, interestingly, the others yielded configurations which showed anisotropic property. Some cases led to positive and negative coupling of axial and shear deformations/stress states. These phenomenon were primarily observed for the random field and hybrid uncertainty cases. This observation may be attributed to the fact that the spatial variability (inhomogeneity) in the micro-structure induces anisotropy and coupling phenomena in the material. Although this requires further investigation and achieving this was unintentional but will prove to be useful for designing material micro-structures with dual/simultaneous properties. • Considering the fact that the design of auxetic metamaterials by topology optimization is a challenging task on its own and even becomes more difficult to obtain robust designs with micro-structural uncertainties, it is worth noting that negative Poisson's ratio was achieved for all the undertaken cases of type III using the proposed framework by adopting appropriate parameters (Table 4). However, one should be careful in choosing the parameters considering their high sensitivity on the resulting topologies (for example, δ of Eq. (12)) Table 5. Optimal solutions for maximum equivalent elastic modulus (Type IV). Results from four case studies are shown. They include, (1) deterministic design, (2) robust design considering parametric uncertainty models for elastic modulus and Poisson's ratio, (3) robust design considering random field model for elastic modulus and (4) robust design considering a combination of random field and parametric uncertainty models for the elastic modulus and Poisson's ratio, respectively. Note that the terms of the equivalent elastic matrix (given by Eq. (4)) represent their mean values. The standard deviation of the matrix terms is consistent with the level of input uncertainty and is observed to be within 5 % bounds. The objective functions of the deterministic and robust optimization have been evaluated by Eqs. (13) and (21), respectively.
Scientific Reports | (2021) 11:15221 | https://doi.org/10.1038/s41598-021-94520-x www.nature.com/scientificreports/ as evident from Table 3. It is also recommended to thoroughly scrutinize the existing equivalent expressions for minimizing the Poisson's ratio and is highly sensitive to the problem type or application in hand. • The effect of uncertainties on the material micro-structures is so high that even considering the accuracy and high precision of additive manufacturing techniques, accounting nominal amount of uncertainty for material design is recommended, at least in the conceptual design stage before the manufacturing. In fact, there is nothing at stake, as it is shown in this study that robustness improves the design in terms of mechanical performance.

Data availability
All relevant data are available from the corresponding author upon reasonable request, and/or are included within the main part and Supplementary Information.