Idealized 3D Auxetic Mechanical Metamaterial: An Analytical, Numerical, and Experimental Study

Mechanical metamaterials are man-made rationally-designed structures that present unprecedented mechanical properties not found in nature. One of the most well-known mechanical metamaterials is auxetics, which demonstrates negative Poisson’s ratio (NPR) behavior that is very beneficial in several industrial applications. In this study, a specific type of auxetic metamaterial structure namely idealized 3D re-entrant structure is studied analytically, numerically, and experimentally. The noted structure is constructed of three types of struts—one loaded purely axially and two loaded simultaneously flexurally and axially, which are inclined and are spatially defined by angles θ and φ. Analytical relationships for elastic modulus, yield stress, and Poisson’s ratio of the 3D re-entrant unit cell are derived based on two well-known beam theories namely Euler–Bernoulli and Timoshenko. Moreover, two finite element approaches one based on beam elements and one based on volumetric elements are implemented. Furthermore, several specimens are additively manufactured (3D printed) and tested under compression. The analytical results had good agreement with the experimental results on the one hand and the volumetric finite element model results on the other hand. Moreover, the effect of various geometrical parameters on the mechanical properties of the structure was studied, and the results demonstrated that angle θ (related to tension-dominated struts) has the highest influence on the sign of Poisson’s ratio and its extent, while angle φ (related to compression-dominated struts) has the lowest influence on the Poisson’s ratio. Nevertheless, the compression-dominated struts (defined by angle φ) provide strength and stiffness for the structure. The results also demonstrated that the structure could have zero Poisson’s ratio for a specific range of θ and φ angles. Finally, a lightened 3D re-entrant structure is introduced, and its results are compared to those of the idealized 3D re-entrant structure.


Introduction
The Poisson's ratio is defined as the ratio of the negative value of transverse strain to the longitudinal strain of a body subjected to uniaxial stress [1]. The Poisson's ratio is usually assumed to be positive (v > 0) for natural materials, which means that the material contracts transversely when stretched longitudinally or expands transversely when compressed longitudinally [2]. In contrast, there are some materials that reveal the opposite behavior of common materials and present negative Poisson's ratio, which means that their transverse dimensions increase (or decrease) when they are subjected to axial tensile (or compressive) loads [3]. These materials have been given different names in the literature, including anti-rubber, dilational [4], and auxetic materials [5]. The term "auxetic" is currently the most widely accepted name used for materials with a negative Poisson's ratio. Auxetic behavior has been observed in a few natural materials such as cristobalite (polymorphic silicones) [6], metals [7], zeolites [8], silicates [9,10], and some biological tissues such as cancellous bone [11], tendons [12], and some animal skins [13,14].
Mechanical metamaterials are man-made materials with unusual mechanical properties and special functionalities not found in nature [15][16][17][18][19]. Negative Poisson's ratio could be considered as one of such special mechanical properties [20][21][22][23][24][25][26][27][28][29][30][31]. Due to recent advances in additive manufacturing techniques (3D printing), which provides the possibility of fabrication of rationally-designed materials with the desired accuracy in nano/microscale [32][33][34][35], the studies concentrated on design and analyzing auxetic metamaterials have been increasing exponentially during the last few years [36][37][38][39][40]. Different geometries with different deformation mechanisms have been introduced for auxetic structures. One of the first auxetic structures was proposed by Roberts [21], the first auxetic foam was presented by Lakes [41], and the first 2D molecular auxetic models (spontaneously forming auxetic phases) were studied by Wojciechowski [22,23]. The term auxetic in the scientific literature was introduced by Evans [24]. The re-entrant hexagonal honeycomb structure is one of the most well-known open-cell auxetic structures [42]. In several studies [43][44][45][46], the effect of relative density and re-entrant angle variations on mechanical properties such as Poisson's ratio and elastic modulus of typical 2D re-entrant structure in the elastic and plastic range have been analyzed for predicting the behavior of the structure. To improve the mechanical properties such as strength and auxeticity of the re-entrant structures and for answering the new requirements in the industry, some researchers have introduced and built novel 3D re-entrant structures inspired by typical 2D re-entrant unit cells and have established analytical and numerical solutions for them [47][48][49][50].
The unit cell presented by Yang et al. in 2015 [50] (Figure 1a) is one of the successful 3D auxetic structures, which has been solved analytically in two directions and evaluated by experimental and numerical results. In another study, Wang et al. (2016) [51] proposed an interlocking assembly manufacturing technique for fabricating the re-entrant lattice structure proposed by   [50] but with a shifted unit cell (Figure 1b), which showed great advantages over additive manufacturing technique. Chen et al. (2017) and Xue et al. (2020) [52,53] applied some modifications to the noted 3D re-entrant auxetic structure by adding some struts to enhance the mechanical properties of the lattice structure, and in particular, to increase the stiffness of this auxetic structure (Figure 1c,d).
Although the 3D reentrant structures studied in the noted works are successful auxetic structures for many applications, they still lack some geometrical characteristics, which could raise problems in many applications. First of all, even though the noted 3D re-entrant unit cells are symmetrical in the lateral directions, but their geometry is different in the vertical (i.e., loading) direction. Therefore, they exhibit dissimilar mechanical properties in the three main orthogonal directions, which means that they are not truly 3D auxetic structures. Second, the struts at the edge of the noted unit cells are shared by adjacent unit cells in the lattice structure. Since the auxetic structures are useful in many applications (such as biomedical implants [54,55], sandwich panels [56], etc.) that require graded microstructures, the noted shared struts cause huge complexities and difficulties in constructing truly graded lattice structures. Hence, there remains a gap to propose an ideal unit cell that could overcome these deficiencies of typical 3D re-entrant structure.
In this study, a specific type of mechanical metamaterial, namely, an idealized 3D reentrant structure [41,57] constructed of three strut types, one type loaded purely axially and two types loaded both flexurally and axially, which are inclined and are spatially defined by angles θ and ϕ (see Figure 2), is studied to obtain analytical relationships for its mechanical properties. Analytical solutions based on two well-known beam theories namely Euler-Bernoulli and Timoshenko theories are derived. Moreover, two numerical methods one based on volumetric elements and one based on beam elements are implemented for validating the results. A definition for strut effective length has been presented in order to simulate what the struts experience in real-life conditions, and analytical results have been presented based on the effective length. Moreover, a lightened auxetic unit cell extracted from the general re-entrant unit cell (presented in Figure 2) is introduced and its mechanical properties are compared to those of the ideal re-entrant structure. Several specimens have also been additively manufactured and tested under compressive loading conditions. Finally, the results of the derived analytical and constructed numerical models and the manufactured specimens are compared to one another. In this study, a specific type of mechanical metamaterial, namely, an idealized 3D reentrant structure [41,57] constructed of three strut types, one type loaded purely axially and two types loaded both flexurally and axially, which are inclined and are spatially defined by angles and (see Figure 2), is studied to obtain analytical relationships for its mechanical properties. Analytical solutions based on two well-known beam theories namely Euler-Bernoulli and Timoshenko theories are derived. Moreover, two numerical methods one based on volumetric elements and one based on beam elements are implemented for validating the results. A definition for strut effective length has been presented  [52], and (d) Yingying Xue et al. (2020) [53].
in order to simulate what the struts experience in real-life conditions, and analytical results have been presented based on the effective length. Moreover, a lightened auxetic unit cell extracted from the general re-entrant unit cell (presented in Figure 2) is introduced and its mechanical properties are compared to those of the ideal re-entrant structure. Several specimens have also been additively manufactured and tested under compressive loading conditions. Finally, the results of the derived analytical and constructed numerical models and the manufactured specimens are compared to one another.  Figure S1) consists of three types of struts with lengths l 1 , l 2 , l 3 and cross-sectional areas A 1 , A 2 , A 3 . The unit cell is fully defined with two unique angles θ and ϕ, the Type-I strut length, and the main cube size (the main cube is the hypothetical cube formed by connecting vertices D i where i = 1 to 8). Angle θ is defined as the angle between Type-II strut and the face diagonals of the main cube, and angle ϕ is defined as the angle between a Type-III strut and the corresponding side of the main cube ( Figure 2b). Each re-entrant unit cell with mentioned strut lengths (l 1 , l 2 , l 3 ) and characteristic angles (θ, ϕ) occupies a cubic volume with side length of 2(l 1 + l 3 cos ϕ − l 2 sin θ) and hence the volume of V total = 8(l 1 + l 3 cos ϕ − l 2 sin θ) 3 . Because the unit cell is made up of six struts of Type I, 24 struts of Type II, and 24 struts of Type III, the volume occupied by all struts inside the unit cell is V strut = 6πr 2 1 l 1 + 24πr 2 2 l 2 + 24πr 2 3 l 3 . Therefore, according to the basic definition of relative density µ = V strut V total , the relative density of a re-entrant unit cell is given by 3 (1)

