Bio-Inspired Morphological Evolution of Metastructures with New Operation Modalities

the

properties with given boundary condition, loading, and constraints. [22] Although topology optimization has been gaining more attention in metamaterial design, [21,[23][24][25] it is heavily dependent on parameter initialization. A number of researchers have developed machine learning and deep learning approaches hybridized with topology and metaheuristic optimization methods for metamaterial design and discovery. [26,27] While these approaches provide acceptable performance, they are based on the data alone to determine the optimal metamaterial structure. Furthermore, they can explore a design space within the limits that they have been calibrated for. The creation of novel forms and designs is a missing part component in these efforts because they cannot be deployed to directly evolve the material structure.
Natural evolution has been a major source of inspiration for scientists for decades. The resulting field, evolutionary computation, has been successfully applied to solve a wide range of engineering problems. [28,29] In the arena of metamaterials, evolutionary computation methods have been applied to optimize only specific properties for single designs in each evolution, [30] such as negative index, [31,32] acoustic gradient index, [33] and large permeability. [34] Besides, recent efforts have been merely focused on using evolutionary computation to predict the mechanical properties of these metastructures by simply evolving the predicting variables. [35] There is a critical shortage in research needed to harness the power of natural evolution for the automated exploration of novel forms of metamaterial systems.
Here, we introduce the concept of "evolving metamaterial (EM)" to directly evolve thousands of metamaterial structures via an evolutionary process akin to biological evolution. Figure 1a,b shows the vision toward the artificial evolution of metamaterial systems, where an evolutionary approach follows the natural evolution process to automatically creature metamaterial prototypes. The proposed EM approach is different from the previous efforts in that it is based on "direct morphological evolution of metamaterial microstructure." It can operate using merely a piece of matter. Therefore, unlike supervised machine learningbased approaches, it does not require a database for training purposes, and it is not limited to exploring a design space only within the training limits. Embedding the principles of evolution in metamaterial design and discovery potentially enables exploring hitherto unknown structures, new modalities of operation, or potential solutions that humans would find inconceivable. In the EM approach, evolution takes place by randomly creating an initial population of metamaterial entities that pass on their genetic material to their offspring through variation, reproduction, and selection. The metamaterial configurations with desired responses emerge during this evolutionary process.

Results and Discussion
Without loss of generality, we focus on the design and discovery of mechanical metamaterial structures. The proposed concept can be arguably applied to other metamaterials with different geometrical designs. The EM method merely requires a piece of matter to launch the design process. For mechanical metamaterials, this piece can be a representative unit cell (RUC). The equivalent constitutive behaviors of an RUC can be readily evaluated using the homogenization theory and periodical boundary conditions (PBCs), [36] as described in the Experimental Section. Here, an RUC is a cellular structure that serves as the "parent metamaterial" in the evolutionary process. Figure 1c shows the EM process for generating RUCs. The process involves population initialization and evolution phases. A large number of initial (parent) RUCs generated during the population initialization phase are fed into the evolutionary process to create new "child" or "offspring" RUCs. EM resembles biological evolution as each RUC in a cellular structure contains the same chromosomes making up the genetic code. Thus, each offspring metamaterial structure inherits the genetic information of its parent cells. This base parent cell should be initially created at random. Stage I in Figure 1c presents an example of how only a single "base parent RUC 1" with a simple shape, e.g., a circle, can be used to launch the EM process to maximize the bulk modulus. The base parent RUC 1 can be copied to create a "copy of base parent RUC 1" during Stage I. The genetic operators such as crossover or mutation can then be applied to the base parent RUC 1 (Parent 1) and its copy (Parent 2) to produce child metamaterials. However, it is more efficient to initially create a population of randomly generated parent metamaterial structures to achieve high diversity, i.e., N base parent RUCs with various basic shapes, as shown in Figure 1c (Stage I). In Stage II, topology optimization is performed on each design based on the parameters assigned in Stage I to obtain the initial population. Then initial population is passed to Stage III to initialize the evolutionary process. Evolutionary operations are performed in this stage to generate new designs. New designs and old designs are passed to Stage IV for evaluation and ranking based on their fitness value, e.g., maximum bulk modulus. The best designs in terms of their fitness are sent back to Stage III for the next cycle of evolution. The iteration is terminated when the desired improvement is achieved, or the defined generation number is reached. Figure S1-S3, Supporting Information, present the flowcharts, pseudocodes, and further description of the EM process for the artificial evolution of metamaterial structures.
The EM method runs a tournament, i.e., best metamaterial structure selection process, based on the defined design criteria. The selection mechanism of EM is mainly inspired by natural selection. In nature, individuals with higher fitness have a higher probability to survive. Therefore, their genes are more likely to be passed to the next generation. To mimic this selection, the wellknown roulette wheel mechanism is also used in EM. Higher probabilities are assigned to individuals with higher fitness (design criteria), and the chance of selection is proportional to their fitness values. Since a roulette wheel is a stochastic operator, individuals with poor fitness still have small probabilities to be selected. Instead of discarding poor solutions, the roulette wheel selection mechanism keeps the diversity of populations. In this step, for example, four metamaterial structures are selected from the initial population of RUCs randomly. They are compared and based on design criteria, two structures are picked as the winners (parents) and two as the losers. The parent metamaterials are transformed using the genetic operators: crossover or mutation.
The genetic operators used in the EM method are similar to those used in traditional genetic algorithms. However, the major difference between these two methods is that the EM algorithm evolves metastructures while the genetic algorithm can only generate coded binary strings of numbers to represent the solution.
A genetic algorithm is generally used for parameter optimization to find the best values for a given set of model parameters. During the mutation procedure, the EM algorithm occasionally selects a portion from an RUC at random and mutates it or changes it. In each cycle of the EM procedure, the loser metamaterial structures are replaced in the tournament with the transformed winner metamaterial structures. The winners of the tournament remain without change. Each of these cycles is called a generation. EM evaluates the evolved metamaterial structures, selects them for reproduction, generates new structures, and finally creates a new generation in all iterations. EM implements many generations until a termination criterion is met. The best metamaterial structure according to the defined design criteria that appeared in any generation defines the output of the EM framework. In the proposed framework, a population of the initial designs can be randomly created from a pool of ten basic Figure 1. Artificial morphological evolution of metamaterial structures. a) Artificially evolved metamaterial design processing steps, including initial creation of a random population of unit cells based on design requirements, applying genetic operators to transform the designs, evolution of the structure prototypes, and design evaluation. b) Transforming the parent metamaterial structures using the crossover and mutation operators. c) Flowchart of the EM process for generating RUCs with maximum bulk modulus. V f , penal_range, r min _range, and c are user-defined topology optimization parameters denoting volume fraction, penalty factor, radius of filter, and optimization objective, respectively. E is the constitutive parameter.
www.advancedsciencenews.com www.advintellsyst.com shapes (circle, oval, triangle, rhombus, square, rectangle, trapezoid, pentagon, hexagon, and octagon) for the RUCs. Since evolution is a relatively slow form of adaptation, the process can be optimized by assigning good starting RUCs in terms of the design criteria for evolution. To this aim, we use topology optimization with randomly initialed parameters inside the EM process. Here, topology optimization serves as a tool to provide starting seed RUCs for higher diversity and then to optimize the generated RUCs. The detailed topology optimization process to achieve this goal is described in Materials and Methods. We consider three design criteria to crate proof-of-concept 2D mechanical metamaterial prototypes: maximum bulk modulus, maximum shear modulus, and minimum Poisson's ratio. The design objectives are not limited to the criteria considered in this study and depend on the targeted application. Figure 2a-c shows the three galleries of randomly generated RUC designs for various objectives using topology optimization. The parameter definitions and settings for the topology optimization are presented in (Table S1 and S2, Supporting Information). A constant value of 0.5 was considered for the volume fraction as a tradeoff between the computational cost and volume constraint. This parameter can range from 0 to 1 depending on the design preference. Each gallery shows 16 different RUCs. These initial RUCs are passed to the evolutionary process, where they are transformed using genetic operators. The evolved new cells are topology optimized again to become the new generation of RUCs. The RUCs in the new generation are evaluated based on their fitness (maximum bulk modulus, maximum shear modulus, minimum Poisson's ratio). The evolutionary process is iterated until requirements are achieved. Since evolutionary computation is a population-based method, it does not only evolve RUC designs that satisfy the design criteria, but also generates thousands of other potential RUC designs that each can be a hitherto unknown structure with unique features. In order to test the proposed EM framework, extensive simulations were conducted. Various combinations of the evolutionary parameters such as initial population, number of generations, crossover rate, and mutation rate were considered. The parameter settings for the evolutionary process are presented in (Table S3- Figure 2d-f ).
The conducted evolutionary process resulted in more than eighty thousand metastructures for the defined design objectives and parameter settings. To validate the evolved designs, 12 optimal RUC designs with maximum bulk modulus (Samples B1 to B4), maximum shear modulus (Samples S1 to S4), and minimum Poisson's ratio (Samples P1 to P4) were selected from the entire population of evolved structures. Figure 3a shows a 3D model of a typical evolved design sample. The mechanical properties of the selected designs were calculated using the homogenization theory presented in the Experimental Section. Figure 3b shows the selected designs and their theoretical values. The dimension of 3D printed samples is 80 mm Â 160 mm Â 10 mm. Compression tests were conducted to determine in-plane bulk modulus and Poisson's ratio for samples B1 to B4 and P1 to P4. Figure 4a-c shows the initial stages and deformed stages of all samples during the tests. Black points were painted on the samples for video gauge tracking. Poisson's ratio was calculated as follows [37,38] where ε transverse and ε longitudinal are the transverse and longitudinal strains, respectively. In-plane bulk modulus was evaluated by compression stress and compression strain defined in Ref. [39] using the following equation where p is pressure, ΔV is the volume change under pressure, V 0 is the initial volume. In-plane shear modulus was calculated using the following equation [40] where b, h, P x , ε transverse , and ε longitudinal are the width, thickness, loading along the x-axis, transverse strain, and longitudinal strain of the sample, respectively. Figure 5a shows the testing results to determine the bulk modulus, shear modulus, and Poisson's ratio for the chosen 2D mechanical metamaterials. More details about the procedure to calculate these parameters are provided in the Experimental Section. Table 1 provides a summary of theoretical and experimental bulk modulus, shear modulus, and Poisson's ratio values. The error rate ranges for the bulk modulus, shear modulus, and Poisson's ratio values are 0.1%-4%, 5%-6%, and 1%-5%, respectively. Overall, there is an acceptable agreement between the experimental and theoretical values. The difference between theoretical and experimental values can potentially stem from the precision of the 3D printing and errors associated with experimental testing equipment.
To investigate the capabilities of the proposed EM concept for discovering 3D mechanical metamaterial, the entire framework has been extended into 3 dimensions. The evolutionary process was examined to explore a suite of proof-of-concept 3D mechanical metamaterial prototypes with maximized bulk modulus. The general working mechanism and design of the 3D prototype are similar to those described for the 2D models. The equivalent constitutive behaviors of an RUC in 3D can also be evaluated using the homogenization theory and PBCs, as described in Materials and Methods. To test the proposed 3D framework, extensive simulations were conducted. Various combinations of the evolutionary parameters such as initial population, number of To show a wide range of potential designs explored during the EM framework, four 3D designs with minimum and maximum bulk moduli have been selected. Figure 6a shows the 3D models   of the discovered architectures with minimum (Samples 1 and 2) and maximum (Samples 3 and 4) bulk modulus among more than three thousand structures that evolved during the evolutionary process. The mechanical properties of the selected designs were calculated using the homogenization theory presented in Materials and Methods. The discovered designs were then fabricated. The dimension of the 3D printed samples is 60 mm Â 60 mm Â 60 mm. Compression tests were conducted to determine the bulk modulus. To ensure the constraint for the compression tests, a specific testing fixture was designed. Figure 6b shows the initial stages and deformed stages of all samples during the tests. The bulk modulus was evaluated using the measured compression stresses and strains, as defined in Ref. [39] and using Equation (2). Figure 6c shows the pressure versus the volume changing ratio and the curve fitting results for 4 samples. Comparing Figure 6a,b, the error rates between the theoretical and experimental bulk modulus values for Samples 1-4 are 1.80%, 3.79%, 1.97%, and 0.84%, respectively.  www.advancedsciencenews.com www.advintellsyst.com

