Auxetic behavior obtained through the large deformations of variants of the rectangular grid

Abstract The deformation behavior of intersecting ligaments forming variants of the square and rectangular grids under mechanical compression was investigated. It was shown that such systems are able to exhibit a negative incremental Poisson’s ratio at relatively large axial compressive strains. Numerical simulations and experimental studies indicated that the extent of auxeticity depends on the relative offset of successive ligaments, the relative lengths of the ligaments as well as on their thickness. It was also shown that there are two distinct modes of deformation, one resembling that of the reentrant hexagonal honeycomb and the other that of the meta-tetrachiral system.

Apart from the work on elastic instabilities, the deformation mechanisms used to attain an auxetic behavior generally involved small deformations of the structures. As a matter of fact, different studies have indicated that the Poisson's ratio is highly dependent on the strain [63][64][65][66][67][68], which can be of concern for certain applications such as impact loading [22][23][24][25]. This has prompted efforts to reduce its variation with deformation with the aim of producing structures that are able to keep their Poisson's ratio constant [63,69,70]. The perception that in some applications the variation of the Poisson's ratio with strain is an undesirable feature might be due to the fact that in most cases the Poisson's ratio tends to increase with increasing strain. Yet, if this were not the case, and the Poisson's ratio can be made to decrease, then there can actually be scope for an improved response of the material. An example is given by impact loading where the indentation resistance would be expected to increase with the induced strain as the Poisson's ratio becomes more negative.
This work looks precisely in this direction by aiming at harnessing large strains to induce auxetic behavior in systems where it is considered not to be present, at least for small deformations. For the purpose, the chosen geometries were those of the square and rectangular grids (Figure 1). Wang and McDowell [71] have shown that the effective Poisson's ratio of the configuration illustrated in Figure 1(a) is positive and depends on the Poisson's ratio of the material and the dimensions of the ligament; while Ghaedizadeh et al. [56] and Ho et al. [57] obtained auxetic behavior from it by making the joints more rigid. On the other hand, the systems shown in Figure 1(b,d) emerge in the context of auxetic systems as a special case of the hexagonal honeycomb that occurs in the transition between the conventional and reentrant honeycombs where the cell wall angle is zero [28,30,72,73]. For this case, from the theoretical standpoint, the Poisson's ratio for the bulk system when loaded axially in the horizontal direction is zero [28,30,72]. (The result does not necessarily hold for finite systems where edge effects would tend to alter the Poisson's ratio [74].) Yet it is clear that for large compressive strains in the horizontal direction, the system will acquire a reentrant configuration, meaning that a negative Poisson's ratio is expected to emerge.
Based on these considerations, this work investigates other square and rectangular honeycombs configurations that can exhibit auxetic behavior with the aim of showing that a negative Poisson's ratio can result from large deformations. The systems considered were first studied through nonlinear FEA. Subsequently a number of prototypes were 3D printed and mechanically tested in order to confirm the numerical results.

The structure
For this study, a more general form of the rectangular grid will be considered whereby the vertical and horizontal ligaments will have different thickness as shown in Figure 2. In principle it can be expected that for large compressions a negative Poisson's ratio would be attained if successive horizontal ligaments are simply offset with respect to one another. This follows from the fact that a compression along the horizontal axis would be expected to cause the vertical ligaments to bend due to the three-point force system (represented as F 1 , F 2 and F 3 in Figure 2) acting on them, thereby inducing a lateral contraction. In view of these considerations, this generalized rectangular grid was investigated further using finite element analyses (FEA) and experimental methods as detailed in the next sections.

