Multiobjective Design Optimization of Reentrant Auxetic Model Using Lichtenberg Algorithm Based on Metamodel

. The optimization of five different responses of an auxetic model was considered: mass; critical buckling load under compression effort; natural frequency; Poisson ’s ratio; and failure load. The Response Surface Methodology was applied, and a new meta-heuristic of optimization called the Multi-Objective Lichtenberg Algorithm was used to find the optimized configuration of the model. It was possible to increase the failure load by 26,75% in compression performance optimization. Furthermore, in the optimization of modal performance, it was possible to increase the natural frequency by 37.43%. Finally, all 5 responses analyzed simultaneously were optimized. In this case, it was possible to increase the critical buckling load by 42.55%, the failure load by 28.70% and reduce the mass and Poisson ’s ratio by 15.97% and 11%, respectively. This paper shows something unprecedented in the literature to date when evaluating in a multi-objective optimization problem, the compression and modal performance of an auxetic reentrant model.


INTRODUCTION
Auxetic structures are those that contract transversely when subjected to compression forces in the axial direction or expand when subjected to traction (Burns,1987). These structures are called auxéticas and several studies have been published on this topic. According to Francisco et al. (2021), the relationship between transverse and longitudinal deformations is called Poisson's ratio and the auxetic structures have a negative Poisson's ratio (NPR).
The NPR makes auxetic models have several improved properties, such as: shear, indentation and fracture resistance, energy absorption, etc. All these properties have been explored by several authors and several applications have been suggested, for example, bulletproof vests, car bumpers, impact absorbers, among others. Two important works showing some applications was done by Duncan et al. (2018) and Francisco et al. (2021). The first carried out a review with several auxetic for sports applications and the second make a review of auxetics applied in energy absorption problems.
In the same line, Dong et al. (2019)  shows that the crushing stress can be underestimated when the nominal stress-strain curve method is applied.
Several studies have been done recently. Wang et al. (2019) developed a novel reentrant auxetic honeycomb with enhanced in-plane impact resistance, Lekesiz et al. (2017) carried out a mechanical characterization of auxetic stainless steel thin sheets reentrant structure, Bronder (2021) focused on a design study for multifunctional 3D reentrant auxetics, Jiang et al. (2020) studied the crashworthiness of novel concentric auxetic reentrant honeycomb with negative Poisson's ratio biologically inspired by coconut palm and Baertsch et al. (2021) carried out a finite element modeling and optimization of 3D printed auxetic reentrant structures with stiffness gradient under low velocity impact.
Auxetic structures are extensively studied submitted to compression, as for example,  studied the behavior of auxetic structures under compression and strong impact. Furthermore, reentrant structures are widely applied due to their ease of construction and because they are widely known. There are several types of reentrant structures as shown in the review by Kelkar et al. (2020) on cellular auxetic structure for mechanical metamaterials.
Due the hexagonal reentrants have excellent performance, Guo et al. (2020)  This manuscript is organized as follows: section 2 presents the background of auxetic hexagonal model, response surface methodology and Lichtenberg algorithm. Section 3 shows the methodological procedure of this work. Section 4 presents the main results and discussion.
Finally, Section 5 draws the conclusions.

