How does the initial cell configuration influence the final topology in a metamaterial generation process?

Abstract

This research aims to evaluate the impact of the initial cell configuration and the limit volume fraction on the generation of mechanical metamaterial cells. The procedure was developed using a methodology based on preprocessing, processing, and post processing to facilitate fast exploration of metamaterial cell space design. The initial cell consisted of a square or cube of material with a central circular or spherical void diameters of 10 mm, 20 mm, and 30 mm. Additionally, the generation process employed three volume fractions limits (30%, 40%, and 50%) and eleven objective functions. These functions intended to generate cells that maximize stiffness in one or multiple directions and cells with maximum compressibility or shear modulus. Some of the obtained cells with tailored mechanical properties exhibit novel geometrical configurations. The results highlight volume fraction as the most significant factor in the generation process, with well-defined metamaterial cells produced using an initial volume fraction of 50% and low to medium void diameters.

Keywords
Mechanical metamaterials; topological optimization; numerical homogenization; cell configuration

Graphical Abstract

1 INTRODUCTION

The metamaterials have gained significant attention in the scientific community due to their unique properties and versatile nature (Schittny et al., 2013; Martin et al., 2013; Kadic et al., 2019). These artificial cellular structures exhibit behavior that deviates from their constituent materials, making them highly advantageous for applications in passive seismic protection (Guevara-Corzo et al., 2022; Martin et al., 2013), non-reciprocal sound propagation (Kadic et al., 2019), camouflage of objects from incident acoustic energy (Cummer et al, 2016), vibration isolation (Ishida et al., 2017; Ungureanu et al., 2016; Miniaci, et al., 2021), unidirectional motion guidance (Bertoldi et al., 2017) and bioengineering and biomedical engineering (Yu et al.,2018). In recent years, Kadic et al. (2012), Fraternali et al. (2018), and Fabbrocino et al. (2015) have successfully applied pentamode metamaterials to structural vibration isolation. Zheng et al. (2016), present scalable metallic mechanical metamaterials that simultaneously exhibit high strength and low density. The approach to design these multiscale metamaterials, was based on assembling microscale filaments following a trajectory defined by a larger length scale structure. Zheng et al. (2014) proposed a material category that maintains nearly constant stiffness even at ultra-low densities. Jenett et al. (2020) introduced a module-based approach for manufacturing heterogeneous metamaterials that exhibit properties comparable to previously reported metamaterials, with their innovative scalable system streamlining the production process. Finally, Zhang et al. (2022) proposed a 3D metamaterial cell consisting of a hollow rhombic dodecahedron and six cylindrical tubes with shock resistance and vibration isolation capacity.

In this context, most studies (Martin et al., 2013; Kadic et al, 2019; Cummer et al, 2016; Yu et al, 2018; Kadic et al, 2012; Fraternali et al, 2018; Zheng et al., 2014; Zhang et al, 2022) have employed pre-defined cell structures. In these approaches, the designer's experience or the use of heuristic procedures was decisive in generating the metamaterial cell topology. However, this has resulted in works focusing more on adapting existing cells rather than generating new ones through exhaustive and systematic exploration. Using inverse design, Sigmund et al. (Sigmund O. , 1994; Sigmund O. , 1995; Bendsøe and Sigmund, 2002) looked at developing cellular networks with extreme properties. In another relevant study, Gibiansky and Sigmund (2000) proposed analytical and numerical methods to obtain two-dimensional, three-phase materials exhibiting extreme conductivity. Milton and Cherkaev conducted seminal research that led to the creation of the theory of analytical materials (Milton and Cherkaev, 1995). Moreover, Milton (2016) provided the first rigorous verification of the existence of auxetic materials without holes or sliding surfaces. He also highlighted the feasibility and potential of metasurfaces and planar lattice metamaterials for practical applications.

Today, the approach to metamaterial design has shifted from adapting or fine-tuning existing structures (cells) to a more pragmatic direction that integrates topology optimization with cutting-edge additive manufacturing. In recent years, Andreassen et al. (Andreassen and Andreasen, 2013; Andreassen et al., 2014) and Chatterjee et al. (2021) have developed and applied robust topological optimization techniques, mixed with numerical homogenization to obtain mechanical metamaterials. In addition, Kollmann et al. (2020) demonstrated the performance of convolutional neural networks to get the optimal design of metamaterials with maximum compressibility or shear modulus. Using couple-stress homogenization techniques, Chen and Huang (2019) reported the obtention of a new chiral metamaterial. The work of Wu et al. (2021) and Li et al. (2021) provide examples of topology optimization for the creation of pentamodal cells using the ground structure technique. The former case employed topological optimization, whereas the latter used a multivariable optimization approach with a genetic algorithm. Although the results reported by Wu et al. (2021) and Li et al. (2021) differ from the original pentamodal cell proposed by Milton and Cherkaev (1995), both authors achieved 3D petamodal lattices from a ground structure evolution via genetic algorithms.

A study by Huang et al. (2011), focused on generating metamaterials used the BESO approach to generate cells that maximized shear modulus and bulk modulus by varying the volume fraction. Although the examples were limited to a single base structure, reliable results were obtained. To complement or improve upon the shortcomings of current methods, Ai and Gao (2019) used the parametric level set method (PLSM) in combination with the Meshfree method to generate cells with maxima in shear modulus, bulk modulus, and auxetic cells. This work extended the processes to cells with two phases. Finally, Wu et al. (2017), proposed and used a robust topological optimization (RTO), an improved version of PLSM, in the generation of auxetic structures. In this context, current generative procedures focus on formulating topological optimization problems employing search methods of local nature. It is evident that, the initial cell configuration (the point to begin the search) and the parameters used in each optimization algorithm play a highly influential role in the final obtained topology.

In these circumstances, research on the generative cell process often overlooks several influential factors. Variations in the initial cell configuration, limiting volume fractions, approximation models, mutation rates, and the choice of different algorithms can all significantly impact the development of a hands-on metamaterial. In this context, this research aims to investigate the generative process of mechanical metamaterials. It focuses on the influence of the initial cell configuration on the generation of elastic mechanical metamaterials. Following a similar approach taken by Andreassen et al. (Andreassen and Andreasen, 2013; Andreassen et al., 2014), Wu et al. (2021) and Li et al. (2021), the generative process developed was based on topological optimization and numerical homogenization for 2D and 3D solids in linear elasticity. The base solid used (as initial cell configuration) was a square (2D) or cubic (3D) element of 50 mm side with an empty circular (2D) or spherical (3D) void inside. Void diameter variations took values of 10 mm, 20 mm, and 30 mm.