Stiffness Matrix
In this subsection, analytical relationships for mechanical properties of the 3D reentrant unit cell such as elastic modulus, Poisson's ratio, and yield stress are derived as functions of elastic properties of the bulk material (E s , G s , σ ys , and v s ) and the geometrical dimensions of the unit cell. The re-entrant unit cell has symmetry with respect to its major planes (see Figure S2). Therefore, studying 1/8 of the unit cell is sufficient for obtaining the mechanical properties of the unit cell and, hence, the corresponding lattice structure. Due to the intrinsic symmetries in the geometry of the re-entrant unit cell and the applied boundary conditions on the unit cell, seven unique group of vertices can be recognized, which are denoted by letters A, B, C, D, E, F, and G, as shown in Figure 2a. For the sake of brevity and as the vertices named by similar letters demonstrate identical behaviors, in many cases, instead of referring to a vertex as X i , we refer to it as X (X representing A, B, C, D, E, F, and G). For instance, in many cases, we use the name E instead of E 1 .
If the origin of the coordinate system is considered to be located in the midst of the unit cell, each re-entrant unit cell has four vertical (XY, YZ, and two bisectors of XY and YZ) and one horizontal (XZ) symmetry planes (as shown in Figure S2). All the vertices located in the symmetry planes are only allowed to translate and rotate in their corresponding planes. Vertices A i and B i (for i = 1, 2) lie in all of the four vertical symmetry planes, so they are only allowed to translate in the Y direction, and they are constrained rotationally. Therefore, they create two degrees of freedom (DOFs) q 8 and q 7 at vertices A i and B i , respectively. Vertices C i (for i = 1 to 8) are located in one of the vertical symmetry planes (XY or YZ), and hence they are only allowed to translate horizontally and vertically and rotate in their corresponding planes, which results in three other DOFs (q 3 , q 4 , q 11 ). Each vertex D i (for i = 1 to 8) is located in a single vertical symmetry plane (one of two bisectors of YZ and YX), and therefore, it is only allowed to translate horizontally and vertically and rotate in the planes it is located in, which leads to three other DOFs (q 5 , q 6 , q 10 ). Vertices E i and F i (for i = 1 to 4) are located in one of two vertical symmetry planes (XY or YZ) and the horizontal symmetry plane XZ, so they are only allowed to translate in the X or Z directions without any rotations. This creates two other DOFs q 2 and q 1 at vertices E and F, respectively. Vertices G i are located in the horizontal symmetry plane XZ and in one of the vertical bisector symmetry planes (bisectors of XY and YZ). Therefore, these points cannot rotate in any directions and are only allowed to translate horizontally in the bisector of the noted planes, which leads to the final DOF q 9 . Overall, the system has a total number of 11 DOFs (Figure 2b).
In this study, two well-known elastic beam theories, namely, Euler-Bernoulli and Timoshenko beam theories, have been implemented for analytical analysis of the re-entrant unit cell. As the beam theories and material models used for the analytical analysis are linear, the overall deformation of each DOF in the system can be considered as the superposition of separate deformations caused by applying individual loads at each DOF. Therefore, the superposition principle could be implemented for obtaining the system of equations of the unit cell. In this approach, each DOF displaces autonomously from other DOFs, and when a DOF is displaced, the other DOFs are kept fixed. By solving the equilibrium equations for the vertices of the unit cell, the resultant forces are obtained. By considering this method, the system of equilibrium equations for the structure could be written as {Q} = [K]{q}. In this equation, {Q} is the force vector comprising the external forces acting on the DOFs, [K] is the stiffness matrix of the system, and the {q} is the displacement vector.
Because the derivation of the stiffness matrix elements is based on beam theories, the general deformation of a cantilever beam ( Figure 3) can be considered as the resultant of four distinct deformations at the free end of the beam, which are as follows: (a) a lateral displacement without any rotation, v; (b) a flexural rotation without any lateral displacement, θ; (c) a longitudinal elongation or contraction, u; (d) an axial twist, ϕ. Therefore, these points cannot rotate in any directions and are only allowed to translate horizontally in the bisector of the noted planes, which leads to the final DOF . Overall, the system has a total number of 11 DOFs (Figure 2b).
In this study, two well-known elastic beam theories, namely, Euler-Bernoulli and Timoshenko beam theories, have been implemented for analytical analysis of the re-entrant unit cell. As the beam theories and material models used for the analytical analysis are linear, the overall deformation of each DOF in the system can be considered as the superposition of separate deformations caused by applying individual loads at each DOF. Therefore, the superposition principle could be implemented for obtaining the system of equations of the unit cell. In this approach, each DOF displaces autonomously from other DOFs, and when a DOF is displaced, the other DOFs are kept fixed. By solving the equilibrium equations for the vertices of the unit cell, the resultant forces are obtained. By considering this method, the system of equilibrium equations for the structure could be written as { } = [ ]{ }. In this equation, { } is the force vector comprising the external forces acting on the DOFs, [ ] is the stiffness matrix of the system, and the { } is the displacement vector.
Because the derivation of the stiffness matrix elements is based on beam theories, the general deformation of a cantilever beam ( Figure 3) can be considered as the resultant of four distinct deformations at the free end of the beam, which are as follows: (a) a lateral displacement without any rotation, ; (b) a flexural rotation without any lateral displacement, ; (c) a longitudinal elongation or contraction, ; (d) an axial twist, . The corresponding forces and moments required to create each of the pure abovementioned deformations are demonstrated in Figure 4 for the Euler-Bernoulli and Timoshenko beam theories. Due to frequent repetition of some terms demonstrated in Figure 4 throughout the whole manuscript, the terms , , , , and are used, which are different for Euler-Bernoulli and Timoshenko beam theories and are listed in Table 1. The corresponding forces and moments required to create each of the pure abovementioned deformations are demonstrated in Figure 4 for the Euler-Bernoulli and Timoshenko beam theories. Due to frequent repetition of some terms demonstrated in Figure 4 throughout the whole manuscript, the terms S i , W i , T i , V i , and U i are used, which are different for Euler-Bernoulli and Timoshenko beam theories and are listed in Table 1.