Summary and Outlook
In summary, we have demonstrated that EM is a viable approach to directly evolve metamaterial designs through a natural evolutionary process. Here, evolution designs a 2D metamaterial artifact which is assessed during evolution with respect to maximum bulk modulus, maximum shear modulus, and minimum Poisson's ratio. In addition, the feasibility of the proposed EM concept in 3D has been illustrated by evolving 3D configurations with various bulk moduli. We showed that the experimental properties of the evolved systems are in acceptable agreement with the theoretical values. A wide range of user-defined parameters allows for more flexibility in adapting various design objectives. One key advantage of the EM concept is that a large population of potential designs is generated during the evolution process. While it is possible to choose specific solutions as the outcome of the EM process, many of the other discovered solutions can possess new and unknown modalities of operation for real-life applications. For example, Sample P2 was not the best metamaterial design based on the defined negative Poisson's ratio criteria, but it is still a reliable option for applications such as energy absorbing at multiple scales (e.g., aircraft, knee pads, highway joints) (see Figure 7a). With the growing interest in metamaterials for various engineering applications, the artificial evolution of metamaterials could open up new avenues toward the more efficient and creative design of these systems. The entire search process can be enhanced by introducing multiple objectives for the evolution. The next frontier stemming from the EM concept is the direct evolution of structural representation at a multi-scale rather than using RUCs as a seeding point. Since genetic operators will act on the entire structural form, the EM approach can result in the creation and exploration of unpredictable and non-repeatable novel forms. This goal can also be achieved by passing more complex genes to the evolutionary process. As an example, minimal surface designs (e.g., Scherk's minimal surfaces) can be used as the initial population in a more open-looped artificial selection process to create a completely new form of hitherto unknown artificial metamaterial systems www.advancedsciencenews.com www.advintellsyst.com ( Figure 7b). However, natural-evolution-inspired methods are computationally intensive and parameter sensitive. Parallel computing and using optimization algorithms to find the optimal parameter settings during the evolutionary process can help tackle these limitations.

Experimental Section
Fabrication and Testing of the Evolved Metamaterial Systems: Fused deposition modeling (FDM) was deployed to 3D print the evolved prototypes. The samples were fabricated using a Raise3D Pro2 3D Printer. Three steps were followed for the fabrication of samples including 3D modeling using SolidWorks, 3D printing of the designs using a Raise3D Pro2 3D Printer, and removal of the 3D printed object supports and extra printed parts. Thermoplastic polyurethane (TPU) (E = 26 MPa, u = 0.48) was used for printing the final designs. To ensure the in-plane constraint for the compression tests, a specific testing fixture was designed, and 3D printed. To evaluate bulk performance for each sample, loading and compressed displacement results were obtained from the testing machine. For shear modulus calculation, ε transverse and ε longitudinal related to one unit cell were obtained. P x and 2bhðε transverse À ε longitudinal Þ were used to obtain the shear modulus, as described in Ref. [40]. To obtain the Poisson's ratio, ε transverse and ε longitudinal were used. ASTM D3518/D3518M-13 standard was followed to carry out to calculate the homogenized in-plane shear modulus via an off-axis 45°t ensile test. Instron 8874 Biaxial Servo-hydraulic Fatigue Testing System was used to conduct all the tests. The loading direction was perpendicular to the printing orientation of the bulk modulus and Poisson's ratio specimens. The filaments were oriented parallel to the load direction for the shear modulus samples. Videos were recorded by a Sony ILCE-7M3 camera model. An open resource video gauge Ncorr [37] was used to track the strain fields along the loading and transverse direction. The proposed EM framework was programmed in MATLAB. High-performance computing was performed for the numerical simulations at the University of Pittsburgh Center for Research Computing (CRC). The objective functions (For 2D: Equation (9)-(11), for 3D: Equation (13)) were directly programmed and computed in MATLAB. Mechanical Characterizations of the Evolved Metamaterial Systems: Within the range of linear elasticity, the periodically patterned metamaterial can be represented by a single cell Y, called RUC. The equivalent constitutive behaviors were evaluated using the homogenization theory and PBCs. [36] The homogenized stiffness tensor E H ijkl of RUC is defined as where |Y| denotes the area (or volume in 3D) of the periodic base cell domain Y. E pqrs is the constitutive parameters, ε oðklÞ pq corresponds to three (2D) or six (3D) linear independent unit strain fields, ε ÃðklÞ pq is the periodical solution to the variation problem.
In finite element analysis, the RUC is discretized into N finite elements and Equation (4) can be approximated by T k e u AðklÞ e (6) where u AðklÞ e is element displacement solutions corresponding to the unit strain fields ε oðklÞ pq , k e is the element stiffness matrix. Modified solid isotropic material penalization (SIMP) [41] method is used to solve the optimization problem. The RUC is discretized into N elements and the same size density design variables are defined as ρ. The element Young's modulus is defined as E e ðρ e Þ ¼ E min þ ρ p e ðE 0 À E min Þ where E 0 is Young's modulus of solid material and E min is Young's modulus of the Ersatz material which is an approximation for void material. [42] p is a penalization factor used to drive the density distribution. The optimization problem is presented as follows min ρ ∶cðE H ijkl ðρÞÞ s:t:∶KU AðklÞ ¼ F kl , k, l ¼ 1, : : : , d v e ρ e jYj ≤ ϑ ∶0 ≤ ρ e ≤ 1, e ¼ 1, : : : , N where K is the global stiffness matrix,U AðklÞ and F kl are the global displacement vector and external force vector for ðklÞ test case. d is the spatial dimension, v e is the element volume, and ϑ is the upper bound of the volume fraction. The homogenized stiffness tensor cðE H ijkl ðρÞÞ is the objective function.
In 2D cases: The maximization of bulk modulus corresponds to the objective function The maximization of the shear modulus corresponds to the objective function c ¼ ÀE 1212 (10) Optimizing designs with negative Poisson's ratio is a challenging task, as it requires imposing other constraints such as bulk modulus and/or isotropy. [43][44][45] In this study, the objective function of minimizing Poisson's ratio is written as Here, β ∈ ð0, 1Þ is a fixed parameter defined by the user, and 0.8 is chosen in this study. The exponential l is the design iteration number.
The sensitivity of the objective function is computed using the adjoint method [46] ∂E H Where k 0 is the element stiffness matrix for an element with unit Young's modulus.