In addition, the optimization process was subjected to volume fractions limits of 30%, 40%, and 50% of the initial solid; lastly, eleven objective functions were studied. The optimization results were post-processed and filtered to generate manufacturable cells without intermediate densities (i.e., more suitable for additive manufacturing). In addition, for the 3D cases, due to the low voxel density, the surface was smoothed to ensure a better-defined cell topology. Finally, the metamaterial cells resulting from the generative process were represented in a simplified model and then analyzed using Ansys® (2023) to validate the outcomes. This paper is organized as follows: Section 2 presents the generation strategy used, Section 3 discusses the topological optimization strategy, Section 4 focuses on the numerical homogenization process, Section 5 describes the initial cell configuration (base solid) and the objective functions employed, Section 6 presents the results of the topological optimization and the developed idealization, and finally, Section 7 provides the conclusions of this research and possible future works.

2 MECHANICAL METAMATERIAL CELL GENERATION STRATEGY

The proposed methodology for the generation of metamaterial cells follows similar schemes to those developed by Esfarjani et al. (2022) and Chatterjee et al. (2021), which are organized into three stages: pre-processing, processing, and post-processing (Figure 1). In the pre-processing step, all essential parameters and data for the numerical model are defined. This involves setting the initial parameters for the finite element analysis (FEA), choosing the topological optimization method, and specifying the study parameters, such as the search space and the objective functions used in the optimization process. Further details are discussed in Section 5.

Figure 1
Flowchart of the proposed mechanical metamaterial cell generation strategy.

The processing stage is focused on the numerical process of metamaterial generation. The phase relies on the use of the Solid Isotropic Material with Penalization (SIMP) technique. However, the objective function inherent in the SIMP algorithm requires the integration of a numerical homogenization process leading to the evolution and creation of a novel cell structure. Additionally, numerical homogenization necessitates an FEA to obtain the mechanical analysis of the solid under multiple loading states.

Lastly, the post-processing stage involves selecting and treating the optimization results, which is divided into four steps. Firstly, the selection step focuses on choosing the best candidate based on superior mechanical behavior and a clearly shape. Subsequently, the filtering phase aims to eliminate intermediate densities generated during the topological optimization process. Surface smoothing is then applied (specifically in the 3D case) to improve the definition of the model, particularly in instances where the voxel size is significant. Finally, the idealization phase simplifies the previously treated result by establishing a manufacturable cell.

3 TOPOLOGICAL OPTIMIZATION STRATEGY

Topological optimization (TO) was developed by Bendsøe and Kikuchi (1988), Bendsøe (1989) and Suzuki and Kikuchi (1991) for simple solid structures with isotropy. However, TO has evolved by developing multiple topological optimization approaches, such as the evolutionary structural optimization (ESO) by Xie and Steven, (1993; 1997), bidirectional evolutionary structural optimization (BESO) by Querin et al. (1998) and Yang et al. (1999), the solid isotropic material with penalization (SIMP) method (Rozvany et al, 1992; Bendsøe and Kikuchi, 1988). Today, techniques such as hybrid methods (HM) (Shao, 2020) and variational topology optimization (VARTOP) (Yago et al., 2022) are available.

Recent studies by Yago et al. (2022) and Shao (2020), highlight the computational efficiency and stability of the SIMP method, which is recognized for its popularity and versatility. The ESO/BESO approaches can generate clearly defined topologies without intermediate densities. The VARTOP (Yago et al., 2022) and HM (Shao, 2020) methods are noted for their efficiency and robustness. Despite these advantages, some approaches may have inherent limitations due to the nature of their algorithms. For instance, when using discrete design spaces in ESO/BESO, stability issues may arise, leading to the creation of disconnected zones in the solid structure. Additionally, the SIMP method requires the use of filters to ensure a well-defined solution. In contrast, the level set (LS) method might exhibit an over-dependence on contour controls.

In the given scenario, the SIMP approach was chosen due to its conceptual simplicity, ease of implementation, and computational stability. The SIMP method seeks an optimal material distribution within a predefined design space, considering specific load cases, boundary conditions, manufacturing constraints, and performance requirements. Unlike shape and size optimization, topology optimization heavily relies on FEA, which enables the division of the design space (domain) into isotropic solid microstructures.

During the TO, the objective is to determine the optimal material state (filled or void) in each finite region or element. The material density distribution within the design domain, denoted as $\rho$, is typically represented discretely using binary values, where one (1) represents a filled state and zero (0) represents an empty state. However, $\rho$ can also be treated continuously, allowing for the consideration of intermediate densities. This continuous approach helps to avoid issues arising from the binary nature of the problem. A minimum value, ${\rho }_{min}$, is often assigned to represent an empty state to prevent the development of singularities in the modeling, while a value of one (1) represents a filled state.

Using a material with a relative density that can change continuously allows the material's mechanical properties to always vary. In this sense, the property that is penalized by $\rho$ would be the Young's modulus. The equation expressing the penalty is as follows:

$E\left({\rho }_e\right)={\rho }_e^P{E}_0$

where $P$ is the penalty factor and $E_0$ is Young's modulus. This penalty allows avoiding intermediate densities, forcing a full or empty distribution. In general, TOs follows the standard formulation as presented by Andreassen et al. (2014):

	$Min\to F\left(x\right),$
St.	$E\left({\rho }_e\right)={\rho }_e^P{E}_0,$	$e=\mathrm{1,2},\dots ,n$
	$\left[K\right]\left\{X^i\right\}=\left\{f^i\right\},$	$i=\mathrm{1,2},\dots ,n$
	$g_j\left(x\right)\ge 0,$	$j=\mathrm{1,2},\dots ,J$
	$\frac{1}{\mathrm{{\rm Y}}}\sum _{e=1}^n\left(v_e{\rho }_e\right)\ge V,$	$e=\mathrm{1,2},\dots ,n$
	${\rho }_{min}\le {\rho }_e\le 1,$	$e=\mathrm{1,2},\dots ,n$