Term
Euler-Bernoulli Theory Timoshenko Theory    To explain the general procedure for deriving the stiffness matrix elements, we could consider a strut PQ with an arbitrary orientation with respect to the global X, Y, and Z coordinate system (aligned with a unit vector u PQ ), which is deformed by a specific DOF q i at its free end. The schematic representation of such a strut and its deformation due to the DOF is illustrated in Figure 5. As shown in the figure, DOF q i causes a longitudinal and a transverse deformation in the strut PQ. The vector of longitudinal deformation of strut PQ can be obtained by projecting the vector q i on the strut PQ unit vector, i.e., q i ·u PQ u PQ . Moreover, the vector of transverse deformation of strut PQ can be obtained as q i − q i ·u PQ u PQ . By multiplying the longitudinal deformation vector of strut PQ by the resultant force value of PQ (i.e., S i ), the resultant longitudinal force vector for strut PQ can be obtained as q i ·u PQ u PQ S i , and by multiplying the transverse deformation vector of strut PQ by the resultant transverse force value of PQ (i.e., T i ), the resultant transverse force vector for this strut can be obtained as q i − q i ·u PQ u PQ T i . The reaction Materials 2021, 14, 993 8 of 36 moment vector, created by the transverse deformation of strut PQ can be found by the crossproduct of the unit vector of PQ (i.e., u PQ ) and the resultant moment vector of the beam (i.e., q i − q i ·u PQ u PQ V i ), which yields − u PQ × q i − q i ·u PQ u PQ V i . The values of S i , T i , and V i for Timoshenko and Euler-Bernoulli beam theories are listed in Table 1. It is worth noting that the procedure of obtaining the forces and moments for a lateral rotational DOF is very similar to the above-mentioned procedure, and it could be easily performed by replacing the terms S i , T i , and V i with the terms W i , V i , and U i , respectively. the resultant force value of PQ (i.e., ), the resultant longitudinal force vector for strut PQ can be obtained as .
, and by multiplying the transverse deformation vector of strut PQ by the resultant transverse force value of PQ (i.e., ), the resultant transverse force vector for this strut can be obtained as − · . The reaction moment vector, created by the transverse deformation of strut PQ can be found by the crossproduct of the unit vector of PQ (i.e., ) and the resultant moment vector of the beam (i.e., − · ), which yields − × − · . The values of , , and for Timoshenko and Euler-Bernoulli beam theories are listed in Table 1. It is worth noting that the procedure of obtaining the forces and moments for a lateral rotational DOF is very similar to the above-mentioned procedure, and it could be easily performed by replacing the terms , , and with the terms , , and , respectively.
In the following, the detailed procedure of derivation of the stiffness matrix elements for three distinct DOFs ( , , and ) are presented. The derivation procedure of other DOFs is similar to the ones presented in the paper, but they are provided in Section S.1.2 and Figure S6.  , where is the unit vector parallel to the first DOF, and is the unit vector parallel to strut . The same system of notation will be used in the rest of the document. By multiplying the elongation vector of strut EF by the resultant force parameter of strut EF (i.e., ), the resultant force vector can be obtained as . .
The equilibrium of forces at point in the direction, i.e., the Z direction, provides In the following, the detailed procedure of derivation of the stiffness matrix elements for three distinct DOFs (q 1 , q 2 , and q 10 ) are presented. The derivation procedure of other DOFs is similar to the ones presented in the paper, but they are provided in Section S.1.2 and Figure S6.
(a) First DOF: Here, the elements of the first column of the stiffness matrix are derived by applying the unit displacement q 1 = 1 and by keeping all the other DOFs fixed (i.e., = 0). This deformation type displaces point F 1 horizontally in the Z direction. This deformation only causes elongation in strut EF and it does not affect the other struts. The vector of elongation of strut EF can be obtained by projecting the vector q 1 on strut EF, i.e., where q 1 is the unit vector parallel to the first DOF, and u E 1 F 1 is the unit vector parallel to strut E 1 F 1 . The same system of notation will be used in the rest of the document. By multiplying the elongation vector of strut EF by the resultant force parameter of strut EF (i.e., S 1 ), the resultant force vector can be obtained as The equilibrium of forces at point F 1 in the q 1 direction, i.e., the Z direction, provides the K 11 element of the stiffness matrix ( Figure 6), which is calculated as Materials 2021, 14, x FOR PEER REVIEW 10 of 36 Figure 6. Deformation of strut EF due to = 1.
(b) Second DOF: = 1 In this subsection, the elements of the second column of the stiffness matrix are derived. Therefore, we set = 1, and the other DOFs are set to zero. This deformation type displaces point E horizontally in the Z direction, which causes strut EF to have a pure contraction and strut ED to have contraction accompanied by transverse displacement. The vector of contraction of strut EF can be obtained by projecting vector on strut EF, i.e., .
. The vector of contraction of strut ED can be obtained by projecting vector on strut ED, i.e., .
, and the vector of transverse displacement of vertex E of strut ED can be calculated as − .
. By multiplying the contraction vector of strut EF by the resultant force parameter of strut EF (i.e., ), the resultant force vector of strut EF can be obtained as .
. Moreover, by multiplying the elongation vector of strut ED by the resultant force parameter of strut ED (i.e., ), the resultant force vector can be obtained as .
. Finally, by multiplying the transverse deformation vector of strut ED by the resultant transverse force parameter of strut ED (i.e., ), the resultant transverse force vector of this strut can be obtained as − .
. It is worth noting that the resultant forces obtained for deformation of struts ED at point E must be multiplied by four because there are four similar struts ED connected to point E. The equilibrium of forces at point E in the direction (i.e., the Z direction), Figure 7b It must be noted that the forces demonstrated in Figure 6 are the only forces that must act on the beam in order for it to deform in the way it is shown. The force the beam applies at each vertex is in the opposite direction of what is shown in Figure 6. The same holds true for all the next derivations. The equilibrium of forces at point E 1 in the q 2 direction (i.e., in the Z direction) provides the K 12 element of the stiffness matrix, which is calculated as Because the displacement of DOF q 1 does not have any effect on other DOFs (q 3 , q 4 , q 5 , q 6 , q 7 , q 8 , q 9 , q 10 , q 11 ), the remaining elements of the first column of the stiffness matrix are zero.
In this subsection, the elements of the second column of the stiffness matrix are derived. Therefore, we set q 2 = 1, and the other DOFs are set to zero. This deformation type displaces point E horizontally in the Z direction, which causes strut EF to have a pure contraction and strut ED to have contraction accompanied by transverse displacement. The vector of contraction of strut EF can be obtained by projecting vector q 2 on strut EF, i.e., The vector of contraction of strut ED can be obtained by projecting vector q 2 on strut ED, i.e., q 2 ·u E 1 D 1 u E 1 D 1 , and the vector of transverse displacement of vertex E of strut ED can be calculated as By multiplying the contraction vector of strut EF by the resultant force parameter of strut EF (i.e., S 1 ), the resultant force vector of strut EF can be obtained as Moreover, by multiplying the elongation vector of strut ED by the resultant force parameter of strut ED (i.e., S 2 ), the resultant force vector can be obtained as Finally, by multiplying the transverse deformation vector of strut ED by the resultant transverse force parameter of strut ED (i.e., T 2 ), the resultant transverse force vector of this strut can be obtained as It is worth noting that the resultant forces obtained for deformation of struts ED at point E must be multiplied by four because there are four similar struts ED connected to point E. The equilibrium of forces at point E in the q 2 direction (i.e., the Z direction), Figure 7b, provides the element K 22 of the stiffness matrix as follows: as .
. Finally, by multiplying the transverse deformation vector of strut ED by the resultant transverse force parameter of strut ED (i.e., ), the resultant transverse force vector of this strut can be obtained as − .
. It is worth noting that the resultant forces obtained for deformation of struts ED at point E must be multiplied by four because there are four similar struts ED connected to point E. The equilibrium of forces at point E in the direction (i.e., the Z direction), Figure 7b The presence of reaction forces and moments at point D due to leads to non-zero values of , , and . By projecting the resultant contraction force vector . and the transverse deformation force vector − . of strut ED in the and directions, and solving the equilibrium of forces at point D in these two directions, the elements and of the stiffness matrix could be obtained. It is important to note that the force vectors at point D have been multiplied by two (as two ED struts are connected to point D): To obtain the element of the stiffness matrix, the equilibrium of moments in the direction at point D must be considered. The reaction moment vector, created by transverse deformation of point E of strut ED can be obtained by the cross-product of unit vector of strut ED, , and the resultant moment vector of the beam at point D, The presence of reaction forces and moments at point D due to q 2 leads to nonzero values of K 25 , K 26 , and K 210 . By projecting the resultant contraction force vector 5 and q 6 directions, and solving the equilibrium of forces at point D in these two directions, the elements K 25 and K 26 of the stiffness matrix could be obtained. It is important to note that the force vectors at point D have been multiplied by two (as two ED struts are connected to point D): To obtain the element K 210 of the stiffness matrix, the equilibrium of moments in the q 10 direction at point D must be considered. The reaction moment vector, created by transverse deformation of point E of strut ED can be obtained by the cross-product of unit vector of strut ED, u E 1 D 1 , and the resultant moment vector of the beam at point D, This vector must be multiplied by two in the equilibrium equation as two struts ED are connected to point D. Finally, by solving the equilibrium of moments in the q 10 direction, the element K 210 of the stiffness matrix could be derived: (c) The 10th DOF: This DOF is a complex deformation in the structure because it applies rotation at vertex D, which is connected to six struts (2 × DC, 2 × DE, DG, and DB). By projecting vector q 10 on the noted struts, the torsion applied to each of the struts can be found (q 10 ·u C 1 D 1 , q 10 ·u E 1 D 1 , q 10 ·u B 1 D 1 , q 10 ·u G 1 D 1 ). Afterward, by multiplying the torsional parameter of each strut (W 2 and W 3 ) by the projected twist, the torque vectors could be obtained as q 10 The extent of flexural rotation of struts DC, DE, DG, and DB can be obtained by subtracting the projected twist vector from vector q 10 , which yields q 10 − q 10 Afterward, by multiplying the resultant bending parameter of each strut (i.e., U 2 and U 3 ) to the noted terms, the lateral force vectors could be calculated. It is worth noting that the terms related to struts CD and ED must be multiplied by two because two of such struts are connected to point D. Finally, by projecting all the resultant axial and lateral moments in the q 10 direction and writing the equation of moment equilibrium in this direction, the element K 1010 can be found: The element K 1011 could be obtained by solving the equilibrium of moments at vertex C. By projecting the reaction torsion and bending moment vectors created by the rotation of vertex D of strut CD in the q 10 direction and solving the equations of moment equilibrium in the same direction, the element K 1011 could be calculated as follows:

Stiffness Matrix Derivation
Finally, using the stiffness matrix elements obtained in Section 2.1.2 and also by considering the loading condition (Q 8 = 2F), the stiffness matrix and system of equations can be constructed as where F is the external load applied on the unit cell.  and multiplying it by the force vector. By having the DOF displacements as functions of the applied external force, unit cell elastic mechanical properties such as elastic modulus, Poisson's ratio, and yield stress can be obtained. The elastic modulus of the unit cell can be calculated according to the basic definition of elastic modulus as follows: where, F UC , L UC , A UC , δ UC , and E UC are the applied external load, unit cell length, unit cell cross-sectional area, unit cell displacement, and elastic modulus of the re-entrant unit cell, respectively. By considering F UC = F, δ UC = 2q 8 , L UC = 2(l 1 + l 3 cos ϕ − l 2 sin θ), and A UC = L UC 2 , the elastic modulus of the re-entrant unit cell can be obtained as Calculation of q 8 from solving Equation (10) and inserting it in Equation (13) provides the elastic modulus relationship of the re-entrant unit cell. Because the final relationships for mechanical properties of the unit cell are very lengthy, the elastic modulus relationship for a specific case (θ = ϕ = 0 • , r 1 = r 2 = r 3 = 0.14b, and l 1 = b) is presented below: Poisson's ratio of the unit cell ν UC can be found by − q 1 q 8 . Hence, by calculating q 1 and q 8 from Equation (10) and substituting them in this relationship, the Poisson's ratio formula for the unit cell can be derived. The Poisson's ratio relationship of the above-mentioned specific case (θ = ϕ = 0 • , r 1 = r 2 = r 3 = 0.14b, and l 1 = b) can be obtained as follows: Evaluating the von Mises stress contours of the FE models with several relative densities showed that vertices B of struts BD are the most critical points (the points with the maximum stress) in the structure. The yielding of the unit cell as a whole occurs when a point in the unit structure reaches the yield stress of the bulk material σ ys . Therefore, the yield stress of the re-entrant unit cell can be obtained by assuming that the structure yielding occurs when the stress at point B of struts BD reaches the yield stress of the bulk material. Because strut BD is influenced by both bending moment and axial force (M DB , F DB ), the maximum stress (σ max ) in point B of this strut could be calculated as follows: where r 2 , I 2 , and A 2 are the maximum distance of the points on the strut cross section from the neutral axis, area moment of inertia, and cross-sectional area of strut DB, respectively. To obtain the bending moment and axial force in strut DB, all of the lateral deformation and axial elongation of this strut applied by DOFs should be considered. DOFs q 5 , q 6 , and q 7 apply axial elongation and lateral displacement to ends B and D of this strut. Moreover, DOF q 10 applies twist to point D. Eventually, the axial forces at point B of strut BD can be obtained as follows: Similarly, the bending moments at point B of strut DB can be obtained as M DB = (−q 5 cos θ + q 7 cos θ + q 6 sin θ)V 2 + (q 10 )U 2 (18) By substituting F DB and M DB from Equations (17) and (18) into Equation (16) and inserting q 5 , q 6 , q 7 , and q 10 from solution of the inverse of Equation (10), the maximum local yield stress in strut DB can be obtained. For the specific case introduced above (θ = ϕ = 0 • , r 1 = r 2 = r 3 = 0.14b, and l 1 = b), the relationship can be simplified as follows: Because the applied stress equal to the structure yield stress of the re-entrant unit cell σ y,UC causes the maximum local yield stress σ max at point B, reaching the bulk material yield stress σ y,s , the normalized yield stress of the unit cell could be obtained by a crossmultiply of structure mean stress σ UC = F UC A UC and σ max on the one hand and the relevant yield stresses on the other hand, which is Because the idealized 3D re-entrant structure includes vertical struts (strut AB), the unit cell can be vulnerable to buckling under compressive loading. Strut AB in the unit cell structure can be considered as a column under a pure compressive load, and the critical elastic buckling force of the column can be obtained using the classic Euler's buckling theory as follows: By considering strut AB as a clamped column based on the periodic boundary conditions, the theoretical K value of the column will be equal to 1 and the critical load of an idealized 3D re-entrant unit cell under compressive load F can be obtained as Because the total load applied to a single unit cell and to strut AB is identical (load F), by considering the 4(l 1 + l 3 cos ϕ − l 2 sin θ) 2 as the cross-sectional area of the unit cell, the critical buckling stress could be obtained as follows:

Effective Lengths
One of the important considerations in deriving the analytical relationships of mechanical properties of structures is that the struts are connected rigidly to each other at the vertices of structures. Ideally, the struts are connected to each other from their very endpoints (vertices demonstrated in Figure 2a), and they do not interact with each other from any other point. However, in practice, when several struts approach a vertex, as the cross-section area of the struts is non-zero, some volumes overlap with each other at the vertices, which in turn decreases the effective lengths of struts participating in the deformation of the structure. This phenomenon becomes of much more significance in auxetic lattices because they have very acute angles at the vertices (see, for example, Figures 2 and 8). In this study, the minimum distances between the joints were considered as the effective length of the struts. This definition for effective length was in accordance with what was observed in the stress and strain distribution of the structure obtained from the finite element (FE) results. In the following, the analytical derivation of effective lengths has been presented.
The effective length , can be calculated by subtracting the intersection length between Type-I and Type-II struts from as follows: The effective length of Type-II strut is calculated by measuring the distance between the intersection of two Type-II struts at point D (struts DB and DE) and the intersection of Type-I and Type-II struts at point B (struts AB and BD), and it can be calculated by subtracting the intersecting lengths from as where is the angle between two Type-II struts. Similarly, the overlapping of Type-III strut occurs in the intersection of Type-II and Type-III struts at point D (struts DC and DG) and the intersection of two Type-II struts at point C. The effective length can then be calculated by subtracting the intersection lengths from as where is the angle between Type-II and Type-III struts. By substituting these effective lengths into the stiffness matrix and solving the system of Equation (10), analytical relationships for mechanical properties of the 3D re-entrant unit cell can be obtained.
It is worth noting that the analytical relationships for angles and can be obtained as follows: = 2 cos + √2 cos sin + √2 sin sin 2 = cos ( cos + 2 sin ) In this study, the minimum distances between the joints were considered as the effective length of the struts. This definition for effective length was in accordance with what was observed in the stress and strain distribution of the structure obtained from the finite element (FE) results. In the following, the analytical derivation of effective lengths has been presented.

Experimental Tests
The effective length l 1,e f f can be calculated by subtracting the intersection length between Type-I and Type-II struts from l 1 as follows: The effective length of Type-II strut is calculated by measuring the distance between the intersection of two Type-II struts at point D (struts DB and DE) and the intersection of Type-I and Type-II struts at point B (struts AB and BD), and it can be calculated by subtracting the intersecting lengths from l 2 as where α 2 is the angle between two Type-II struts. Similarly, the overlapping of Type-III strut occurs in the intersection of Type-II and Type-III struts at point D (struts DC and DG) and the intersection of two Type-II struts at point C. The effective length can then be calculated by subtracting the intersection lengths from l 3 as where α 1 is the angle between Type-II and Type-III struts. By substituting these effective lengths into the stiffness matrix and solving the system of Equation (10), analytical relationships for mechanical properties of the 3D re-entrant unit cell can be obtained.

Experimental Tests
A Quantum Generous 3D printer (Persia Company, Iran) with poly-lactic acid (PLA) filaments was used for manufacturing the test specimens. Because the mechanical properties of PLA could be different after the printing process as compared to its properties in the original filament form, in order to obtain the mechanical properties of the bulk material, three types of cylindrical specimens were manufactured and tested under compressive loading. The cylinders were printed in the horizontal (0 • ), vertical (90 • ), and inclined (45 • ) orientations with respect to the build plate direction. Four specimens were manufactured for each cylinder type (12 specimens in total). All the cylindrical specimens had a nominal length of 20 mm and a diameter of 15 mm.
Three sets of re-entrant unit cell specimens with five different relative densities µ = 0.1, 0.15, 0.2, 0.25, and 0.3 were manufactured (15 specimens in total), see Figures S4 and S5. The dimensions of all the re-entrant unit cell specimens were 7.71 × 7.71 × 7.71 cm 3 . All the specimens (unit cells and cylinders) were manufactured with a layer thickness of 100 µm and infill density of 100%. Because the forces required for collapsing the manufactured specimens were in the range of 750 N to 3700 N, to avoid the inaccuracies caused by the deformation of the test machine itself, all the tests were performed with a large compression test machine (150 kN) for which the deformation in the machine was negligible. More information on the effect of deformation of the mechanical test bench can be found in [59]. The static compression tests were performed using an STM-20 (STANTAM, Iran) mechanical testing bench with a 20 kN load cell (Figure 8). The displacement rate was set to 0.1 mm/min. The normalized elastic modulus and normalized yield stress of specimens were obtained from load-displacement graphs, and the Poisson's ratio of the unit cells were extracted from the side-view images.