Auxetic Reentrant Structures
In 80's, Gibson et al. (1982) was one of the pioneers to study the hexagonal reentrant auxetic structures. The mechanical properties of these types of structures were verified via analytical and experimental methods. The author developed models that describe the desired properties in terms of bending, elastic buckling and plastic collapse. The Figure 1a shows the structure proposed by the author.
The design parameters of the structure proposed by Gibson et al. (1982) are shown in Figure 1a and are the length of the horizontal bar (b), length of the oblique bar (l) and angle between bars (θ). This angle can be positive, forming a conventional hexagonal structure, or negative, forming a structure with auxetic behavior. When the structure is subjected to compression, the angle tends to increase, reducing the diameter of the structure as well.
An auxetic structure is formed by several equal unit cells, thus, evaluating only one cell it is possible to describe the behavior of the complete structure. In this way, Gibson et al. (1982) developed analytical models that could explain the deformations of a single cell when subjected to stress. The authors considered that each part of the unit cell is a bar and described Equations 1 to 4 that describe, respectively, modulus of elasticity in the direction of force, modulus of elasticity in the direction transverse to force, Poisson's ratio, and the shear modulus of the structure.
The modulus of elasticity and Poisson's ratio in direction 1 and 2 can be related by In addition, Gibson et al. (1982) made some considerations to develop the analytical model: i) the bars deform only along the axis ii) there is no angular movement of the bars and iii) the structure is under stretching efforts. Thus, considering that the cell is subjected to tensile stresses in the vertical direction, the force exerted due to the applied stress σ2 is given by Equation 6 and the P component of w that acts along the diagonal bar is given by Equation 7 It is also possible to find all the parameters of the structure when it is subjected to stresses in the horizontal direction. This procedure will not be presented in this work, as the objective is not to develop the analytical model of the problem. The Equations shown in this work show that the Poisson's ratio depends on the length of the horizontal bar (b) of the cell, the length of the diagonal bar (l) and the angle between them (θ). Thus, this manuscript aims to optimize the properties of an auxetic structure with hexagonal reentrant cells varying the 3 parameters mentioned.

Response Surface Methodology
Design of experiments is a statistical method for optimally performing experimental tests. This technique has several arrangements, among them the response surface methodology where x = (x1, x2, …, xk) is the vector of control variables, β is the Equation coefficients and ε is the model error. The β coefficients are found by analytical methods (three-level factorial design, box-Behnken design, central composite design, and Doehlert design). The authors of this work adopted the central composite design (CCD). Considering k as the number of controllable factors, we have that the CCD has 2 k factorial points, 2k axial points and one central points.
The present work consists of 2 3 factorial points, 6 axial points and one central points with 5 replicates. An analysis will also be done through the analysis of variance (ANOVA) to verify which factors are important for the process and which can be eliminated. The model fit is checked by the adjusted R² parameter, the closer to 100%, the better.

Multi-objective Lichtenberg Algorithm Optimization (MOLA)
The Multi-objective Lichtenberg Algorithm has trajectory and population characteristics and was created by . The single-objective version of the algorithm was proposed by the same author   All the values must be adopted take into account the complexity of the problem.

Numerical Modelling using Finite Element Method
The structure was numerically modeled in a software based on the finite element method that encompasses everything from the creation of the structure to the creation of the mesh and the results. The chosen beam element is based on Timoshenko beam theory that is an extension of Euler-Bernoulli beam theory that allows the shear-deformation effects and has six degrees of freedom at each node. Figure 3 shows the geometry, node locations, and coordinate system of the element that will be used in this work. PLA was adopted to carry out this work.
The 3D printing manufacturing is a process that changes the properties of the original material due to the build-orientation selected for fabrication and anisotropic nature of the fused deposition modeling (FDM) process. In this way, to find the properties of the material used, 5 specimens were built (Figure 4a) using the same process and the same parameters that will be used to build the auxetic structure, the property data found will be used to feed the numerical model according to ASTM D695 (Figure 4b). Resultant materials properties are listed in Table   1 and the values are used in the Finite Element Analysis.  The structure has a rectangular cross section of 110.00 × 49.52 mm² and a depth of 30 mm. Furthermore, the construction of the reentrant cell depends on 3 parameters: angle (θ), oblique length (l) and base length (b), as shown in Figure 1. Numerical modeling was done according to experimental tests to be a faithful representation of reality. The numerical analyzes were according to ASTM C365 which is a standard test method for flat-wise compression properties of sandwich cores and the crosshead speed is equal to 5mm/min. A study of buckling in the structure was made using the eigenvalue associated with the structure stiffness matrix, i.e., the maximum buckling load is given by the eigenvalue multiplied by the applied load. Equation 15 and 16 shows the relationship between eigenvalue (Λ) and critical buckling load (λ).
where [K] is the global stiffness matrix of the structure and [KG] is the global geometric stiffness matrix and the value is proportional to the initial force Fi. Thus, the eigenvector θ correspond the buckling mode and the eigenvalue Λ is associated to the buckling load (λ).