where $\left[K\right]$ is the stiffness matrix, $\left\{{\chi }^i\right\}$ is the displacement vector, $\left\{f^i\right\}$ is the load case, $E_0$ is the Young modulus of the material and finally the design variable ${\rho }_e$ represents the state in which the material is found which is restricted between two values $\left({\rho }_{min}\le {\rho }_j\le 1\right)$ during the optimization process. The SIMP algorithm relies heavily on the application of the penalty factor in the calculation of the global stiffness matrix, which is modulated by the function:

$K\left(\rho \right)=\sum _{e=1}^n\left[{\rho }_{min}+\left(1-{\rho }_{min}\right){\rho }_e^P\right]K_e$

Additionally, the SIMP algorithm performs a sensitivity analysis to evaluate the impact of varying material densities on the objective function. Mathematically, the sensitivity analysis is expressed as the derivative of the objective function with respect to the densities, resulting in:

$\frac{dC}{d{\rho }_e}=-P{\left({\rho }_e\right)}^{P-1}C$

C represent the compliance, a measure of the overall flexibility or softness of a structure, usually quantified by the deformation energy. P is the penalty factor, which is generally three (3) for 2D cells and five (5) for 3D structures, and it enforces intermediate densities to converge toward extreme values. The TO code used in this research was developed in FORTRAN based on the work developed by Andreassen et al. (2011) and Liu and Tovar (2014). In this work, the optimality criterion (OC) was used in most of the generative processes, as ilustrated in section 6. In addition, in specific cases where isotropy conditions needed to be considered, the method of moving asymptotes (MMA) was implemented in accordance with the proposal of Svanberg (1987), as shown in section 6. Finally, to solve the system of equations, we chose to use the MA87 library (HSL, 2023) and the visualization of results was done using the software ParaView (Ayachit, 2005). The FORTRAN routines and Makefile for the compilation are available in the GitHub repository (Guevara-Corzo et al., 2023).

4 NUMERICAL HOMOGENIZATION PROCEDURE

Numerical homogenization processes have been used in the generation of cellular structures by Andreassen and Andreasen (2013) or Wu et al. (2021). According to the theory used to develop the homogenization process, as presented by Andreassen et al. (2014) and Yvonnet (2019), the homogenized stiffness tensor $E_{ijkl}^H$ of a periodic material is calculated by the equation:

$E_{ijkl}^H=\frac{1}{V}{\int }_0^v{E}_{pqrs}\left({\epsilon }_{pq}^{0\left(ij\right)}-{\epsilon }_{pq}^{\left(ij\right)}\right)\left({\epsilon }_{rs}^{0\left(kl\right)}-{\epsilon }_{rs}^{\left(kl\right)}\right)dV$

being V the cell volume, $E_{pqrs}$ the local tensor, ${\epsilon }_{pq}^{0\left(ij\right)}$ the preset macroscopic displacement field and ${\epsilon }_{pq}^{\left(ij\right)}$ the local variation of displacements. The latter is determined by the equation:

${\epsilon }_{pq}^{\left(ij\right)}=\mathrm{}{\epsilon }_{pq}\left({\chi }^{ij}\right)=\frac{1}{2}\left({\chi }_{pq}^{ij}-{\chi }_{qp}^{ij}\right)$

where ${\chi }^{ij}$ is the displacement field obtained by solving the system of equations with a prescribed macroscopic strain Writing the process in a way that facilitates its implementation. The homogenised stiffness tensor is also presented in simplified form in terms of the mutual energy of the element (Chatterjee et al., 2021). Using Voigt notation (e.g., for 2D we write $11\leftarrow 1$, $22\leftarrow 2$ and $12\leftarrow 3$) the stiffness tensor can be expressed as (Chatterjee et al., 2021):

$E_{ijkl}^H=\left[\begin{array}{ccc}{E}_{1111}^H& E_{1122}^H& E_{1112}^H\\ E_{2211}^H& E_{2222}^H& E_{2212}^H\\ E_{1211}^H& E_{1222}^H& E_{1212}^H\end{array}\right]=\left[\begin{array}{ccc}{Q}_{11}& Q_{12}& Q_{13}\\ Q_{21}& Q_{22}& Q_{23}\\ Q_{31}& Q_{32}& Q_{33}\end{array}\right]$

Where the sum of the mutual energies of the elements $Q_{ij}$ can be calculated with the equation:

$Q_{ij}=\frac{1}{V}\sum _{i=1}^n{U}_i^T{K}_e{U}_j$

Where $K_e$ is the stiffness matrix for each element and U is a matrix containing the displacement field for the different states of deformation (as shown in the example in Figure 2). The displacement fields are obtained by solving the equilibrium for the multiple deformation states.

Figure 2
Representation of Three Load Cases Used in Numerical Homogenization: a) Undeformed state; b) ε1 = [1,0,0]T ; c) ε2 = [0,1,0]T and d) ε3 = [0,0,1]T

For this, we obtain equivalent load cases calculated with the following equation:

$f^i=\sum _{i=1}^n{\int }_0^{V_e}{B}_e^T{C}_e{\epsilon }^{i}d{V}_e$

Where the matrix $B_e$ is the strain–displacement matrix, $C_e$ is the constitutive matrix for the element and ${\epsilon }^i$ is the strain state vector (see Figure 2).

5 MODELING: INITIAL CELL CONFIGURATION (BASE SOLID) AND OBJETIVE FUNCTIONS

The proposed base solids are a 2D square cell and a 3D cubic cell, both with 50 mm sides. The void within these structures is either circular or spherical in shape depending on the dimensionality of the structure. The void diameter (VD) was varied to analyze the TO behavior, performing analyses for structures with VD of 10 mm, 20 mm, and 30 mm (See Figure 3). Finally, the initial cell was assigned the mechanical properties of structural steel, with a Young's modulus of $E=2.1\mathrm{x}{10}^5$ MPa and a Poisson's ratio of $\upsilon =0.3$.

Figure 3
Initial Cell Configurations studied. (a) – (c): 2D cases with 10, 20, and 30 mm circular void diameters, respectively; (d) – (e): 3D cases initial cell cross-section with 10, 20, and 30 mm spherical void diameters, respectively.

In the TO process, a penalty factor $P$ of 3.0 and 5.0 was used for the 2D and 3D cases, respectively. A limit volume fraction (VF) of 30%, 40% and 50% was set for all analyses, and no stress constraints were applied. Additionally, all TO processes were limited to a maximum of one hundred iterations to ensure stability and a density vector change limit of 1E-3, which does not significantly affect the result. A filter radius of 1.5 times the full element size was employed, and the maximum finite element size was set to 1.0 mm (2D) and 2.0 mm (3D) for pixels and voxels, respectively (see Fig. 3). A classical finite element analysis formulation based on the displacement approximation was used. The meshing was performed using Ansys® and subsequently exported to the developed FORTRAN code (Guevara-Corzo et al., 2023).