Numerical Modeling
Two FE models were developed in ANSYS software (Pennsylvania, United States) to calculate the mechanical properties of the re-entrant lattice structures. Because analytical formulas of mechanical properties of unit cells have been derived based on Euler-Bernoulli and Timoshenko beam theories, in the first FE model type, Timoshenko beam elements (3D Quadratic finite strain beam, i.e., beam element Type 189 in ANSYS) were used for discretizing the struts of the unit cell. Timoshenko beam theory takes into account the transverse shear deformation. Each strut was divided into five beam elements. The struts were considered to be rigidly connected to each other at the vertices. The linear elastic material properties of the PLA obtained from the experimental tests were adopted. The elastic modulus, Poisson's ratio, and yield stress of the beam elements were set to 1.93 GPa, 0.33, and 69 MPa, respectively. Periodic boundary conditions were applied to the FE models. The lowermost nodes of the unit cell at point A 2 Figure 2a) were fully constrained in all the transitional and rotational directions. The top node of the unit cell at point A 1 (Figure 2a) was displaced downward, but they were not also allowed to have any horizontal displacement or rotation in any direction ( Figure S3). The second FE model type was developed using volumetric elements. This FE model is expected to better take into account the overlapping effect at the joints on the overall mechanical properties of the structure. Five 3D models of the re-entrant unit cell with geometries almost identical to the geometry of the manufactured specimens (relative densities of 0.1, 0.15, 0.2, 0.25, and 0.3) were constructed. An adaptive fine volumetric mesh was implemented for discretizing the models using tetrahedral 10-node elements (SOLID187 elements in ANSYS) to have a fine mesh at the vertices and efficient larger mesh at the middle of the struts. The finest and most accurate options in ANSYS adaptive meshing, i.e., fine resolution (resolution = 6), slow transition, and fine span angle were chosen. The minimum and maximum element sizes for the lowest-density structure (µ = 0.1) were 0.3912 mm and 1.78 mm, respectively. The minimum and maximum element sizes for the highest-density structure (µ = 0.3) were 0.5786 mm and 2.45 mm, respectively. The number of elements in the lowest-density and the highest-density structures was 470,701 and 437,269, respectively. Similar to the beam elements, the elastic modulus and Poisson's ratio of 1.93 GPa and 0.33, respectively, were considered for the bulk material properties. To apply the loading and boundary condition, the lowermost face of the model (lower cross-sectional face of vertex of the strut in Figure 2) was fixed in all the directions and the top face of the model at the point was displaced downward. The material properties assigned to the volumetric elements were the same as those used in the first FE model.
In order to calculate the elastic modulus of the FE models, the sum of reaction force(s) at the lowermost node(s), F UC , due to the vertical displacement of the top node(s), δ UC , was measured and substituted in E UC = F UC L UC A UC δ UC , where A UC is the cross-sectional area of the FE model in the plane normal to the loading direction, and L UC is the length of the FE structure in the direction parallel to the loading direction.
For obtaining the Poisson's ratio, the elongation of the structure in the lateral direction (X or Z direction) was divided by the contraction of the structure in the Y direction (loading direction).
Eventually, for obtaining the numerical normalized yield stress, the maximum stresses (von Mises stress) in strut DB of the unit cell (σ max ) was extracted and inserted into the equation F UC A UC .σ max .

Results
In this section, all the obtained analytical, numerical, and experimental results are presented and compared to each other. To have a better understanding of the main factors influencing the overall performance of the auxetic lattice structure, the effect of each main geometrical parameter of the structure such as relative density, angle θ, and angle ϕ are studied individually. Moreover, a lightened unit cell, obtained by eliminating the struts l 3 from the general 3D re-entrant unit cell, is introduced and its mechanical properties are presented. More details on the results can be found in Section S.2 of the supplementary materials ( Figures S7-S27).

Effect of Relative Density
According to analytical relationships of mechanical properties of the re-entrant unit cell obtained in Section 2.1, elastic modulus, Poisson's ratio, and yield stress of the structure depend on seven geometrical main parameters, namely, r 1 , r 2 , r 3 , l 1 , l 2 or l 3 , θ, and ϕ. To validate the analytical relationships, the mechanical properties of a specific case has been considered in which all the struts have the same radii (r 1 = r 2 = r 3 ), angles θ and ϕ are considered to be equal to 22.5 • , and l 1 = b where b is shown in Figure 2b. It is worth noting that all the additively manufactured specimens have also been constructed according to the above-mentioned geometrical parameters. Effects of variations in the relative density of 3D re-entrant unit cell with respect to angles θ and ϕ is shown in Figure S7.
The analytical (based on Euler-Bernoulli and Timoshenko beam theories), numerical (based on beam and volumetric elements), and experimental results of this case are presented and compared in Figure 9. The main take-out of this figure is that both the Euler-Bernoulli and Timoshenko analytical models after consideration of the effective length, the numerical models based on volumetric elements, and the experimental data points have good agreement with respect to each other. On the other hand, the results of both the Euler-Bernoulli and Timoshenko analytical models without consideration of the effective length and the numerical model based on beam elements (rather than volumetric elements) are close to each other but quite far from the experimental data points. As it can be seen in Figure 9, the analytical models considering the overlapping effect and the numerical models made from volumetric elements have higher levels of stiffness and yield stress, and lower levels of Poisson's ratio in comparison with the analytical models without considering the overlapping effect and their equivalent numerical models, which are based on beam elements. the analytical models based on Euler-Bernoulli and Timoshenko beam theories are 14.84% and 0.91%, respectively. It is worth noting that the maximum negative Poisson's ratio for this case is very high and equal to −0.7791. These numerical/analytical differences for normalized yield stress are 13.76% and 5.32% for Euler-Bernoulli and Timoshenko beam theories, respectively. The maximum differences between the elastic modulus results of numerical models made from volumetric elements and the analytical models, which consider the overlapping effect, are 54.97% and 45.75% for Euler-Bernoulli and Timoshenko theories, respectively (Figure 9b). These differences for Poisson's ratio are 85.7% and 88.24% for Euler-Bernoulli and Timoshenko beam theories, respectively (Figure 9a). Finally, the noted differences for normalized yield stress are 32.31% and 6.85% for Euler-Bernoulli and Timoshenko beam theories, respectively (Figure 9c). The mean values of the mechanical properties of the PLA filament obtained from compression tests on cylindrical specimens are presented in Table 2. All the auxetic specimens failed from the struts that were expected to have the maximum local stress in the unit cell (strut DB; Figure 10 shows the specimens from the top view). The maximum differences between the experimental results and volumetric element numerical results for normalized elastic modulus, Poisson's ratio, and normalized yield stress are 22.79%, 16.58%, and 15.59%, respectively. The maximum differences between the experimental results and Timoshenko analytical (by considering the overlapping effect) results for nor- As for the case of analytical relationships that neglect the overlapping effect at the vertices, it is evident that as compared to the analytical results based on the Euler-Bernoulli beam theory, the analytical results based on Timoshenko beam theory have much better agreement with the beam-based FE model results for all relative densities (Figure 9). By increasing the relative density, the difference between the numerical and analytical results for both beam theories increase, and the maximum difference between the numerical and analytical results (at a relative density of µ = 0.4) for normalized elastic modulus are 7.91% and 0.62% for Euler-Bernoulli and Timoshenko beam theories, respectively. As for the Poisson's ratio, the maximum difference between the beam-based FE model and the analytical models based on Euler-Bernoulli and Timoshenko beam theories are 14.84% and 0.91%, respectively. It is worth noting that the maximum negative Poisson's ratio for this case is very high and equal to −0.7791. These numerical/analytical differences for normalized yield stress are 13.76% and 5.32% for Euler-Bernoulli and Timoshenko beam theories, respectively.
The maximum differences between the elastic modulus results of numerical models made from volumetric elements and the analytical models, which consider the overlapping effect, are 54.97% and 45.75% for Euler-Bernoulli and Timoshenko theories, respectively (Figure 9b). These differences for Poisson's ratio are 85.7% and 88.24% for Euler-Bernoulli and Timoshenko beam theories, respectively (Figure 9a). Finally, the noted differences for normalized yield stress are 32.31% and 6.85% for Euler-Bernoulli and Timoshenko beam theories, respectively (Figure 9c).
The mean values of the mechanical properties of the PLA filament obtained from compression tests on cylindrical specimens are presented in Table 2. All the auxetic specimens failed from the struts that were expected to have the maximum local stress in the unit cell (strut DB; Figure 10 shows the specimens from the top view). The maximum differences between the experimental results and volumetric element numerical results for normalized elastic modulus, Poisson's ratio, and normalized yield stress are 22.79%, 16.58%, and 15.59%, respectively. The maximum differences between the experimental results and Timoshenko analytical (by considering the overlapping effect) results for normalized elastic modulus, Poisson's ratio, and normalized yield stress are 28.52%, 25.95%, and 18.14%, respectively.  Figure 11a shows the von Mises stress distribution in the unit cell's struts for relative density of µ = 0.2. From this figure, it can be inferred that the maximum von Mises stress occurs in the same struts expected in analytical and numerical results with beam elements (struts BD). The stress distribution can be better visualized in the sliced unit cell in Figure 11b. The important point to note regarding the numerical results of this FE model is that the maximum stress locations are at the central point of vertex B (and experimental test results). This, again, signifies the importance of considering the effective lengths rather than vertex-to-vertex lengths in the analytical models. The stress and strain distribution for all other unit cells with various relative densities are presented in Sections S.2.4 and S.2.5 of Supplementary Material (Figures S15-S27).
unit cell (strut DB; Figure 10 shows the specimens from the top view). The maximum differences between the experimental results and volumetric element numerical results for normalized elastic modulus, Poisson's ratio, and normalized yield stress are 22.79%, 16.58%, and 15.59%, respectively. The maximum differences between the experimental results and Timoshenko analytical (by considering the overlapping effect) results for normalized elastic modulus, Poisson's ratio, and normalized yield stress are 28.52%, 25.95%, and 18.14%, respectively.   Figure 11a shows the von Mises stress distribution in the unit cell's struts for relative density of = 0.2. From this figure, it can be inferred that the maximum von Mises stress occurs in the same struts expected in analytical and numerical results with beam elements