Finite element analysis
In order to study the deformation of the structure, nonlinear simulations were carried out using the finite element software ANSYS V R Academic Research Mechanical APDL Release 13, as this is a widely used technique. The minimum repeating unit cell needed to represent the system, shown in Figure 2, was adopted and periodic boundary conditions were used in order to save on the computational time. For the same reason, the length of the vertical ligament l 2 was kept constant at 25 mm while that of the horizontal ligament l 1 was varied between 10 mm and 50 mm in increments of 2.5 mm. This allowed the aspect ratio of the rectangle, given by to vary between 0.4 and 2 while minimizing the number of structures being simulated. For the sake of simplicity, the thicknesses of the horizontal and vertical ligaments (t 1 and t 2 respectively) were taken to be equal and henceforth referred to as t. The thickness was then assigned values of 1 and 2.5 mm so that was set to 0.04 and 0.1. Regarding the offset d, this was given an initial value of 1 mm which was then increased by 1 mm each time until it reached the midpoint of the vertical ligament. It should be noted that in order for the horizontal ligament to reach the midpoint of the vertical one, the last value of d was not necessarily an integer multiple of 1 mm since the maximum value of the offset is (l 2 þt)/2. Hence, the associated non-dimensional parameter, took values between 0.04 and (1 þ g)/2. For the bulk of the simulations, the Young's modulus, E S , was taken to be 1.6 GPa while the Poisson's ratio S was set to 0.45. These values were based on those obtained from various dog bones made from the material meant to be used for the experimental tests. The values needed for the software to model the nonlinear behavior of the material, i.e.  the stress-strain data for the structural material, were also obtained in the same way. Subsequently, in order to further investigate if and how the results obtained are dependent on the material properties, additional simulations were carried out by changing the Poisson's ratio to 0 and À0.4 while keeping the Young's modulus constant at 1.6 GPa and using the same nonlinear data. Furthermore, additional materials were considered, namely aluminum AA 6061-T651, data for which was obtained from Aakash et al. [75,76] and steel 4130 data for which was obtained from Majzoobi and Alavinia [77] as quoted in Seifi and Salimi-Majd [78]. For these simulations t was set to 1 mm (g ¼ 0.04), l 2 was kept constant at 25 mm while l 1 was assigned values 10, 25, 35 and 50 mm (a ¼ 0.4, 1.0, 1.4, and 2.0 respectively). Three different offsets were considered with d taking values 13, 10 and 7 mm (v ¼ 0.52, 0.40, and 0.28 respectively). Further details are provided in the Supplementary Material. PLANE183, a higher order 2D element with 8 or 6 nodes each of which has two degrees of freedom [79] was the chosen element type. This element type has a quadratic displacement behavior that is suitable for modeling irregular meshes, so that it is frequently used in this type of analysis [80][81][82]. Furthermore, plane stress was assumed. Periodic boundary conditions were applied, these being analogous to the ones adopted by Cauchi et al. [73] and discussed at length by Mizzi et al. [83]. This choice of boundary conditions allowed the modeling of the structure as a 2D infinite system. Subsequently mesh convergence tests were conducted in order to determine the element size that yielded results that were overlaid on those obtained with a finer mesh. In consequence, an element size equal to 0.2 mm was adopted for the simulations presented here.
A compressive force was then applied along the Ox 1 direction and a nonlinear solution obtained. Analysis of the results indicated that the transverse strain e 2 varied nonlinearly with the axial strain e 1 . For this reason, a specific application software was written using the programming language MATLAB in order to determine the incremental Poisson's ratio ( 12 ¼ À de 2 de 1 ) around specific intervals of the values of e 1 , namely 0, 0.01 and 0.02. This was done by carrying out a linear regression on the points within a specified region of values of values of e 1 with the possible inclusion of the nearest neighboring points.

Production of the samples and mechanical testing
In order to investigate further the behavior of the structures, for the main analysis, two physical prototypes were manufactured. The dimensions were chosen on the bases of the results of the FEA simulations (see Section 3) and were meant to be representative of the two modes of deformation that emerged from the analysis. They were produced using an Ultimaker 3 FFF (fused filament fabrication) 3D printer with Ultimaker PLA (polylactic acid) filament being chosen for the structural material. During the manufacturing process, the extruder temperature was set to 205 C while that of the bed temperature to 60 C. The infill density was set at 100% and the layer height was chosen to be 0.2 mm. The printing adhesive MagigooV R 3D was applied on the printing bed so as to prevent the structure from warping.
Based on the FEA analysis and the specifications of the 3D printer, the following dimensions were chosen for the two structures: Structure 1 (Figure 3(a)) had an offset For both structures the length of the ligaments was l 1 ¼ l 2 ¼ 25 mm (a ¼ 1) and their thickness t ¼ 1.2 mm (g ¼ 0.048). The out of plane thickness was set to 15 mm. Based on the printing bed dimensions, three-unit cells in the Ox 1 direction and six in the Ox 2 directions were used. Furthermore, for the purpose of distributing the applied force over the whole length of the structure, a 10 mm bar was included along the faces perpendicular to the loading direction. In addition to these two structures, another two prototypes were manufactured to complement the main analysis. These had mostly the same dimensions outlined above with the only exception being the offset which was set equal to 10.2 mm (v ¼ 0.408; Structure 3) and 7.2 mm (v ¼ 0.288; Structure 4). Additional information is provided in the Supplementary Material.
The printed structures were mechanically tested under compression with a Testometric universal loading machine (M350-20CT) and a 1000 N load cell (Serial Number: 31931). A set of two white markers was applied in the loading direction while two sets (each containing two) of white markers were used in the transverse direction. These defined a central unit cell which was used for the calculation of the Poisson's ratio. The structures were compressed till they buckled. Their deformation was recorded with a calibrated Messphysik-Videoextensiometer camera. This together with the pattern recognition feature found within the Videoextensiometer software allowed the determination of the changes in the lengths, yielding one set of readings for the axial direction and two sets of readings for the transverse direction.
To further analyze the results, nonlinear regression was used to fit a polynomial to the data. The order of the polynomial was chosen so that the function provided an adequate fit to the data being analyzed, while at the same time care was taken to avoid overfitting. Differentiating the polynomial then provided an estimate of the incremental Passion's ratio in the region being studied. In the case of the FEA results, the smooth variation between the points allowed the use of interpolation. For this purpose, the interpolation function found in Mathematica 12.1, with polynomials of order three was adopted.