Eleven objective functions were used to evaluate the influence of initial configuration and VF on the cell generation process in multiple cases. The choices of these functions were focused to generate cells with maximum stiffness in one or more directions (i.e., functions F1 and F2 in the 2D scenario and functions F1, F2, F3, F5 and F6 in the 3D scenario) and generating structures that maximize shear or bulk modulus, represented by functions F3 and F4 in 2D and functions F4 and F7 in the 3D scenarios, respectively. These approaches were taken from the work of Chatterjee et al. (2021) and Kollmann et al. (2020). The objective functions used in the TO are presented in Table 1.

Table 1
Objective Functions used for 2D and 3D scenarios. OF: objective function

2D OF	3D OF
$F_1=E_{1111}^H$	$F_1=E_{1111}^H$
$F_2=E_{1111}^H+E_{2222}^H$	$F_2=E_{1111}^H+E_{2222}^H$
$F_3=E_{1111}^H+E_{2222}^H+E_{1122}^H+E_{2211}^H$	$F_3=E_{1111}^H+E_{2222}^H+E_{3333}^H$
$F_4=E_{1212}^H$	$F_4=E_{1111}^H+E_{2222}^H+E_{3333}^H$ $\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{E}_{1122}^H+E_{2233}^H+E_{1133}^H$
	$E_{2211}^H+E_{3322}^H+E_{3311}^H$
	$F_5=E_{1212}^H$
	$F_6=E_{1212}^H+E_{2323}^H$

6 OPTIMIZATION AND IDEALIZATION RESULTS

The results of the processing and post-processing stages are presented below. For the 2D case, Figures 4, 5, 6 and 7 show the results of the topological optimization process, values of the objective function, idealization, and validation, respectively. Additionally, tables 2 and 3 present the values of the unnormalized and normalized elasticity tensor. Similarly, Figures 8, 9, 10 present the results of the optimization process for the 3D case. Figures 9 and 10 consider the isotropy constraint using the OC and MMA algorithms respectively. Finally, Figures 11, 12, 13 and 14 present the results obtained compared with mechanical metamaterial cells reported in the literature, along with the values of the objective function, the idealization, and the evolution of the topological optimization process.

Figure 4
Final Topology for 2D Cases. VD: void diameter; VF: volume fraction; OF: objective function.

Figure 5
Value of the OF (F1 to F4) and volume fractions (30% to 50%) for the 3D Cells: a) Void Diameter of 10 mm; b) Void Diameter of 20 mm, and c) Void Diameter of 30 mm.

Figure 6
Selection and Post-processing for 2D with 50% of VF, and 20 mm VD.

Figure 7
Evolution of the OF during the 2D TO process for the selected cases. a) F1. b) F2. b) F3 and d) F4. All with VF=50% and VD=20 mm

Figure 8
Final Topology for 3D Cases. VD: void diameter; VF: volume fraction; OF: objective function.

Figure 9
Comparison of F3 and F7 using isotropy restriction and OC algorithm.

Figure 10
Comparison of F3 and F7 using isotropy restriction and MMA algorithm.

Figure 11
Comparison of four 3D final topologies obtained in this work without isotropy constrains with those reported by Lu et al. (2022)

Figure 12
Value of the OFs (F1 to F7) and VF (30% to 50%) for the 3D Cells: a) VD of 10 mm; b) VD of 20 mm, and c) VD of 30 mm.

Figure 13
Selection and Post-Processing for the 3D Results.

Figure 14
Evolution of the OF during the 3D TO process. a) F1 (VF=50%, VD=20mm). b), c), d), e) and g) F2 to F5 and F7 (VF=50%, VD=10mm). f) F6 (VF=40%, VD=20mm).

According to the results for the 2D cells (see Figure 4), the F1 function influenced all final topologies to develop a double bar pattern in the X-direction. The bars´ robustness varies with the available volume fraction; they are slender for a low VF and robust for a high VF. For the F2 function, a crossbar topology was obtained in the X and Y direction, where it was one bar in each direction for the cells of 10 mm VD with 30% VF; 20 mm VD with 30% VF and 30 mm VD with 30% and 40% VF and double cross-bars for the remaining cases. For the F3 function, most cases develop a rhomboid pattern. However, as exceptions, the final topologies with a 30 mm VD and 30% VF generated a double diagonal bar pattern, while those with 10 mm and 20 mm VD and 50% VF resulted in a circular topology. Finally, the function F4 also developed a rhomboid pattern, highlighting that the topology with 10 mm VD and a VF of 50% set a double bar rhomboidal material distribution; the exception was the topology resulting from 30 mm VD with 30% VF which generated a double diagonal bar.

Note that some results obtained with the F2, F3, and F4 functions are consistent with works by Huang et al. (2011) or Chatterjee et al. (2021) who studied the generation of cellular structures with a different approach to the one used in this research. In addition, the TO employed in this work has generated topologies that have not been reported in the literature, such as those obtained with the F3 function, which has a VF of 50% and a VD of 10 mm and 20 mm.

On the other hand, Figure 5 indicates that the OF is directly proportional to VF. Conversely, the OF generally shows an inverse proportionality to VD. A 10% increase in VF from 30% to 40%, may result in an 8% to 45% rise in the OF. In opposition, altering VD can lead to a maximum reduction of 17% in the OF. Remarkably, an increase in OF was noted when VF was 40% for F2, 30% for F3 and both 30% and 50% for F4, with a maximum increase of 6% in OF.

For this first 2D instances, all the TO processes developed generated clearly defined topologies that have the potential to be applied as microgrid-like metamaterials. On the other hand, although all results were consistent, only those with a VF of 50% with a VD of 20 mm were selected to continue with the post-processing stage (see Figure 6). In addition, Figure 7 shows the behavior of the selected cases during the optimization process. It starts with the initial solid, followed by a stage of reduction of the VF until the assigned limit is reached. At this point, the TO process begins the generation of a clearly defined solid and finally the process stabilizes generating the final topology. In all cases a stable behavior is observed, reaching the final topology in less than 30 iterations. This pattern is also repeated in the 3D case, as can be seen in Figure 14.