Effect of Angles and
Angles and are the most critical parameters in the 3D re-entrant structure because these parameters play the main role in creating the negative Poisson's ratio effect ( Figure S9). Therefore, analyzing the effect of these parameters could be very beneficial for predicting the most extreme values of the mechanical properties, in particular, the negative Poisson's ratio. Using 3D surface plot demonstration, Figure 12 presents variations in the mechanical properties of 3D re-entrant unit cell with respect to angles and each varying in the range of −25° to +25°. Moreover, the special cases of = 0°, = 0°, and = are marked in the 3D plots with solid, dashed, and dashed-dot lines, respectively. All the plots are presented for = = = 0.14 and with considering a constant volume for the circumscribed cube of the unit cell. Although one might assume that the mechanical properties of the unit should have symmetry with respect to angle , the results shown in Figure 12 demonstrate that the mechanical properties are not symmetrical with respect to = 0°. Due to the constant volume assumed for circumscribed cube of unit cell, the length changes based on angle , and is shorter for < 0 than for > 0. Therefore, the results are not symmetrical with respect to = 0°, and we have higher levels of normalized elastic modulus and yield stress for < 0 as compared to > 0. According to Figure 12a

Effect of Angles θ and ϕ
Angles θ and ϕ are the most critical parameters in the 3D re-entrant structure because these parameters play the main role in creating the negative Poisson's ratio effect ( Figure S9). Therefore, analyzing the effect of these parameters could be very beneficial for predicting the most extreme values of the mechanical properties, in particular, the negative Poisson's ratio. Using 3D surface plot demonstration, Figure 12 presents variations in the mechanical properties of 3D re-entrant unit cell with respect to angles θ and ϕ each varying in the range of −25 • to +25 • . Moreover, the special cases of θ = 0 • , ϕ = 0 • , and θ = ϕ are marked in the 3D plots with solid, dashed, and dashed-dot lines, respectively. All the plots are presented for r 1 = r 2 = r 3 = 0.14b and with considering a constant volume for the circumscribed cube of the unit cell. Although one might assume that the mechanical properties of the unit should have symmetry with respect to angle θ, the results shown in Figure 12 demonstrate that the mechanical properties are not symmetrical with respect to θ = 0 • . Due to the constant volume assumed for circumscribed cube of unit cell, the length l 1 changes based on angle θ, and l 1 is shorter for θ < 0 than for θ > 0. Therefore, the results are not symmetrical with respect to θ = 0 • , and we have higher levels of normalized elastic modulus and yield stress for θ < 0 as compared to θ > 0. According to Figure 12a,c, the maximum normalized elastic modulus of 0.002234 occurs at θ = −25 • and ϕ = 2 • , and the maximum normalized yield stress of 0.008441 occurs at θ = −25 • and ϕ = 4.5 • . The behavior of Poisson's ratio with respect to variations of angles θ and ϕ is relatively different in comparison with elastic modulus and yield stress. As shown in Figure 12b, the maximum positive and negative values of Poisson's ratio occur at θ = ϕ = −25 • and θ = −ϕ = 25 • , respectively. The graph also demonstrates that by proper combination of θ and ϕ, the structure can be either auxetic or non-auxetic and with the desired extent. Moreover, the graph shows that to have a near-zero Poisson's ratio, the angle θ needs to remain in the range of 4 • < θ < 8.5 • , while angle ϕ can take any value from −25 • to +25 • . To understand better the effect of angles θ and ϕ on the Poisson's ratio of re-entrant structure in other relative densities, the variations of Poisson's ratio for four relative densities (r/b = 0.1, r/b = 0.2, r/b = 0.3, and r/b = 0.4) are presented in Figure 13. According to this figure, the range of Poisson's ratio changes between −0.6 and +0.

Effect of Relative Density on the Mechanical Properties of the Special Case of = = 0°
A critical case can be considered for the structure in which angles and are considered to be zero. The effect of variation in relative density of such structure is plotted in Figure 14. The figure shows that in case of neglecting the overlapping effect, as compared to the analytical results based on Euler-Bernoulli beam theory, the analytical results based on Timoshenko beam theory have better agreement with numerical results of the FE models based on beam elements for all values of relative density. According to Figure 14a, the re-entrant structure does not have a negative Poisson's ratio, and the value of Poisson's ratio for all relative densities is positive and varies from 0.127 to 0.074. Considering the overlapping effect decreases the values of Poisson's ratio, and the minimum value of Poisson's ratio for this case reaches 0.009158. A critical case can be considered for the structure in which angles θ and ϕ are considered to be zero. The effect of variation in relative density of such structure is plotted in Figure 14. The figure shows that in case of neglecting the overlapping effect, as compared to the analytical results based on Euler-Bernoulli beam theory, the analytical results based on Timoshenko beam theory have better agreement with numerical results of the FE models based on beam elements for all values of relative density. According to Figure 14a

Elimination of Strut l 3
According to the contours of mechanical properties of 3D re-entrant unit cell shown in Figure 13, the angle ϕ had negligible effect on the sign of Poisson's ratio in all values of r b . This gives an inspiration to create a new lightened re-entrant unit cell by eliminating strut l 3 from the general 3D re-entrant unit cell. The lightened unit cell with negative, zero, and positive values of θ are shown in Figure 15. The analytical relationships for this unit cell could be easily obtained by setting r 3 = 0 in the general relationships (i.e., Equations (14), (15) and (20)). The resulting mechanical properties of the lightened re-entrant unit cell are presented in Figure 16. The results show that the maximum normalized elastic modulus and yield stress for the lightened unit cell occur at θ = +25 • for all values of relative densities. The maximum negative and positive Poisson's ratio for the lightened unit cell is obtained as −0.6322 and +0.352, respectively. The variations of mechanical properties of the lightened unit cell with respect to relative density in the special condition of θ = 0 have been presented with solid lines in Figure 16. Similar to the general unit cell, the lightened structure has positive Poisson's ratio for all values of relative densities in θ = 0, and the value of Poisson's ratio for this case varies from 0.1688 to 0.09837 (Figure 16b).
Comparison of the general idealized 3D re-entrant unit cell and lightened 3D re-entrant unit cell can provide good insights for designers to decide which design to choose from. Figure 17 shows the mechanical properties contours for general and lightened unit cells in various values of θ and r/b. For quantitative comparison, the comparisons are performed for the mean values of both relative density and r/b, i.e., µ = 0.25 and r/b = 0.25. As for the case of µ = 0.25, the value of normalized elastic modulus in θ = 25 • were calculated as 0.006095 and 0.009178 for the general re-entrant and lightened structures, respectively. As for the case of r/b = 0.25, the values of normalized elastic modulus were obtained as 0.005367 and 0.004184 for the general re-entrant and lightened structures, respectively ( Figure 17 and Figure S10a). Moreover, the normalized yield strength at µ = 0.25 were 0.007666 and 0.01157 for the general re-entrant and lightened structures, respectively ( Figure S11c), while at r/b = 0.25, they were 0.006837 and 0.005808, respectively ( Figure 17 and Figure S10c). Therefore, based on results, we can conclude that the lightened re-entrant structure has better strength and stiffness at the same relative densities in comparison with the general re-entrant structure. In contrast, the lightened re-entrant unit cell is weaker than the general re-entrant structure at the same ratio of r/b (Figure 17). This comparison could be drawn for Poisson's ratio as well. At the same relative density of µ = 0.25 and θ = 25 • , Poisson's ratio of the general and lightened re-entrant structures were obtained as −0.2191 and −0.2673, respectively. These values at the conditions of r/b = 0.25 and θ = 25 • were −0.2318 and −0.3419 for the general and lightened re-entrant structures, respectively ( Figure 17 and Figure S10b). Therefore, eliminating strut l 3 always increases the extent of either the positive or negative values of Poisson's ratio for both the cases of identical relative densities and r/b ratio. modulus and yield stress for the lightened unit cell occur at = +25° for all values of relative densities. The maximum negative and positive Poisson's ratio for the lightened unit cell is obtained as −0.6322 and +0.352, respectively. The variations of mechanical properties of the lightened unit cell with respect to relative density in the special condition of = 0 have been presented with solid lines in Figure 16. Similar to the general unit cell, the lightened structure has positive Poisson's ratio for all values of relative densities in = 0, and the value of Poisson's ratio for this case varies from 0.1688 to 0.09837 ( Figure   16b).  (e) (f) Comparison of the general idealized 3D re-entrant unit cell and lightened 3D re-entrant unit cell can provide good insights for designers to decide which design to choose from. Figure 17