Results and discussion
The typical variations of the transverse strain with the longitudinal strain when the system is compressed are illustrated in Figures 4 and 5 for different l 1 (a), t (g) and d (v). From these it can be noted that for small axial strains (i.e. for values less than around 1%) the change in the transverse strain is very small so that e 2 stays close to zero. It can be further noted that in this region the value of e 2 appears to be slightly positive, meaning that the incremental Poisson's ratio is also positive (see also Figure 6).
As the structure is compressed further, the value of the transverse strain begins to decrease so that it eventually becomes negative. Thus, in this region, the incremental Poisson's ratio becomes negative since the gradient of the graph is positive, meaning that the Poisson's ratio is going from a positive near zero value to a negative value. It should be noted that this occurs as soon as e 2 begins to decrease even if its value is still greater than zero. Figures 4 and 5 also indicate that when l 1 is smaller than l 2 (a ¼ 0.4) the variation of e 2 with e 1 becomes less prominent as d (v) attains larger values. This means that 12 becomes less negative with increasing d. As the magnitude of l 1 increases as compared to l 2 the spread in the variation of e 2 with e 1 for different offsets decreases, so that the values of the incremental Passion's ratio approach each other. Eventually,   when l 1 is relatively much larger than l 2 (a ¼ 2.0), the initial trend appears to be reversed, meaning that 12 becomes more negative as d increases. It can also be noted by comparing the two figures that when the thickness of the ligaments increases, the values of e 2 become more positive.
In order to quantify these observations, an application program was written to determine the incremental Passion's ratio around e 1 ¼ 0 (by taking the region of value of e 1 ranging between 0 and 5 Â 10 À4 ), e 1 ¼ 0.01 (by taking the region of value of e 1 ranging between 0.01 and 0.0105) and e 1 ¼ 0.02 (by taking the region of value of e 1 ranging between 0.02 and 0.0205) using the methodology discussed in Section 2.2. The results corresponding to the structures considered in Figures 4 and 5 are shown in Figure 6 (the corresponding graphs using nondimensional variables can be found in the Supplementary Material). As can be noted, the two structures have rather similar behaviors. For small axial strains 12 is positive. The Poisson's ratio also exhibits a general trend of increasing slightly with decreasing offset and increasing length l 1 (a). Deviations from this behavior can be observed at smaller values of d (v) and large values of l 1 whereby as d increases, the incremental Poisson's ratio first exhibits a moderate increase before it starts decreasing. Figures 6(c,d) indicate that at around 1% axial strain the incremental Poisson's ratio is predominantly negative. The only exception occurs for small l 1 and larger thicknesses. In fact, 12 appears to attain a maximum at small l 1 . Then, as l 1 starts to increase the incremental Poisson's ratio decreases and continues to do so for larger d while for smaller offsets it attains a minimum before it starts increasing again. It can further be noted that the minimum in the value of 12 shifts toward smaller values of d and/or larger values of l 1 as the thickness of the ligaments is increased. This is also accompanied by a decrease in the range of values the incremental Poisson's ratio attains.
Considering the deformation at a larger axial strain of 2% it can be noted that the behavior of the incremental Poisson's ratio is similar to that observed when e 1 ¼ 0.01. The differences lie mainly in the position of the minimum which is shifted toward smaller values of l 1 and/or larger values of d. In addition, more negative values are attained.
To better understand the origin of the variations in the incremental Poisson's ratio, the simulated deformed structures were studied in detail. As illustrated in Figure 7, the structure was observed to deform through a combination of two distinct modes. When the offset is small, the bending of the ligaments is analogous to that exhibited by the meta-tetrachiral [84] with both the vertical and the horizontal ligaments being bent (Figure 7(a)). On the other hand, when the offset ligament is at the center of the other ligament, the deformation is analogous to that of the reentrant honeycomb, with only the ligament in the Ox 2 direction that bends (Figure 7(c)). In-between these two extremes, the deformation appears to be a combination of these bending mechanisms (Figure 7(b)). The reason for the variation in deformation can be readily explained in relation to the fact that when d 6 ¼ (l 2 þ t 1 )/2 (or equivalently v 6 ¼ (1 þ g)/2), the vertical ligament is split in unequal parts. As a consequence, the shorter of the two is likely to be able to bend less in response to the three-point force system. On similar lines, increasing the thickness of the ligaments alters the relative stiffness making it harder to deform the structure. In this regard, it can be mentioned that having different t 1 and t 2 is also expected to alter the relative stiffness. However, its effect on the deformation mechanism is still being studied.
Once the main deformation mechanisms have been identified, they can be related to the observed variation illustrated in Figure 6. It can be noted that at small offsets, when the meta-tetrachiral character dominates, the incremental Poisson's ratio has the tendency of attaining higher values for small deformations, i.e. for e 1 close to zero. On increasing the axial strain this trend seems to be inverted for small values of l 1 leading to more negative values for 12 . However, the added advantage for this configuration decreased with increasing l 1 possibly due to a leverage effect which is more efficient in the reentrant mode than in the meta-tetrachiral configuration. Hence, the optimal configuration to maximize auxeticity depends not just on the offset but also on the length of the horizontal ligament.
Additional simulations were also carried out to investigate the effect of changing the Poisson's ratio of the material or using a completely different material on the results presented above. Analysis of the results, which are found in the Supplementary Material, indicates that the deviations are relatively small. Slight differences emerged mainly for smaller offsets and for larger values of l 1 . It is most likely that the observed deviations are due to differences in the high-strain nonlinear behavior of the specific structural material. In fact, differences were mainly observable for aluminum and steel which had a distinct stress-strain relation. This suggests that the observed auxetic behavior was not much affected by changes in the mechanical properties of the material being used. (See the Supplementary Material for details).
In order to extend further the investigation, for the main analysis two prototypes were manufactured and mechanically testing as described in Section 2.3. The experimental results together with those obtained from the numerical simulations are illustrated in Figure 8. As can be noted, Structure 1 was confirmed to deform like a reentrant hexachiral honeycomb for large deformations while Structure 2 like a meta-tetrachiral. The deformed shapes observed in the experiment are in general agreement with those obtained numerically. The same can be said of the variations of the strains determined from the two methods.
For the purpose of obtaining the incremental Poisson's ratio, a second order polynomial was fitted to the experiential data of Structures 1 and 2, while interpolation using polynomials of order three was used on the FEA results. This provided a way of estimating 12 throughout the deformation process. The results, illustrated in Figure 9, indicate that there is relatively good agreement. This is particularly the case of Structure 1, most likely due to the rather smooth variation observed in the experimental data (Figure 8(g)).
The differences in the result obtained, particularly those of Structure 2, should be interpreted with care. The fitting of the polynomial to the experimental measurements, can never be perfect and, in fact, can only be considered to be an approximation to the observed behavior. In the case of Structure 2, the measured data is somewhat noisy. This makes it all that important to use some form of smoothing to get the average behavior. At the same time, it makes it hard to attain an optimal polynomial fit. In fact, increasing the order of the polynomial was observed to lead to overfitting, with the curve for the experimental results of the incremental Poisson's ratio increasingly oscillating around that representing the FEA. These aspects need to be considered when interpreting the results.
The reason for the noisy variation of the experimental data of Structure 2 can also be of interest. One possible explanation is that this is due to the fact that the structure was more prone to shear collapse than Structure 1. In fact, Structure 2 experienced a structural failure at an axial strain of around À0.02 (as determined visually from the e 2 against e 1 graph) while Structure 1 experienced a shear collapse when e 1 was around À0.05. Yet, both structures had similar incremental Poisson's ratios when they failed. However, this was not the case of Structures 3 and 4 (see Supplementary Material for details) which at the point when they failed had a value of 12 that was larger than that exhibited by Structures 1 and 2. This shows how small modification can affect the minimum value of the incremental Poisson's ratio attained at large deformations. On the other hand, it can be  mentioned that good agreement was attained between the numerical and experimental results of Structures 3 and 4, which gives further confidence in the analysis.
The observed tendency of the structures to experience shear collapse at some point is one of the main limitations of the square and rectangular grids being considered here. Given that relatively large deformations are involved, this eventuality was to be expected. It should be further noted that the same type of collapse is expected to occur if the system is loaded in an off-axis direction. Based on these considerations, the failure modes of the structures need to be studied to greater depth before it can be used for practical applications. The investigation should also look into the strain distribution inside the structure so as to gain further insight on the behavior of the system. It can be expected that with improved design these limitations can be mitigated. However, such additional investigations were beyond the scope of the work which is focused on illustrating the concept and the potential of the methodology.
The results being reported here are very significant as it does not seem that harnessing large deformation to produce auxetic behavior was ever considered. A potential application immediately follows. Given that two beneficial features of auxetic structures are their ability to resist indentation [3][4][5] and better energy absorption characteristics [6], they have been proposed and investigated for possible use as core materials in sandwich panels to mitigate the effects of impact loading [22][23][24][25]. Yet, the reported variation of the Poisson's ratio with the applied strain, which is expected to be large in the case of impact loading, can limit the benefits of using an auxetic material. On the other hand, it can be inferred from Figures 4-6 that the Poisson's ratio for the proposed system can actually become more negative with strain. This could lead to the development of structures that are able to improve their energy absorption and indentation resistance characteristics with increasing strain rather than simply losing it or keep it constant. The added advantages of this paradigm shift can be noteworthy.
It should also be mentioned that the regular square grid (Figure 1(a)) is one of the designs that has emerged as a prominent candidate for use in sandwich structures [85][86][87] intended to mitigate impacts. In these investigations, the configuration considered had the direction of the holes perpendicular to that of the plates. These studies have shown that with this alignment the square grid is able to exhibiting high energy absorbing, high crushing strength, and substantial stretch resistance [85,87]. On the other hand, this work has suggested that using a nonzero offset between successive horizontal ligaments might introduce added benefits. It also indicates that adopting a different alignment, whereby the plates are positioned along the vertical ligaments allowing the rectangular grid to have an auxetic behavior at larger deformation, might also lead to a sandwich panel which is effective in mitigating impact loading. This configuration is of specific relevance in view of that fact that a number of studies indicated that rectangular honeycombs perform well in forced convective heat transfer applications [71,88,89]. Hence, the proposed design might allow the rectangular honeycombs to be used as a protection to heated bodies while at the same time form part of the cooling element. One such deployment could be their use as an outer surface layer of space vehicles in order provide cooling during reentry [88] while at the same time provide additional protection against any accidental impact with, for example, space debris. These considerations might open the way for the rectangular honeycombs to be employed in multifunctional applications. It is however important to remark that these reflections require further investigations to assess their validity, particularly with respect to the effects of thermal loading and high velocity impact on the mechanical response of these metamaterials.