Response Surface Design
The objective of this work is to find the optimal set of parameters for the reentrant cell so that the structure has the best possible performance. As already mentioned, the cell parameter set is formed by the length of the oblique bar, the length of the horizontal bar and the angle between them (the parameters are represented in Figure 1). and 30mm, respectively. The author adopted 7 cells in the horizontal and 3 in the vertical direction. Table 2 shows the parameters adopted by Guo et al. (2020). Based on the values in Table 2, the 3 parameters that will serve as input to the RSM will be varied. Thus, the length of the oblique bar of the unit cell will vary from 6 to 10mm, the angle will vary between 50 and 70°, the length of the base will vary between 14 and 18. The central composite design (CCD) using RSM was developed considering the bounds for the three variables as related before as shown in Table 3.

Model Validation
Numerical analysis was performed to represent reality. The authors took all possible precautions so that the numerical model had results that were faithful to the experimental tests.
Thus, the ASTM C365 standard was adopted as a rule for the analysis, this is a standard test method for flat-wise compression properties of sandwich cores. For example, the crosshead speed in compression analysis was equal to 5mm/min. The results obtained by the numerical and experimental analyzes are shown in Table 4 (θ = 60°, b = 16mm and l = 8mm). Figure 6 shows the deformation of the part during the experimental test. (e) t = 80s (f) t = 100s (g) t = 120s (h) t = 140s It is possible to notice that the errors obtained between the experimental tests and the numerical model are low, as shown in Table 4. This difference obtained is normal in this type of analysis, as there are some experimental errors and some errors in the numerical calculations that influence the results. Furthermore, the deformation of the structure was analyzed, and the final results are shown in Figure 7. Note the similarity that exists between Figures 6g and 6h with Figures 7a and 7b, showing that the numerical model is in accordance with the experimental analysis.