Finally, the idealized cases (see Figure 6) were loaded into the Material Designer module of Ansys® to obtain the homogenized elasticity tensor. Unfortunately, this module does not allow the study of 2D cellular structures, so an idealized surface extrusion was performed, the numerical homogenization was done using the designer material and finally the tensor values related to the extrusion direction were omitted. The results are presented in the tables 2 and 3.

Table 2
Value of the ANSYS Homogenized Elastic Tensor (HET) for 2D Cases. HET is in Pa

HET	$F_1$	$F_2$	$F_3$	$F_4$
$E\left[\mathrm{1,1}\right]$	8.89E+10	4.60E+10	2.53E+10	1.75E+10
$E\left[\mathrm{2,1}\right]$	4.14E+09	4.75E+09	2.02E+10	1.68E+10
$E\left[\mathrm{2,2}\right]$	2.24E+10	4.60E+10	2.53E+10	1.75E+10
$E\left[\mathrm{3,3}\right]$	3.73E+08	9.76E+08	4.16E+09	1.48E+10

Table 3
Elastic Tensor Values Normalized (TVN) for 2D Cases. TVN was calculated with respect to the highest tensor value for each case.

HET	$F_1$	$F_2$	$F_3$	$F_4$
$E\left[\mathrm{1,1}\right]$	1.00	1.00	1.00	1.00
$E\left[\mathrm{2,1}\right]$	0.05	0.10	0.80	0.96
$E\left[\mathrm{2,2}\right]$	0.25	1.00	1.00	1.00
$E\left[\mathrm{3,3}\right]$	0.00	0.02	0.16	0.85

According to the results obtained by Ansys® (see Table 2) it is evident that the tensor values influenced by each OF are higher with respect to the other tensor components, except for the result with the OF F4, where E[3,3] was not the highest tensor value, but it was the highest E[3,3] of all results by far. In addition to the normalized values (see Table 3), it is possible to visualize that, the components with the highest ratio are indeed the homogenized components. The results for the case of 3D structures are presented below (Figure 8).

The topologies obtained with F1 developed a similar cylindrical pattern, with six being hollow and three solid. However, the cases with a VF of 40% and a VD of 30 mm, as well as a VF of 30% with VDs of 20 mm and 30 mm developed a separate structure. With F2, a dual plane pattern was obtained which, depending on the VF, could be either slender or robust, except for the topology obtained with 30% VF and 10mm VD which only developed one plane. F3 generated a topology characterized by a robust plane in one direction. However, the structures with a VF of 40% and 50% with a VD of 10 mm and 20 mm as well as the structure with a VF of 30% with a VD of 30 mm developed transverse elements, in some cases flat and in others following a cylindrical pattern. With F4, a similar topology was generated for all cases, characterized by the shape of three hollow cylinders aligned with the three axes. This pattern closely resembles that obtained by Huang et al. (2011) or the ´P shape´ reported by Lu et al. (2022). However, for structures with a 50% VF and a VD of 10 mm and 20 mm, and for a 40% VF with a VD of 20 mm, the topology was generated with a plan resembling reinforcement. The topologies generated with F5 (which are similar to the F2) are essentially the same as those generated by F2, in that they both develop a double plane structure. However, there is an exception: when generated with a 30% VF and a VD of 10 mm, only a single structure forms. Additionally, with VDs of 20 mm and 30 mm, a distinct cross-shaped pattern emerges within the plane. F6 produces well-defined topologies in all cases, which varied according to their initial configurations. The topologies generated with a 10 mm VD featured elements supported at the cell´s center. With a 50% VF and 20 mm VD, the resulting a topology was like the empty octet. The remaining topologies are similar to the octahedron reported by Lu et al. (2022) but without the horizontal elements. Finally, the F7 function varied according to their initial configuration. For the 20 mm and 30 mm VD, truncated octahedron or P shape cells were generated for the 20 mm and 30 mm VD initial cells. The structures with VD of 10 mm and VFs of 40% and 50% yielded designs with a central support and an octet, respectively, while the remaining cells generated a cross-shaped topology.

In situations involving F3 or F7 functions, additional constraints may be necessary. For instance, one could apply the isotropy constraint suggested by Andreassen et al. (2014), or employ of the Zener ratio as proposed by Wu et al. (2021). These constraints limit the algorithm´s search space to promote symmetric cell generation. The optimization model can enforce this by introducing a penalty factor to the objective function or by embedding an inequality constraint. By repeating the search indicated in figure 8, for the cases where F3 and F7 were employed, and using the isotropy constraint proposed by Andreassen et al. (2014), the results obtained were as follows (figure 9 and 10):

The additional constraints were successful in achieving their objective in all evaluated scenarios with F7, except for the case of 10 mm DV and 40% VF, where symmetric cell formation was observed. However, this strategy was inadequate when the F3 function was used; satisfactory results were only achieved with 10 mm VD. This finding emphasizes the significance of the initial conditions and the optimization algorithm used. It is important to note that when using the MMA algorithm, which is more robust but slightly more computationally expensive, satisfactory results were obtained in all cases. This resulted in the formation of symmetrical cells, although with less well-defined topologies at 30 mm VD.

The cells derived from the 3D TO process were well-defined. Additionally, previously unreported cases were observed, such as the reinforced P-shape obtained with F4 with both 10 mm and 20 mm VD, and the pattern produced with the F6 function with a VD of 20 mm and a VF of 50%. Additionally, it is evident that the proposed methodology is effective in the generation process, since structures close to the Hashin-Shtrikman limits are obtained, as reported by Berger et al. (2017). In Figure 11, we present four of the final cells obtained in this work, featuring OFs F3, F4, and F7, alongside the idealized topologies reported by Lu et al. (2022) and Huang et al. (2011).

According to Figure 12 and in the same way as in the 2D case, the VF is directly proportional to the OF value, obtaining even more marked increases ranging from 23% to 82% of OF value. However, for the 3D case, the OF Values does not present a clear tendency variation with VD where each case must be analyzed independently (Figure 12).

From the sixty-three generated topologies (see Figure 8), the following seven cases were considered for the post-processing stage: The structure with VF of 50% with 20 mm VD with OF F1, with 10 mm VD with OFs F2, F3, F4, F5 and F7 and finally with VF of 40% with 20 mm VD with OF F6. The results of the post-processing stage can be seen in the Figure 13.