Why Do Timoshenko Relationships Have Better Accuracy?
In this paper, the relationships based on Timoshenko beam theory have been derived for elastic modulus, Poisson's ratio, and yield strength of the re-entrant lattice structure.  Figure S13 presents the effect of effective lengths on mechanical properties of both the general and lightened 3D reentrant structure for a specific value of θ = ϕ = 25 • . According to this figure, it is clear that the influence of effective length in both structures is similar, and considering the effective lengths due to overlapping of the struts at the vertices increases the stiffness and strength of the structure. In addition, the values of Poisson's ratio of both general and lightened structures are decreased due to the higher resistance to deformation by considering the effective lengths. In general, based on Figure S13, it can be concluded that the effective length's influence on the mechanical properties of the lightened structure as compared to its influence on the general re-entrant unit cell are similar. Additional data and figures for comparison between the general and lightened re-entrant structures can be found in Section S.2.3 of the Supplementary Materials (Figures S11-S14).

Why Do Timoshenko Relationships Have Better Accuracy?
In this paper, the relationships based on Timoshenko beam theory have been derived for elastic modulus, Poisson's ratio, and yield strength of the re-entrant lattice structure. Based on the presented results in all cases of studies and also in the all relative densities of structure, the Timoshenko beam theory presents better and more precise results for mechanical properties of the lattice structure in comparison with the Euler-Bernoulli beam theory. The main reason for the accuracy of the Timoshenko beam theory is that this theory takes into account shear deformation and rotational bending effects, making it suitable for describing the behavior of thick beams. For this reason, the results of the analytical model based on Euler-Bernoulli theory deviates from numerical results in higher values of relative densities, while the Timoshenko beam theory results overlap with the numerical results in all values of relative density. It worth noting that taking into account the shear deformation effect lowers the flexural stiffness of the struts and hence the stiffness of the structure, which yields to larger deflections of the struts under a static load. This leads to higher Poisson's ratio and lower elastic modulus values for Timoshenko beam theory in comparison with Euler-Bernoulli beam theory.

The Most Influencing Parameters on the Structure's Poisson's Ratio
In general, angles θ and ϕ are the main design parameters in creating a negative Poisson's ratio effect in the re-entrant structure. According to the results (Figures 12 and 13), the magnitude of θ determines the sign of Poisson's ratio (positive or negative). Increasing θ leads to more negative values of Poisson's ratio in the unit cell. It is important to note that even though the value of ϕ has some effects on the amount of Poisson's ratio, but it does not have any effects on the sign of Poisson's ratio. It is worth mentioning that there exists a specific case (for specific values of θ and ϕ) in the re-entrant structure where the Poisson's ratio of structure becomes zero.
Considering the effective length in the analytical calculations decreases the effects of θ and ϕ, which is a representation of what would be experienced in real-life conditions (as also confirmed by numerical models based on volumetric elements and experimental results). More importantly, the effect of variation in ϕ on Poisson's ratio is almost eliminated when the effective length is taken into consideration (see Section S.2.2 of Supplementary).
According to stress distribution in the re-entrant structure, struts BD are the most influential elements in determining the sign and value of the transverse deformation of the unit cell. Applied load at point A of the unit cell is directly transferred to point B of strut BD via strut AB and this causes strut BD to have normal tensile stress, as shown in Figure 18. However, strut CD is under compressive normal stress due to the load applied to the structure. This makes strut CD demonstrate inverse behavior and counteract regarding the extent of negative Poisson's ratio. Therefore, the effect of angle θ on Poisson's ratio is much greater than that of ϕ, either effective lengths are considered or not. tive densities. Thus, Figures S11−S13 could be a better reference for comparison of the stiffness and strength of the general and lightened re-entrant unit cells. According to these figures, we can conclude that at the same relative densities, the lightened structure could be a better option in comparison with the general structure due to its higher elastic modulus and yield stress. In addition, the results show that the lightened structure could provide higher values of negative or positive Poisson's ratio either in the same relative density or / ratio. All of the above-mentioned comparisons are based on the assumed condition of = in the general 3D re-entrant structure. Therefore, it should be emphasized that the general re-entrant structure can present a much wider range of mechanical properties for ≠ as compared to lightened unit cell due to its additional parameters and . Moreover, with a proper choice of and values, the general structure can present ideal modulated mechanical properties with the desired combination of mechanical properties (Poisson's ratio, elastic modulus, and yield strength). For instance, as it can be seen in Figure 12b, angle has a relatively low effect on Poisson's ratio ( is the governing parameter here), while it has a huge effect on elastic modulus and yield strength. Therefore, as compared to lightened structure, the general re-entrant unit cell has the advantage of giving the opportunity of designing a structure with a desirable combination of elastic modulus and Poisson's ratio.

Some Points Regarding Effective Length (Definition, Effect, Discreparancies, etc.)
All the struts in the analytical solution of the structures were considered as beams, and they were considered to be rigidly connected to each other at the vertices. However, in practice, when the struts reach each other at the joints, they overlap with each other and, in fact, the lengths that contribute to the deformation of struts and the structure as a whole are decreased. Generally, there is no accurate definition for the effective length, and in this study, the minimum distance between the two sides of struts that are not in the overlapping regions has been considered as struts' effective length. This definition has been chosen based on the stress and strain distribution observations in the struts of the volumetric numerical models. In other words, the effective lengths have been defined by the regions of the struts where stress and strain values start to increase abruptly. By considering the effective length, the analytical results give much better agreement with experimental test results which can be very helpful in analyzing the structure in real-life conditions. The other parameter effective on the structure's Poisson's ratio is relative densitythe higher the relative density, the lower the negative Poisson's ratio will be. An ideal model (analytical beam model without consideration of vertex overlapping) has a negative Poisson's ratio in all values of relative densities. However, in more realistic models, and particularly in higher relative densities, (very small) positive Poisson's ratios can be observed (Figure 9). This is due to the positive Poisson's ratio of the bulk material (PLA), which takes more dominance in the structural deformation of the system in higher relative densities.