Conclusion
This work investigates the in-plane auxetic potential of rectangular honeycombs. It was shown that, for the configurations considered, at low axial strains rectangular honeycombs exhibit a small positive incremental Poisson's ratio. However, as the axial strain increases, the incremental Poisson's ratio becomes negative at some point. The behavior of the system was found to depend on the offset as well as on the relative lengths and the thickness of the ligaments. Careful examination of the simulated deformations revealed that for small values of the offset, the ligaments of the rectangular honeycomb bend like those of the meta-tetrachiral. On the other hand, for large values of the offset, the rectangular honeycomb behaves like the reentrant hexagonal honeycomb. These results were corroborated through the mechanical testing of 3D printed prototypes. Hence, even though the starting point was a normal rectangular grid, it was shown that small modifications can lead to a highly altered behavior in terms of deformation and hence mechanical properties.
A particular feature that is of relevance is the fact that for the studied structures the Poisson's ratio can become more negative as the strain increases. The increase in auxetic character with large deformations can be of particular relevance in applications requiring the mitigation of impacts as it might lead to improved performance. Also, of note is the fact that rectangular girds can find usage as multifunctional systems due to their good performance in heat transfer under forced convection conditions. It is thus hoped that this work will spur further research in the behavior and potential use of rectangular honeycombs.