After the post-processing stage, the numerical analysis was performed using the Ansys®. The results of the numerical analysis (homogenized elasticity tensor and normalized values) are presented in the tables 4 and 5. Similarly to the 2D case, the values of the tensor influenced by each objective function are larger with respect to the other tensor components. For the structures with OFs F5, F6, and F7, where the shear related components are sought to be maximized, these are not the largest in magnitude within the element, but they stand out when compared to the other structures. In addition, an examination of the ratios reveals the influence of the objective function. It is also noteworthy that optimizing these functions, which are related to the tensor components associated with shear, indirectly influences the tensor components linked to compression.

Table 4
Value of the ANSYS Homogenized Elastic Tensor (HET) for 3D Cases. HET is in Pa.

	$F_1$	$F_2$	$F_3$	$F_4$	$F_5$	$F_6$	$F_7$
$E\left[\mathrm{1,1}\right]$	1.07E+11	1.00E+11	6.80E+10	2.42E+10	1.00E+11	3.58E+09	7.42E+09
$E\left[\mathrm{2,1}\right]$	1.77E+10	2.79E+10	1.54E+10	1.17E+10	2.79E+10	3.38E+09	3.59E+09
$E\left[\mathrm{3,1}\right]$	1.77E+10	7.96E+09	1.54E+10	1.17E+10	7.96E+09	7.03E+07	3.59E+09
$E\left[\mathrm{2,2}\right]$	3.30E+10	8.89E+10	6.80E+10	2.38E+10	8.89E+10	7.16E+09	7.42E+09
$E\left[\mathrm{3,2}\right]$	2.61E+10	4.14E+09	1.54E+10	1.15E+10	4.14E+09	3.38E+09	3.59E+09
$E\left[\mathrm{3,3}\right]$	3.30E+10	2.24E+10	6.80E+10	2.38E+10	2.24E+10	3.58E+09	7.42E+09
$E\left[\mathrm{4,4}\right]$	2.25E+10	3.09E+10	1.41E+10	9.10E+09	3.09E+10	2.81E+09	3.29E+09
$E\left[\mathrm{5,5}\right]$	1.42E+10	3.74E+08	1.41E+10	9.04E+09	3.74E+08	2.81E+09	3.29E+09
$E\left[\mathrm{6,6}\right]$	2.25E+10	8.42E+09	1.41E+10	9.06E+09	8.42E+09	9.28E+07	3.29E+09

Table 5
Elastic Tensor Values Normalized (TVN) for 3D Cases. TVN was calculated with respect to the highest tensor value for each case.

	$F_1$	$F_2$	$F_3$	$F_4$	$F_5$	$F_6$	$F_7$
$E\left[\mathrm{1,1}\right]$	1.00	1.00	1.00	1.00	1.00	0.50	1.00
$E\left[\mathrm{2,1}\right]$	0.17	0.28	0.23	0.48	0.28	0.47	0.48
$E\left[\mathrm{3,1}\right]$	0.17	0.08	0.23	0.48	0.08	0.01	0.48
$E\left[\mathrm{2,2}\right]$	0.31	0.89	1.00	0.98	0.89	1.00	1.00
$E\left[\mathrm{3,2}\right]$	0.24	0.04	0.23	0.48	0.04	0.47	0.48
$E\left[\mathrm{3,3}\right]$	0.31	0.22	1.00	0.99	0.22	0.50	1.00
$E\left[\mathrm{4,4}\right]$	0.21	0.31	0.21	0.38	0.31	0.39	0.44
$E\left[\mathrm{5,5}\right]$	0.13	0.00	0.21	0.37	0.00	0.39	0.44
$E\left[\mathrm{6,6}\right]$	0.21	0.08	0.21	0.37	0.08	0.01	0.44

Additionally, an ANOVA analysis was performed using Minitab software (Mathews, 2004) with the data obtained at the end of the optimization process (presented in Figures 5 and 12) to determine the influence of the factors studied. The analysis was considered separately for 2D and 3D, taking FO, VD and VF as factors. An analysis of terms of up to second order (interaction of two (2) factors) was considered in the response with a significance level of 0.05. The results are presented in Figure 15.

Figure 15
Pareto Diagrams of standardized effects. a) 2D cases; b) 3D cases.

The factors with the greatest influence on the value of the OF (in both the 2D and 3D cases) are the functions used, followed by the VF and the OF-VF interaction. The behavior of the remaining terms varies depending on the case. In the 2D case, the term with the greatest influence is the VD, ending with OF-VD and VD-VF interactions, which are below the threshold cutoff line of 2.18 (indicating that these interactions are not significant on the final response). To conclude, for the 3D case, the terms that follow are the OF-VD, VD-VF interactions, and, finally, the VD that have little influence, as they are under the threshold limit line of 2.06.

The results presented here are validated by comparison with previous studies by Huang et al. (2011), Yu et al. (2018), Vangelatos et al. (2021), Wang et al. (2022). Specifically, in 3D for the solution of the OF F7, the octet-tissue unit cell (Deshpande et al., 2001; Vangelatos, et al., 2021; Wang et al., 2022) or stretch-dominated octet-truss unit cell (Yu et al., 2018), known for its high strength-to-weight ratio (Yu et al., 2018; Vangelatos, et al., 2021) was obtained. For the solution of the 3D OF F4, a cell like the Schwars P-shell (Wang et al., 2022) was found; this design was also reported by Huang et al. (2011) as optimal for maximum bulk material design. However, the collected literature does not report cases such as the cell obtained from the OFs F3 (see Figure 6), F4 and F6 (see Figure 13).

7 CONCLUSIONS

Based on the results obtained, the following conclusions can be inferred.

1. The influence of the initial cell configuration was evident. Both the VF and the VD significantly affected the values of the OFs, which are directly proportional to the VF. However, in certain cases like F1 for the 2D case and F2 and F5 for the 3D cases, the influence of VD and VF was minimal, generating topologies with insignificant variation. When working with more complex functions, such as those maximizing the bulk modulus and the shear modulus, the initial configuration´s influence became more pronounced, even resulting in a topology not previously reported in the literature reviewed in this research.

2. The use of an initial cell with a central void to generate metamaterial cells is an effective and robust approach, providing a solid foundation for developing well-defined cell topologies and ensuring a reliable performance across both low and high VF. Notably, the best results were achieved at high VF (50%) for both 2D and 3D models, as well as at low and medium VDs (10mm and 20mm).