Differences between General and Lightened 3D Re-Entrant Structures
Based on the normal and bending stress distribution in the struts of the structure shown in Figure 18, it was predictable that eliminating struts CD from the system could contribute to an increase in the value of unit cell Poisson's ratio, even though it would lead to a decrease in the strength and stiffness of the structure. According to comparison results presented in Figure 17, the stiffness and strength of the lightened structures resulted from eliminating struts CD are lower than the general re-entrant structure in the same values of r/b. However, it must be mentioned that the mechanical properties of the lattice structures and metamaterials are usually compared with each other at similar relative densities. Thus, Figures S11−S13 could be a better reference for comparison of the stiffness and strength of the general and lightened re-entrant unit cells. According to these figures, we can conclude that at the same relative densities, the lightened structure could be a better option in comparison with the general structure due to its higher elastic modulus and yield stress. In addition, the results show that the lightened structure could provide higher values of negative or positive Poisson's ratio either in the same relative density or r/b ratio.
All of the above-mentioned comparisons are based on the assumed condition of θ = ϕ in the general 3D re-entrant structure. Therefore, it should be emphasized that the general re-entrant structure can present a much wider range of mechanical properties for θ = ϕ as compared to lightened unit cell due to its additional parameters ϕ and l 3 . Moreover, with a proper choice of θ and ϕ values, the general structure can present ideal modulated mechanical properties with the desired combination of mechanical properties (Poisson's ratio, elastic modulus, and yield strength). For instance, as it can be seen in Figure 12b, angle ϕ has a relatively low effect on Poisson's ratio (θ is the governing parameter here), while it has a huge effect on elastic modulus and yield strength. Therefore, as compared to lightened structure, the general re-entrant unit cell has the advantage of giving the opportunity of designing a structure with a desirable combination of elastic modulus and Poisson's ratio.

Some Points Regarding Effective Length (Definition, Effect, Discreparancies, etc.)
All the struts in the analytical solution of the structures were considered as beams, and they were considered to be rigidly connected to each other at the vertices. However, in practice, when the struts reach each other at the joints, they overlap with each other and, in fact, the lengths that contribute to the deformation of struts and the structure as a whole are decreased. Generally, there is no accurate definition for the effective length, and in this study, the minimum distance between the two sides of struts that are not in the overlapping regions has been considered as struts' effective length. This definition has been chosen based on the stress and strain distribution observations in the struts of the volumetric numerical models. In other words, the effective lengths have been defined by the regions of the struts where stress and strain values start to increase abruptly. By considering the effective length, the analytical results give much better agreement with experimental test results which can be very helpful in analyzing the structure in real-life conditions.

Failure of Idealized 3D Re-Entrant Structure
The idealized 3D re-entrant structure could be considered as a bending dominant lattice structure due to its architecture with several sloped struts. Figure 18 shows the difference between the levels of normal and bending stress components in the struts during a compressive loading condition. According to this figure, it is obvious that the bending stress in both contracted and stretched regions of the struts plays the main role and it has higher magnitudes than normal stress. Therefore, the failure of the structure occurs in the struts with the highest bending stress levels (i.e., strut BD), which is also validated by experimental test results as shown in Figure 10. One would argue that other failure scenarios such as buckling of the struts could come into play, especially in lower values of structure's relative density. The analytical relationship of critical stress for buckling of the idealized 3D re-entrant structure was derived and presented in Section 2.1.4 (Equation (23)). By calculating the values of buckling stress ( Figure S8), it is found that the level of buckling stress for the unit cell is several orders of magnitude higher than yield stress (shown in Figure 9c). Hence, the buckling failure is very unlikely to happen in the idealized 3D re-entrant structure, especially in the range of relative densities applicable for fabricating structures by additive manufacturing.

Limitations and Suggestions
Although the exact analytical solutions based on Euler-Bernoulli and Timoshenko beam theories presented in this study were capable of providing acceptable results for mechanical properties of idealized 3D re-entrant unit cells and their corresponding lattice structure, these relationships are only valid for linear elastic region and small deformations of unit cells. Some of the applications of lattice and cellular metamaterials need large deformations and nonlinear behavior of unit cells to present some extraordinary functionalities [60,61]. For this aim, the material model and the theory for the analytical solution need to be considered as nonlinear elastic or hyper elastic to obtain accurate and realistic results. Hence, the nonlinear solution and analytical relationships, which are valid for nonlinear regions, will be a very valuable step that could be taken for further studies.
For some unit cell types, due to complexity of unit cell shape, asymmetry in the unit cell, asymmetry in the loading condition, size of the structure (being constructed by a large number of cells), and nonlinearity in the deformation (post-yielding behavior), it is not possible to derive analytical solution and/or to perform computational calculations using FE models with micro-structural topology cost-effectively. In such instances, the Cauchy continuum technique [62] could be very beneficial because it considers the behavior of each cell similar to the behavior of continuous materials (without pores). Considering the lattice structures as a continuous material, of course, will bring in some errors in the results because it might not consider some microstructural features. Fortunately, we were able to perform both the analytical and numerical analyses on the more accurate micro-structural models due to symmetry in the unit cell, its deformation being in the linear elastic range, and the possibility of considering the deformation of a unit cell as being representative of the deformation of a lattice structure.

Conclusions
In this paper, a new complex 3D re-entrant structure with negative Poisson's ratio was studied and analytical relationships for its elastic mechanical properties (elastic modulus, yield strength, and Poisson's ratio) were derived. Two different FE models, one based on beam elements and the other based on volumetric elements, were implemented for validating the analytical results. Furthermore, five sets of specimens with different relative density values were additively manufactured and tested under compressive loading conditions, and the obtained experimental data were used for the analysis of the behavior of the structure in actual conditions. The results of the analytical models, which were based on effective length, had good agreement with the volumetric numerical model and experimental results. The effect of various parameters on the mechanical properties of the structure was also studied, and the results demonstrated that the angle θ has the highest effect on the sign of the Poisson's ratio, while the angle ϕ has the lowest effect on the Poisson's ratio. However, the struts corresponding to angle ϕ provide strength and stiffness for the structure. The results also showed that the structure could have zero Poisson's ratio for a specific range of θ and ϕ angles.
In the end, a lightened re-entrant structure was introduced, and its results were compared to those of the general re-entrant structure with θ = ϕ. The results showed that at the same relative densities, the lightened structure could be a better option in comparison with the general structure due to its higher mechanical stiffness and strength. In addition, the results showed that the lightened structure could provide higher values of negative or positive Poisson's ratio. It should be emphasized that the general re-entrant structure can present a much wider range of mechanical properties for θ = ϕ as compared to the lightened unit cell due to its additional geometrical parameters. With a proper choice of θ and ϕ values, the general re-entrant structure can present ideal modulated mechanical properties with the desired combination of mechanical properties (Poisson's ratio, elastic modulus, and yield strength). In summary, the general re-entrant unit cell presented in this research has the advantage of giving the opportunity of designing a structure with a desirable combination of elastic modulus and Poisson's ratio.
The wide range of mechanical properties such as Poisson's ratio and elastic modulus provided by the idealized 3D re-entrant and lightened structures create valuable possibilities to utilize these structures in lattice meta-implants with more compatible stress and strain distributions. In addition, unlike many 3D re-entrant unit cells, the idealized 3D re-entrant unit cell studied here can effectively be used in lattice structures with graded Poisson's ratio distribution due to the fact that no strut is shared with adjacent cells. Hence, the results of this study can be efficiently used in designing and manufacturing metaimplants with the graded distribution of Poisson's ratio to enhance the micromotion at the bone-implant interfaces.
Supplementary Materials: The following are available online at https://www.mdpi.com/1996 -1944/14/4/993/s1, Figure S1: The side view of re-entrant unit cell and its surrounded square, Figure S2: The symmetry planes of 3D re-entrant unit cell, Figure S3: Applied repetitive boundary condition to a single reentrant unit cell, Figure S4: SolidWorks CAD models of additive manufactured specimens, Figure S5: 3D-printed specimens in compression test, Figure S6: Deformations of (a) strut EF, (b) strut CD, (c) strut DG, and (d) strut BD due to q 5 = 1, Figure S7: Variation of relative density with respect to θ and ϕ, Figure S8: Buckling critical stress of general idealized 3D re-entrant unit cell with and without considering effective length, Figure S9: Mechanical properties variations of structure with ϕ = 22.5 • and r/b = 0.14 with respect to θ: (a) elastic modulus, (b) Poisson's ratio, and (c) yield stress, Figure S10: Mechanical properties variations of structure with θ = 22.5 • and r/b = 0.14 with respect to ϕ: (a) elastic modulus, (b) Poisson's ratio, and (c) yield stress., Figure S11: Variation of mechanical properties for general and lightened unit cell with respect to r/b and θ parameters: (a) normalized elastic modulus, (b) Poisson's ratio, and (c) normalized yield stress, Figure S12: Variation of mechanical properties for general and lightened unit cell with respect to µ and θ parameters: (a) normalized elastic modulus, (b) Poisson's ratio, and c) normalized yield stress, Figure S13: Mechanical properties curves of general and lightened 3D re-entrant structures with ϕ = 22.5, θ = 22.5: (a) Poisson's ratio, (b) elastic modulus, and (c) yield stress, Figure S14: Contour plots of mechanical properties for general (left) and lightened (right) unit cell with respect to µ and θ parameters: (a,b) normalized elastic modulus, (c,d) Poisson's ratio, and (e,f) normalized yield stress, Figure S15: Von Mises equivalent stress (Beam 189 elements), Figure S16: Von Mises equivalent strain (Beam 189 elements), Figure S17: 1st principal stress (Beam 189 elements), Figure S18: 2nd principal stress (Beam 189 elements), Figure S19: 3rd principal stress (Beam 189 elements), Figure S20: 1st principal strain (Beam 189 elements), Figure S21: 2nd principal strain (Beam 189 elements), Figure  S22: 3rd principal strain (Beam 189 elements), Figure S23: Normal stress distribution in the beam elements of models made of BEAM 189 element types, Figure S24: Bending stress distribution in the beam elements of models made of BEAM 189 element types, Figure S25: Von Mises equivalent stress for FE models based on solid elements, Figure S26: Von Mises equivalent strain for FE models based on solid elements. Figure S27