RSM Analysis
The objective functions that will be used in the optimization were obtained from the response surface methodology, five different models that represent each of the responses of the hexagonal reentrant auxetic structure will be analyzed and shown in Table 5. The five responses are: mass (M), Poisson's ratio (ν), Critical buckling load (λ), natural frequency (ωn) and failure load (Qu). All these models were analyzed to verify the fit to the data and the results are shown in Table 6.  It is noticed that the adjustments found were high, all above 95%. Thus, the models shown in Table 5 will be used to perform the multi objective optimization of the model. As shown in Figure 8a, only the angle is not statistically significant for mass. This is because the amount of material deposited does not change with varying the angle. Furthermore, in Figure 8b it is possible to notice that the angle has no influence on the natural frequency either, this is because ωn is strongly influenced by the mass and, as the angle has no significance for the mass, neither for the natural frequency. In Figure 8c, it is possible to see that all factors are important for the critical buckling load. Furthermore, it can be noted that the length of the oblique bar is the factor that most influences this response, and the angle is the factor that has the least influence.
Along the same lines, Figure 8d shows that the length of the horizontal bar of the unit cell has no influence on the failure load and the length of the oblique bar has a very large influence on this result. Thus, it is possible to notice that the modification of length l is what has the most influence for the failure load and for critical buckling load. However, the length of the horizontal bar has no significance for the failure load. Finally, Figure 8e shows that the angle has a great influence on the Poisson's ratio and the length of the horizontal bar has no significance for this response. This happens because the angle plays a fundamental role in the deformation of the structure, i.e., the variation of the angle generates a greater or lesser deformation of the model.
It is important to emphasize that the angle has a great importance for the Poisson's ratio, that is, the modification of ϴ causes great interference in the results of ν. However, this same significance of the angle is not seen for the other responses. It is not significant for mass and natural frequency and has little significance (relative to other factors) for failure load and critical buckling load. On the other hand, the length of the oblique bar is the factor with the greatest influence on all answers, except for the Poisson's ratio, which has the angle at this position.
Furthermore, the length of the horizontal bar b has discrete influence on some responses. It is not the most influential, but it does influence mass, natural frequency, and critical buckling load. And it has no influence on failure load and Poisson's ratio. This happens because the load is being applied perpendicularly to this bar, meaning that it does not play an important role in the failure resistance of the structure or in the deformation of the model. Figure   9 shows how each of these factors impact responses.   Figure 9 shows how the factors influence the response, for example, in Figure 9a it is possible to notice that increasing the length of the horizontal and vertical bar also increases the structure's mass. This is because the construction of a structure with a longer length requires a greater amount of material. On the other hand, in Figure 9b, it is possible to notice that the increase in the length of the oblique bar decreases the natural frequency, this is because ωn is inversely proportional to the square root of the mass, that is, the increase in mass causes a decrease in the natural frequency, as if you can notice in that case. It is possible to notice that the modification of some factors impacts on several answers at the same time, which is why a multi objective analysis of the problems is very important. Furthermore, the alteration of several simultaneous responses when a factor is modified may indicate a correlation between them. Thus, an analysis was carried out to assess whether there is (or not) correlation between the responses and the results are shown in Table   7. Figure 10 shows the dendrogram to assess the similarity between the responses. From Figure 10, the greatest degree of similarity is between critical buckling load, failure load and natural frequency. These three responses are farther from Poisson's ratio and mass. The quantification of the correlation between the responses is shown in Table 4. The lower value corresponds to the p-value of the correlation test, thus, if the value is less than 0.05, it means that there is a correlation between the response. The higher number is the Pearson coefficient, the closer to unity (negative or positive) the stronger the correlation. From the results shown in Table 5, the mass is negatively correlated with the critical buckling load, failure load and natural frequency. This means that the increase in mass causes these responses to decrease. This is mainly because of the length of the oblique bar. The increase in the length l causes the mass to increase, however, this same increase causes a decrease in the three responses mentioned above. The Poisson's ratio is only correlated with the failure load, i.e., the greater the failure load, the greater the Poisson's ratio. This happens because increasing the angle and the length of the oblique bar, both Qu and v decrease.
Furthermore, the critical load is positively correlated with the failure load and natural frequency and negatively with the mass. That is, increasing the critical buckling load means increasing the failure load, increasing the natural frequency and decreasing the mass. This result is very important, as it shows that by modifying one response (critical buckling load), all the others will improve together, since the focus is on a structure with greater resistance to compression (greater λ and greater Qu) and far from resonance (greater ). The full optimization analysis will be done in the next topic.

Multi Objective Design Optimization
The multi objective Lichtenberg Algorithm was used to find the set of optimal parameters of the model studied in this work. The five answers were divided into three cases: compression performance, modal performance and multi objective optimization. Several nondominated solutions (optimal points) have been found and will be demonstrated via the Pareto's fronts. The criterion chosen to determine the best point among the best was TOPSIS.