3. The TO process developed runs automatically with minimal user intervention. It systematically, rapidly, and stably generates topologies that are comparable to those reported in the literature (refer to Figure 8). In the postprocessing stage, new cells were obtained with the OFs F3 (see Figure 6), F4 and F6 (see Figure 13); these cells are not documented in the collected literature. However, during the TO process, cells with exceptionally low densities are produced (see Figure 8 with OFs F2 and F5 for any VD and VF), resulting in idealized cells that may seem disconnected. Despite this appearance, they remain structurally unaltered during the filtering process, and their mechanical performance is unaffected.

4. For objective functions F3, F4, and F7, additional constraints must be implemented to obtain results consistent with the problem; however, these functions, being the most complex, render the OC method inadequate for cell generation. Therefore, the use of more robust algorithms such as MMA, GCMMA or similar is highly recommended. Although these algorithms may be more computationally expensive, they provide accurate and reliable results.

Although the research was limited to a linear elastic approach, it is important for future research to extend the optimization/generation process to more complex models, such as geometric nonlinearity with or without specific physical non-linearities. Additionally, it is important to extend this process to metamaterials designed to manipulate dynamic properties. Pre-established forms are often used in this context, yet the methodology for their generation warrants further exploration. For instance, future work could focus on phenomena like the manipulation of natural frequencies to intentionally create specific bands gaps.

Finally, evaluating the behavior of the generated cells in experimental tests is important. Exploring potential applications is also crucial, particularly the standout potential of the 3D results. These results reveal cells clearly constituted by planes, suggesting they could be explored and applied in origami or kirigami-type structures, which have been successfully used in impact absorption.

ACKNOWLEDGEMENTS

All authors would like to thank Universidad Industrial de Santander-UIS, INME-UIS, and MSU-UNIFEI research groups for their technical support during the development of this research. Additionally, the first author is grateful for the financial support provided by Universidad Industrial de Santander.

REFERENCES

• Ai, L., & Gao, X. (2019). Topology optimization of 2-D mechanical metamaterials using a parametric level set method combined with a meshfree algorithm. Composite Structures, 229.
• Andreassen, E., & Andreasen, C. (2013). How to determine composite material properties using numerical homogenization. Computational Materials Science, 83, 488-495.
• Andreassen, E., Clausen, A., Schevenels, M., Lazarov, B., & Sigmund, O. (2011). Efficient topology optimization in MATLAB using 88 lines of code. Structural and Multidisciplinary Optimization, 43(1), 1-16.
• Andreassen, E., Lazarov, B., & Sigmund, O. (2014). Design of manufacturable 3D extremal elastic microstructure. Mechanics of Materials, 69(1), 1-10.
• Ansys®. (2023). Academic Research Mechanical, Release 23.1.
• Ayachit, U. (2005). The ParaView Guide: A Parallel Visualization Application. 4
• Bendsøe, M.P. (1989). Optimal shape design as a material distribution problem. Structural optimization, 1, 193-202.
• Bendsøe, M.P., & Sigmund, O. (2002). Topology optimization: theory, methods, and applications. Springer.
• Bendsøe, M.P., & Kikuchi, N. (1988). Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering, 71(2), 197-224.
• Berger, J., Wadley, H., & McMeeking, R. (2017). Mechanical metamaterials at the theoretical limit of isotropic elastic stiffness. Nature, 543(7646), 533-537.
• Bertoldi, K., Vitelli, V., Christensen, J., & van Hecke, M. (2017). Flexible mechanical metamaterials. Nature Reviews Materials, 2(11), 1-11.
• Chatterjee, T., Chakraborty, S., Goswami, S., Adhikari, S., & Friswell, M. (2021). Robust topological designs for extreme metamaterial micro-structures. Scientific Reports, 11(1).
• Chen, W., & Huang, X. (2019). Topological design of 3D chiral metamaterials based on couple-stress homogenization. Journal of the Mechanics and Physics of Solids, 131, 372-386.
• Cummer, S., Christensen, J., & Alù, A. (2016). Controlling sound with acoustic metamaterials. Nature Reviews Materials, 1(3), 1-13.
• Deshpande, V., Fleck, N., & Ashby, M. (2001). Effective properties of the octet-truss lattice material. Journal of the Mechanics and Physics of Solids, 49(8), 1747-1769.
• Esfarjani, S., Dadashi, A., & Azadi, M. (2022). Topology optimization of additive-manufactured metamaterial structures: A review focused on multi-material types. Forces in Mechanics, 7
• Fabbrocino, F., Amendola, A., Benzoni, G., & Fraternali, F. (2015). Seismic application of pentamode lattice. International Journal of Earthquake Engineering, 62-70.
• Fraternali, F., Amendola, A., & Benzoni, G. (2018). Innovative seismic isolation devices based on lattice materials: A review. Ingegneria Sismica, 35(4), 93-113.
• Gibiansky, L., & Sigmund, O. (2000). Multiphase composites with extremal bulk modulus. Journal of the Mechanics and Physics of Solids, 48(3), 461-498.
• Guevara-Corzo, J., Begambre-Carrillo, O., Garcia-Sanchez, J., & Quintero-Ramirez, C. (2023). Mechanical Metamaterials Generation Code in Fortran based on SIMP and Numerical Homogenization Retrieved from https://github.com/jeffreyguevaracorzo/NH-TOPOPT-FORTRAN
  » https://github.com/jeffreyguevaracorzo/NH-TOPOPT-FORTRAN
• Guevara-Corzo, J., Begambre-Carrillo, O., Garcia-Sanchez, J., & Sanchez-Acevedo, H. (2022). Passive seismic protection systems with mechanical metamaterials: A current review. Structural Engineering and Mechanics
• HSL. (2023). A collection of Fortran codes for large scale scientific computation. Available.
• Huang, X., Radman, A., & Xie, Y. (2011). Topological design of microstructures of cellular materials for maximum bulk or shear modulus. Computational Materials Science, 50(6), 1861-1870.
• Ishida, S., Suzuki, K., & Shimosaka, H. (2017). Design and Experimental Analysis of Origami Inspired Vibration Isolator with Quasi-Zero-Stiffness Characteristic. Journal of Vibration and Acoustics, 139(5), 1-5.
• Jenett, B., Cameron, C., Tourlomousis, F., Rubio, A., Ochalek, M., & Gershenfeld, N. (2020). Discretely assembled mechanical metamaterials. Science advances, 6(47).
• Kadic, M., Bückmann, T., Stenger, N., Thiel, M., & Wegener, M. (2012). On the practicability of pentamode mechanical metamaterials. Applied Physics Letters, 100(19).
• Kadic, M., Milton, G., van Hecke, M., & Wegener, M. (2019). 3D Metamaterials. Nature Reviews Physics, 1(3), 198-210.
• Kollmann, H., Abueidda, D., Koric, S., Guleryuz, E., & Sobh, N. (2020). Deep learning for topology optimization of 2D metamaterials. Materials & Design, 196.
• Li, Z., Luo, Z., Zhang, L., & Wang, C. (2021). Topological design of pentamode lattice metamaterials using a ground structure method. Materials and Design, 202.
• Liu, K., & Tovar, A. (2014). An efficient 3D topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 50(6), 1174-1196.
• Lu, C., Hsieh, M., Huang, Z., Zhang, C., Lin, Y., Shen, Q., Zhang, L. (2022). Architectural design and additive manufacturing of mechanical metamaterials: a review. Engineering, 17, 44-63.
• Martin, A., Kadic, M., Schittny, R., Bückmann, T., & Wegener, M. (2013). Phonon band structures of three-dimensional pentamode metamaterials. Physical Review B - Condensed Matter and Materials Physics, 86(15), 2-6.
• Mathews, P. G. (2004). Design of Experiments with MINITAB. Quality press.
• Milton, G. (2016). Analytic materials. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472(2195).
• Milton, G., & Cherkaev, A. (1995). Which Elasticity Tensors are Realizable? Engineering Materials and Technology, 117(4), 483-493.
• Miniaci, M., Kherraz, N., Cröenne, C., Mazzotti, M., Morvaridi, M., Gliozzi, A., Pugno, N. (2021). Hierarchical large-scale elastic metamaterials for passive seismic wave mitigation. EPJ Applied Metamaterials, 14, 8.
• Querin, O., Steven, G., & Xie, Y. (1998). Evolutionary structural optimization (ESO) using a bi-directional algorithm. Engineering Computations, 15(8), 1031-1048.
• Rozvany, G., Zhou, M., & Birker, T. (1992). Generalized shape optimization without homogenization. Structural optimization, 4, 250-252.
• Schittny, R., Bückmann, T., Kadic, M., & Wegener, M. (2013). Elastic measurements on macroscopic three-dimensional pentamode metamaterials. Applied Physics Letters, 103(23).
• Shao, G. (2020). Comparison of BESO and SIMP to do structural topology optimization in discrete digital design, and then combine them into a hybrid method. Architectural Intelligence: Selected Papers from the 1st International Conference on Computational Design and Robotic Fabrication (CDRF 2019).
• Sigmund, O. (1994). Materials with prescribed constitutive parameters: An inverse homogenization problem. International Journal of Solids and Structures, 31(18), 2313-2329.
• Sigmund, O. (1995). Tailoring materials with prescribed elastic properties. Mechanics of Materials, 20(4), 351-368.
• Suzuki, K., & Kikuchi, N. (1991). Shape and topology optimization by a homogenization method. Computer Methods in Applied Mechanics and Engineering, 93, 291-318.
• Svanberg, K. (1987). The method of moving asymptotes—a new method for structural optimization. International journal for numerical methods in engineering, 24(2), 359-373.
• Ungureanu, B., Achaoui, Y., Enoch, S., Brûlé, S., & Guenneau, S. (2016). Auxetic-like metamaterials as novel earthquake protections. EPJ Applied Metamaterials, 2(17).
• Vangelatos, Z., Sheikh, H., Marcus, P., Grigoropoulos, C., Lopez, V., & Flamourakis, G. (2021). Strength through defects: A novel Bayesian approach for the optimization of architected materials. Science advances, 7.
• Wang, Y., Zhang, X., Li, Z., Gao, H., & Li, X. (2022). Achieving the theoretical limit of strength in shell-based carbon nanolattices. Proceedings of the National Academy of Sciences, 119(34), 1-11.
• Wu, J., Luo, Z., Li, H., & Zhang, N. (2017). Level-set topology optimization for mechanical metamaterials under hybrid uncertainties. Computer Methods in Applied Mechanics and Engineering, 319, 414-441.
• Wu, S., Luo, Z., Li, Z., Liu, S., & Zhang, L. (2021). Topological design of pentamode metamaterials with additive manufacturing. Computer Methods in Applied Mechanics and Engineering, 377.
• Xie, Y., & Steven, G. (1993). A simple evolutionary procedure for structural optimization. Computers and Structures, 49(5), 885-896.
• Xie, Y., & Steven, G. (1997). Basic Evolutionary Structural Optimization. Springer London.
• Yago, D., Cante, J., Lloberas-Valls, O., & Oliver, J. (2022). Topology optimization methods for 3D structural problems: a comparative study. Archives of Computational Methods in Engineering, 29(3), 1525-1567.
• Yang, X., Xie, Y., Steven, G., & Querin, O. (1999). Bidirectional evolutionary method for stiffness optimization. AIAA Journa, 37(11), 1483-1488.
• Yu, X., Zhou, J., Liang, H., Jiang, Z., & Wu, L. (2018). Mechanical metamaterials associated with stiffness, rigidity and compressibility: A brief review. Progress in Materials Science, 94, 114-173.
• Yvonnet, J. (2019). Computational Homogenization of Heterogeneous Materials with Finite Elements. Springer.
• Zhang, L., Bai, Z., & Chen, Y. (2022). Dual-functional hierarchical mechanical metamaterial for vibration insulation and energy absorption. Engineering Structures, 271.
• Zheng, X., Lee, H., Weisgraber, T. H., Shusteff, M., DeOtte, J., Duoss, E. B., Jackson, J. A. (2014). Ultralight, ultrastiff mechanical metamaterials. Science, 344(6190), 1373-1377.
• Zheng, X., Smith, W., Jackson, J., Moran, B., Cui, H., Chen, D., Weisgraber, T. (2016). Metamaterials, Multiscale metallic. Nature materials, 15(10), 1100-1106.

Edited by

Editor: Marco L. Bittencourt

Publication Dates

• Publication in this collection
  21 June 2024
• Date of issue
  2024

History

• Received
  03 Nov 2023
• Reviewed
  16 Feb 2024
• Accepted
  23 Apr 2024
• Published
  03 May 2024