Case I: Compression performance
This section is focused on evaluating the compression performance of the auxetic structure. The objective is to maximize the critical buckling load, maximize failure load and minimize the Poisson's ratio. It is important to point out that there are many great points in a multi-objective optimization, and it is impossible to improve all the answers simultaneously. In this case, increasing the buckling load and failure load means increasing the value of the Poisson's ratio. All results obtained by MOLA are shown at the Pareto's front (Figure 11). All points shown in Figure 11 are considered optimal and the best choice depends on the application of the structure. It is possible that a structure with a lower Poisson's ratio or a structure with a higher failure load is needed and there are different configurations for each of these options. In the present work, the TOPSIS criterion was adopted to decide the optimal configuration of the optimal ones and the result is shown in Table 6 and Table 7.  It is noticed that the MOLA generated results very consistent with the finite element analysis, i.e., the relative error between the responses was small. It is important to emphasize that the results found are quite expressive, showing that the structure is resistant to buckling and withstands a high load before failure, in addition to having a low Poisson's ratio. Figure 12 shows the of the initial and the optimized structure. Finally, it is important to note that a small change in structure causes a large change in the results, which is why optimization analysis is so important these days. Furthermore, Table   11 shows the results obtained in the structure initially proposed and the structure optimized in this section. Realize that it was possible to improve all the answers in the structure. The focus of this section is to optimize the modal performance of the structure. The objective is to minimize the mass and Poisson's ratio, while the natural frequency is maximized.
All the great results found are shown on the Pareto's front shown in Figure 13. All the responses found are optimal configurations of the structure, however, each configuration has different results. There is no gain without loss, so improving one requirement means making the other worse. Using the TOPSIS decision criterion, we have the configuration shown in Table 12 and the results are shown in Table 13.  The results obtained by the MOLA algorithm were close to those obtained by the finite element method. This shows that the model generated was very faithful to the data and the results are reliable. The model on initial and optimized structure are shown in Figure 14.
(a) (b) Figure 14 -The model of (a) initial and (b) optimized structure for modal performance.
The optimal configuration is 2mm less both in the length of the oblique bar and in the length of the horizontal bar and angle between bars a little larger. Thus, the model was chosen according to the TOPSIS criterion and the comparison between the results of the initial and the optimized structures is shown in Table 14. Note that all answers have been improved in relation to the initial configuration.

Case III: Multi objective optimization
The purpose of this section is to optimize the five responses simultaneously, that is, one must minimize mass, Poisson's ratio and maximize failure load, critical buckling load and natural frequency. It is not possible to construct the Pareto's front in 5 dimensions, however there are also several optimal points and, using the TOPSIS criterion, the configuration shown in Table 15 was reached and the results are shown in Table 16.  It can be noted that the results obtained by MOLA and by finite element analysis are very close, showing reliability in the results. The comparison of the results of the optimal configuration and the initial configuration is shown in Table 17. It is possible to notice that all the responses were improved with the optimization process in relation to the initially proposed structure. The smallest improvement was the Poisson's ratio, which was 11% lower, and the biggest improvement was the critical buckling load, which was more than 42% higher. Finally, this section shows how an optimization process can generate a major improvement in the product without impairing its performance. the Figure 15 show the initial and the optimized structure.
(a) (b) Figure 15 -The model of (a) initial and (b) optimized structure.

CONCLUSIONS
The optimization of an auxetic reentrant structure was performed using the response surface methodology in conjunction with the finite element method and the Multi objective Lichtenberg Algorithm. In addition, a comparison was made with an auxetic reentrant model proposed by another author to assess the improvement in structure.
The models generated by RSM proved to be fully accurate and reliable. The results obtained by the optimization of the Equations were close to the results obtained by the finite element analysis. The R 2 (adj) parameters of all models were above 95%, and this was an important factor for the continuity of the work.
The model created for the finite element analysis was validated with experimental data and generated acceptable results. The RSM was able to adjust well to the data which means that the model generated reliable data for each proposed configuration. The FEM proved to be an especially important technique in this project, as it was possible to carry out several experiments at a low cost.
The MOLA was able to find optimal points that proved to be reliable, since both the results generated by the algorithm and the ones generated by the finite element model were close. MOLA is an algorithm recently developed that can be exploited in several problems of mechanical optimization.
Finally, the results obtained in this work show the best possible structures considering more than one objective simultaneously. With the optimization analysis, it was possible to increase the critical buckling load by 42.82% and the failure load in compression performance by 26.75%. Furthermore, in the optimization of modal performance, it was possible to increase the natural frequency by 37.43% and decrease the mass by 15.97%. Finally, all 5 responses analyzed simultaneously were optimized. In this case, it was possible to increase the critical buckling load by 42.55%, the failure load by 28.70% and reduce the mass and Poisson's ratio by 15.97% and 11%, respectively.