Improving the Accuracy of Analytical Relationships for Mechanical Properties of Permeable Metamaterials

: Permeable porous implants must satisfy several physical and biological requirements in order to be promising materials for orthopaedic application: they should have the proper levels of stiffness, permeability, and fatigue resistance approximately matching the corresponding levels in bone tissues. This can be achieved using designer materials, which exhibit exotic properties, commonly known as metamaterials. In recent years, several experimental, numerical, and analytical studies have been carried out on the inﬂuence of unit cell micro-architecture on the mechanical and physical properties of metamaterials. Even though experimental and numerical approaches can study and predict the behaviour of different micro-structures effectively, they lack the ease and quickness provided by analytical relationships in predicting the answer. Although it is well known that Timoshenko beam theory is much more accurate in predicting the deformation of a beam (and as a result lattice structures), many of the already-existing relationships in the literature have been derived based on Euler–Bernoulli beam theory. The question that arises here is whether or not there exists a convenient way to convert the already-existing analytical relationships based on Euler–Bernoulli theory to relationships based on Timoshenko beam theory without the need to rewrite all the derivations from the start point. In this paper, this question is addressed and answered, and a handy and easy-to-use approach is presented. This technique is applied to six unit cell types (body-centred cubic (BCC), hexagonal packing, rhombicuboctahedron, diamond, truncated cube, and truncated octahedron) for which Euler–Bernoulli analytical relationships already exist in the literature while Timoshenko theory-based relationships could not be found. The results of this study demonstrated that converting analytical relationships based on Euler–Bernoulli to equivalent Timoshenko ones can decrease the difference between the analytical and numerical values for one order of magnitude, which is a signiﬁcant improvement in accuracy of the analytical formulas. The methodology presented in this study is not only beneﬁcial for improving the already-existing analytical relationships, but it also facilitates derivation of accurate analytical relationships for other, yet unexplored, unit cell types.


Introduction
Recently, partially or fully porous load-bearing implants have been proposed to replace the traditional solid implants for repairing large bony defects. While metallic foams manufactured by conventional techniques such as powder metallurgy [1], investment casting [2], and space-holder [3,4] have found their way in this field [5], they all lack a good controllability over the microstructural geometry of the implants, and hence their static mechanical properties [6], fatigue resistance [7], and biological response [8]. A recent

Euler-Bernoulli Beam Theory
For a beam under no distributed load, the Euler-Bernoulli beam equation can be written as where is the deflection. The solution to this differential equation can be expressed as: where constants, − , are determined by applying boundary conditions. For a cantilever Euler-Bernoulli beam with a point load, , at its end, we have On the other hand, for a cantilever beam with a concentrated moment, , at its end, the displacement and rotation are as follows: In beams where the angle of the free end does not change during the deformation (e.g., the beams of the lattice structure considered in this study), the rotations produced by the lateral load, , and moment, , must be equal and opposite, from which the value of can be identified: While force, , tends to increase the deflection, moment, , tends to reduce it. The total deflection created by force, , and moment, , is then

Euler-Bernoulli Beam Theory
For a beam under no distributed load, the Euler-Bernoulli beam equation can be written as where w is the deflection. The solution to this differential equation can be expressed as: where constants, c 0 − c 3 , are determined by applying boundary conditions. For a cantilever Euler-Bernoulli beam with a point load, P, at its end, we have δ = Fl 3 3E s I and θ = Fl 2 2E s I On the other hand, for a cantilever beam with a concentrated moment, M, at its end, the displacement and rotation are as follows: In beams where the angle of the free end does not change during the deformation (e.g., the beams of the lattice structure considered in this study), the rotations produced by the lateral load, F, and moment, M, must be equal and opposite, from which the value of M can be identified: While force, F, tends to increase the deflection, moment, M, tends to reduce it. The total deflection created by force, F, and moment, M, is then Rewriting Equation (6) as a function of F gives (see Figure 2a) = 3 − 2 2 = 12 (6) Rewriting Equation (6) as a function of gives (see Figure 2a) = 12 (7) According to Equation (5), we will have = . Similarly, the axial force required to displace the end of a rod for is / (see Figure 2c). The equations for a cantilever beam with rotation but with no displacement in the end can be obtained in a similar manner and the results are shown in Figure 2b.

Timoshenko Beam Theory
Now, we try to find out how the moments and forces shown in Figure 2a,b would change due to change in the beam theory. The Timoshenko beam theory takes into account shear deformation and rotational inertia effects, making it suitable for describing the behaviour of short beams. For a homogenous beam of constant cross-section, the Timoshenko beam governing equations are as follows: where is the angle of rotation of the normal to the mid-surface of the beam and is the shear coefficient factor. The coefficient, , is a dimensionless quantity, dependent on the shape of the cross-section, which is introduced to account for the fact that the shear stress and strain are distributed not-uniformly over the cross-section [38]. In a linear elastic Timoshenko beam, the bending moment, , and the shear force, , are related to the angle of rotation, , and the deflection, , by (a) A cantilever beam with lateral displacement but with no rotation in the end (Figure 3a): Since the distributed load (force per length), ( ), is zero, the solution to the first differential equation of Equation (8) can be expressed as: By applying the boundary condition of no rotation at the root of the cantilever beam, Equation The two other constants ( , ) can be found from the boundary conditions in = 0 and = . At the root ( = 0), the deflection equals zero, hence: Forces and moments required to cause (a) lateral displacement, δ,, with no rotation at the free end of the beam, (b) rotation, θ, with no lateral displacement at the free end of the beam, (c) pure axial extension, and (d) pure twist in the free end of an Euler-Bernoulli beam [32,36].
According to Equation (5), we will have M = 6E s I l 2 δ. Similarly, the axial force required to displace the end of a rod for u is AE s u/l (see Figure 2c). The equations for a cantilever beam with rotation but with no displacement in the end can be obtained in a similar manner and the results are shown in Figure 2b.

Timoshenko Beam Theory
Now, we try to find out how the moments and forces shown in Figure 2a,b would change due to change in the beam theory. The Timoshenko beam theory takes into account shear deformation and rotational inertia effects, making it suitable for describing the behaviour of short beams. For a homogenous beam of constant cross-section, the Timoshenko beam governing equations are as follows: where ϕ is the angle of rotation of the normal to the mid-surface of the beam and κ is the shear coefficient factor. The coefficient, κ, is a dimensionless quantity, dependent on the shape of the cross-section, which is introduced to account for the fact that the shear stress and strain are distributed not-uniformly over the cross-section [38]. In a linear elastic Timoshenko beam, the bending moment, M xx , and the shear force, Q x , are related to the angle of rotation, ϕ, and the deflection, w, by (a) A cantilever beam with lateral displacement but with no rotation in the end (Figure 3a): Since the distributed load (force per length), q(x), is zero, the solution to the first differential equation of Equation (8) can be expressed as:

From Euler-Bernoulli to Timoshenko
According to the relationships obtained for Euler-Bernoulli and Timoshenko beam theories in Sections 2.1 and 2.2, the analytical relationships for elastic modulus and Poisson's ratio of any structure based on Euler-Bernoulli theory can be converted into rela- By applying the boundary condition of no rotation at the root of the cantilever beam, Equation (10) gives C 2 = 0. Similarly, applying the boundary condition of no rotation (ϕ = 0) at the end of the beam (x = l) to Equation (10) gives C 1 = − C 0 l 2 . Substituting C 1 and C 2 in the second line of Timoshenko beam theory (Equation (8)) gives the deflection function of the beam as The two other constants (C 0 , C 3 ) can be found from the boundary conditions in x = 0 and x = l. At the root (x = 0), the deflection equals zero, hence: Moreover, the deflection at the end of the beam equals δ, which gives: By substituting C 0 , C 1 , C 2 , C 3 in the first line of Equation (8) and in Equation (11), the Timoshenko beam theory for a cantilever beam with lateral displacement but with no rotation could be obtained as follows: Now, it is possible to find the moments and forces shown in Figure 3a. As mentioned above, the bending moment, M xx , and the shear force, Q x , are related to the angle of rotation, ϕ, and the deflection, w, by Equation (9). Therefore, by substituting M = M 0 and Q x = F at x = l in respectively the first and second lines of Equation (9), we have: (b) A cantilever beam with rotation but with no lateral displacement in the end ( Figure 3b): Since the distributed load, q(x), and the boundary condition of the beam at the clamped side (x = 0) for this case is similar to the previous case (Case a), the relationship for angle of rotation, ϕ, is the same as that given in Equation (10), and also C 2 = 0. By considering w = 0 at both x = 0 and x = l, the constants C 1 and C 3 can be found as: The beam at the free side (x = l) has a rotation with the known value of ϕ = θ, which after being substituted in Equation (10) gives: Therefore, the Timoshenko beam theory governing equations for a cantilever beam with rotation but with no lateral displacement at the end could be obtained as follows: By considering M xx = M 0 and Q x = F at x = l in Equation (9), the relationship of bending moment and force can be found as

From Euler-Bernoulli to Timoshenko
According to the relationships obtained for Euler-Bernoulli and Timoshenko beam theories in Sections 2.1 and 2.2, the analytical relationships for elastic modulus and Poisson's ratio of any structure based on Euler-Bernoulli theory can be converted into relationships based on Timoshenko beam theory by making the replacements suggested by Table 1 in the stiffness matrix or in the derivation formulas. To evaluate the effectiveness of the conversion approach presented in Table 1, six well known strut-based lattice structures including BCC, hexagonal packing, rhombicuboctahedron, diamond, truncated cube, and truncated octahedron (Figures 4 and 5) were considered. In the following, the procedure of converting Euler-Bernoulli formulas [26,[30][31][32]35,39] into equivalent Timoshenko ones is described for each of the six unit cell types.

From Euler-Bernoulli to Timoshenko
According to the relationships obtained for Euler-Bernoulli and Timoshenko beam theories in Sections 2.1 and 2.2, the analytical relationships for elastic modulus and Poisson's ratio of any structure based on Euler-Bernoulli theory can be converted into relationships based on Timoshenko beam theory by making the replacements suggested by Table 1 in the stiffness matrix or in the derivation formulas. To evaluate the effectiveness of the conversion approach presented in Table 1, six well known strut-based lattice structures including BCC, hexagonal packing, rhombicuboctahedron, diamond, truncated cube, and truncated octahedron (Figures 4 and 5) were considered. In the following, the procedure of converting Euler-Bernoulli formulas [26,[30][31][32]35,39] into equivalent Timoshenko ones is described for each of the six unit cell types.         The commonly known geometries for the six noted unit cell types are presented in Figure 4, and the actual geometries used for both the analytical and numerical analyses in this study are illustrated in Figure 5. The reason of altering the unit cell shape for the BCC case is obvious as the molecular structure of BCC ( Figure 4a) is composed of spheres rather than struts ( Figure 5a). In the case of diamond, the unit cell shown in Figure 4b is rotated, and the unit cell in a specific direction (for which all the struts have similar angles with respect to the horizontal plane) has been considered for analysis ( Figure 5b). In the other cases, the unit cell position has been shifted in order to avoid neighbouring cells having adjacent struts. If the neighbour cells have adjacent side-by-side edges, the analytical relationship obtained for the unit cell does not represent that of a lattice structure. More explanations regarding this can be found in Section 2.1.2 of [32]. Although in this paper, all the analytical relationships for the lattice structures and unit cells are presented in a normalized manner and, hence, the dimensions of the unit cells do not affect the results, an equal volume of 5 × 5 × 5 mm 3 was considered for all the unit cell types. Moreover, in each unit cell, the strut radius was increased from very small values (r ∼ = 0) to high values up to the point where relative density of the unit cell reached µ = 0.5.
The conversion procedure from Euler-Bernoulli beam theory into Timoshenko beam theory for each unit cell is described in the following: (a). BCC The relationships for elastic modulus and Poisson's ratio of BCC unit cell have been presented in Equations (A18) and (A19) of Appendix B (extracted from [31]). In these equations, the axial extension ( AE s l ) and lateral bending ( 12E s I l 3 ) terms of the Euler-Bernoulli theory can be identified easily. By substituting the term 12E s I l 3 with 1 l 3 12Es I + l κ AGs , the relationships for elastic modulus and Poisson's ratio can be converted.

(b). Diamond
For transforming relationships of elastic modulus and Poisson's ratio from Euler-Bernoulli into Timoshenko beam theory, the basic Equations (A20) and (A21) of Appendix B (extracted from [26]) for elastic modulus and (A22) and (A23) of Appendix B (extracted from [26]) for Poisson's ratio were considered. Axial extension ( AE s l ) and lateral bending ( 12E s I l 3 ) terms of Euler-Bernoulli theory can be identified easily. By substituting the term (c). Hexagonal packing For this unit cell, since the analytical relationships for mechanical properties of structure have not been derived based on Euler-Bernoulli theory in the literature [35], the analytical relationships for Euler-Bernoulli and Timoshenko beam theories have both been derived in this study. The detailed derivations can be found in Appendix C.

(d). Rhombicuboctahedron
Since the final relationships for elastic modulus and Poisson's ratio of this unit cell have been presented in r l terms, and the stiffness matrix of unit cell contains AE s l , 12E s I l 3 , G s J l , and 6E s I l 2 , recognition of terms provided in Table 1 in the final Equations (A24) and (A25) of Appendix B (extracted from [30]) is not possible. Therefore, by substituting the term 12E s I l 3 with  [30]) and solving the system of equations, the relationships based on Timoshenko theory were obtained.
(e). Truncated cube The procedure for truncated cube unit cell is very similar to what was described above for the rhombicuboctahedron unit cell. The final stiffness matrix of the unit cell in Equation (A27) of Appendix B (extracted from [32]) contains AE s l and 12E s I l 3 terms. Therefore, by substituting the term 12E s I l 3 with , the stiffness matrix based on Timoshenko beam theory can be derived. Afterwards, by solving the system of equations, the mechanical properties relationships based on Timoshenko beam theory were obtained.

(f). Truncated octahedron
The relationships for elastic modulus and Poisson's ratio of truncated octahedron unit cell have been presented in Equation (A28) and (A29) of Appendix B (extracted from [39]). Since the deformation of this unit cell includes axial extension and lateral bending, the terms AE s l and 12E s I l 3 can be extracted from these equations. By substituting the term 12EI , the relationships presented for elastic modulus and Poisson's ratio can be transformed into corresponding Timoshenko ones. It is worth noting that the relationships for normalized yield stress for five of the six noted unit cells (BCC, hexagonal packing, rhombicuboctahedron, diamond, and truncated octahedron) based on Euler-Bernoulli and Timoshenko beam theories have been obtained separately in another study [40].

Numerical Analysis
The 3D representation of the unit cells used for numerical modelling and analysis are demonstrated in Figure 5. The actual FE models of the unit cells along with the boundary conditions are demonstrated in Appendix D. The struts of the unit cells were discretized using Timoshenko beam elements (element type BEAM189 in ANSYS), and each strut was discretized using five beam elements. The beam elements were rigidly connected to each other at their shared vertices, and they were not allowed to rotate in any direction at the connecting point. The mechanical properties of the titanium alloy Ti-6Al-4V-ELI were used for modelling the behaviour of the matrix material in the FE models. A linear elastic material model with elastic modulus and Poisson's ratio of 113.8 GPa and 0.342, respectively, was implemented. Since the BEAM189 element uses linear interpolation and takes transverse shear deformation into account, it is expected that the numerical results will be closer to the Timoshenko analytical solution. In the FE models of BCC, diamond, truncated cube, and truncated octahedron, a single unit cell with periodic boundary condition was analysed under compressive loading ( Figure A2a,b,e,f in Appendix D). For the hexagonal packing and rhombicuboctahedron topologies, lattice structures consisting of 11 × 11 × 11 unit cells were used for numerical modelling. This was due to the complexity of modelling the repetitive boundary conditions in these two unit cell types. More specifically, in these two cases, each side of the unit cell is composed of two strut types (rather than one strut type in the case of other unit cell types), which have non-symmetrical displacements under compressive loading (demonstrated as strut types A and B in Figure 5c,d).
In all the FE models, the lowermost nodes of the structure were fixed in the direction parallel to the loading direction and were not allowed to rotate in any direction. For the case of FE models made out of single unit cells ( Figure A2a,b,e,f in Appendix D), the side vertices of the unit cell were constrained rotationally as they were symmetrically connected to the (imaginary) adjacent unit cells (repetitive boundary condition). In all the FE models, a downward displacement was applied on the uppermost node(s) of the structure to induce axial deformation. Moreover, the uppermost nodes were not allowed to rotate in any directions.
Mechanical properties of the FE models have been calculated based on the basic definition of elastic modulus, Poisson's ratio, and yield stress: was used for calculating numerical elastic modulus, where L UC is the structure length in the direction parallel to loading direction, A UC is the cross-sectional area of the structure in the direction perpendicular to the loading direction, δ UC is the downward displacement applied to the uppermost nodes, and F UC is obtained by summing the reaction forces of the lowermost nodes.
was used for obtaining Poisson's ratio. In this formula, δ 1 and L 1 are the downward displacement applied to the uppermost nodes and unit cell's length in the direction parallel to loading direction, respectively. Parameters δ 2 and L 2 are respectively the lateral displacement of the side nodes and the structure length in the direction perpendicular to loading direction. • Yield stress: The formula σ y σ ys = F UC A UC σ max was used to calculate normalized yield stress. In this formula, σ max is the maximum von Mises stress experienced in the most critical point of the structure. The critical points of each unit cell can be seen in Section 4.1.
In the cases where a lattice structure was implemented for numerical modelling (rhombicuboctahedron and hexagonal packing), all the above-mentioned terms denoted by UC should rather be denoted by Lattice, and the lattice structure dimensions should be used for calculations.

Results
According to the initial results, by not considering the shear deformation effect in the beam theory, the forces and moments required to create a particular deformation in a single strut could be predicted by 15-20% higher for r/l as large as 0.15. The complete results are presented in Appendix E.
The transformed elastic modulus, Poisson's ratio, and yields stress relationships for the six geometries are presented in Tables 2-4 , and W = G s J l and have been used in the noted tables for the relationships based on the Timoshenko beam theory. For four of the geometries (BCC [31], rhombicuboctahedron [30], truncated cube [32], and truncated octahedron [39]), the original Euler-Bernoulli relationships have been presented in Tables 2 and 3 as well. In the case of hexagonal packing, the original Euler-Bernoulli relationships [35] for elastic modulus and Poisson's ratio have been improved and adjusted (the description on how and why can be found in Appendix C), and the improved relationships are presented in Tables 2 and 3. For the case of diamond unit cell, the original relationships for elastic modulus and Poisson's ratio have been conserved but stated as a function of r/l rather than of relative density, µ. As for the yield stress (Table 4), the analytical Euler-Bernoulli relationships presented in our other paper [40] have been presented.
In order to compare the already-existing [26,[30][31][32]35,39] or improved Euler-Bernoulli analytical relationships with the newly transformed Timoshenko analytical relationships, the results of the analytical relationships for both the Euler-Bernoulli and Timoshenko beam theory have been compared with their numerical and experimental counterparts in Figures 6-8. The relative density relationships for the six geometries are presented in Table 5.
For all the geometries, the newly transformed Timoshenko relationships have exceptionally good agreement with the numerical results for all the mechanical properties: elastic modulus, Poisson's ratio, and yield stress (Figures 6-8). As for the previously obtained formulas obtained in the literature or newly adjusted Euler-Bernoulli formulas, the maximum difference between Euler-Bernoulli analytical elastic modulus and corresponding numerical values for BCC, diamond, hexagonal packing, rhombicuboctahedron, truncated cube, and truncated octahedron (at relative density of µ = 0.5) are, respectively, 21.34%, 57.71%, 20.21%, 14.52%, 14.98%, and 45.54%. However, the corresponding differences between the Timoshenko analytical relationships and the numerical values for the noted geometries are, respectively, 1.13%, 2.21%, 8.29%, 2.97%, 0.43%, and 3.15%.
As for the yield stress (Figure 8), the analytical Euler-Bernoulli and Timoshenko relationships are identical for BCC, diamond, and truncated cube unit cells. The reason is explained in Section 4.3. However, for three other geometries, namely hexagonal packing,   [39] Similarly, the maximum differences between Euler-Bernoulli analytical Poisson's ratio and corresponding numerical values for BCC, diamond, hexagonal packing, rhombicuboctahedron, truncated cube, and truncated octahedron are, respectively, 3.07%, 27.69%, 28.00%, 21.54%, 73.58%, and 40.51%. The corresponding differences between the Timoshenko analytical relationships and the numerical values for the noted geometries are, respectively, 0.133%, 0.826%, 0.899%, 6.24%, 1.36%, and 0.498%, which shows a significant improvement in the accuracy of the analytical relationships. As can be seen in Figures 6-8, the numerical/analytical discrepancy for the case of Timoshenko beam theory is significantly less (around 1/10) than that of Euler-Bernoulli beam theory.
As for the yield stress (Figure 8), the analytical Euler-Bernoulli and Timoshenko relationships are identical for BCC, diamond, and truncated cube unit cells. The reason is explained in Section 4.3. However, for three other geometries, namely hexagonal packing, rhombicuboctahedron, and truncated octahedron, the analytical relationships based on Euler-Bernoulli differ from the relationships obtained based on Timoshenko beam theory. Nevertheless, for all the cases, the Timoshenko analytical yield stress curve has exceptionally good agreement with numerical results (Figure 8). For BCC, diamond, and truncated cube unit cells (the geometries that have the same yield stress analytical relationships for Euler-Bernoulli and Timoshenko beam theories), the maximum difference between analytical and numerical values is less than 0.25%. As for the three other geometries, the maximum difference between Euler-Bernoulli analytical normalized yield stress and corresponding numerical values for hexagonal packing, rhombicuboctahedron and truncated octahedron are, respectively, 21.25%, 3.02%, and 7.52%. However, such differences for the analytical relationships based on Timoshenko beam theory and numerical values are, respectively, 9.49%, 0.09%, and 0.162%, which shows a significant improvement in the accuracy of the analytical relationships.
As for the proximity of the analytical/numerical elastic modulus and yield stress values to the experimental data, the experimental data points have been presented in Figures 6 and 8 for rhombicuboctahedron, truncated cube, and diamond unit cells only, as in the literature, there are experimental measurements for such unit cells only [41]. As for the elastic modulus, converting the analytical relationships from Euler-Bernoulli to Timoshenko has led to closer proximity of analytical and experimental values (Figure 6b,d,e). As for the yield stress, for the three geometries for which experimental data points are available (diamond, rhombicuboctahedron, and truncated cube), the analytical relationships based on Euler-Bernoulli and Timoshenko are identical or almost overlapping, and both are in good agreement with experimental data points (Figure 8b,d,e).

Unit Cell's Behaviour
The main mechanism of deformation in the struts of lattice structures and porous materials is flexure, stretching/contraction, or a combination of them. In the unit cells with a high fraction of vertical struts (struts aligned with the loading condition), the main reason for collapse is the axial normal stresses generated in their struts. Therefore, the deformation in these structures is stretching-dominated and the lattice structure collapses due to the generation of unbearable axial stress in the struts. As can be seen in Figure 9, the hexagonal packing and truncated cube unit cells could be considered as stretching-dominated structures due to the high presence of vertical struts in their architecture, while the other structures namely BCC, diamond, truncated octahedron as well as rhombicuboctahedron could be considered as bending-dominated structure, as they are mostly made up of oblique struts. The main mode of deformation in a unit cell (stretching-dominated or bendingdominated) determines the general deformation of the unit cell and, therefore, its stiffness and yield strength. According to Figure 9, the critical points of the hexagonal packing and truncated cube unit cells are located in the vertical struts due to their stretching-dominated behaviour. On the other hand, the critical points of BCC, diamond, truncated octahedron, and rhombicuboctahedron unit cells are located at the end of oblique struts due to their bending-dominated behaviour. In addition, it can be seen in Figures 6 and 8 that among all the unit cell types, the hexagonal packing and truncated cube structures have the highest stiffness and yield strength, especially in lower values of relative density due to their stretching-dominated behaviour. It is worth noting that since the rhombicuboctahedron has vertical struts, it has an in-between behaviour and gives higher stiffness and strength in comparison with other bending-dominated unit cells. More figures that illustrate the axial and bending stresses in the struts of the unit cells can be found in Appendix D.

Unit Cell's Behaviour
The main mechanism of deformation in the struts of lattice structures and porous materials is flexure, stretching/contraction, or a combination of them. In the unit cells with a high fraction of vertical struts (struts aligned with the loading condition), the main reason for collapse is the axial normal stresses generated in their struts. Therefore, the deformation in these structures is stretching-dominated and the lattice structure collapses due to the generation of unbearable axial stress in the struts. As can be seen in Figure 9, the hexagonal packing and truncated cube unit cells could be considered as stretching-dominated structures due to the high presence of vertical struts in their architecture, while the other structures namely BCC, diamond, truncated octahedron as well as rhombicuboctahedron could be considered as bending-dominated structure, as they are mostly made up of oblique struts. The main mode of deformation in a unit cell (stretching-dominated or bending-dominated) determines the general deformation of the unit cell and, therefore, its stiffness and yield strength. According to Figure 9, the critical points of the hexagonal packing and truncated cube unit cells are located in the vertical struts due to their stretching-dominated behaviour. On the other hand, the critical points of BCC, diamond, truncated octahedron, and rhombicuboctahedron unit cells are located at the end of oblique struts due to their bending-dominated behaviour. In addition, it can be seen in Figures 6  and 8 that among all the unit cell types, the hexagonal packing and truncated cube structures have the highest stiffness and yield strength, especially in lower values of relative density due to their stretching-dominated behaviour. It is worth noting that since the rhombicuboctahedron has vertical struts, it has an in-between behaviour and gives higher stiffness and strength in comparison with other bending-dominated unit cells. More figures that illustrate the axial and bending stresses in the struts of the unit cells can be found in Appendix D.

Why the New Relationships Give Much Better Accuracy?
In this paper, new relationships based on Timoshenko beam theory have been derived for elastic modulus, Poisson's ratio, and yield strength of several topologies. In addition to implementing Timoshenko beam theory, some adjustments in deriving analytical relationships based on Euler-Bernoulli beam theory have been implemented for hexagonal packing unit cell (see Appendix C). Timoshenko beam theory takes into account shear deformation and rotational bending effects, making it suitable for describing the behaviour of thick beams. As a result, the new relationships based on Timoshenko beam theory give much better accuracy even at high relative densities. This improvement is significant for analytical/experimental agreement and exceptional for analytical/numerical agreement. To give a more physically tangible understanding, it must be noted that taking into account the shear deformation effect increases the flexibility of the beam, which effectively leads to larger deflections of the struts (and therefore the whole lattice structure) under an imposed load. This leads to respectively lower and higher elastic modulus and Poisson's ratio of the structure for Timoshenko theory as compared to Euler-Bernoulli theory. It is worth noting that in the bending-dominated unit cells, the discrepancy between the Euler-Bernoulli and Timoshenko results is much greater as compared to that for stretch-dominated unit cells ( Figure 6). In other words, the inability of the Euler-Bernoulli beam theory to predict the effective elastic moduli of open-cell structures based on BCC, diamond, and truncated octahedron unit cells can be attributed to their main mode of deformation (bending-domination) and the importance of considering the shear deformation effect.

Yield Strength
As mentioned in the Results section, the normalized yield stress relationships based on Timoshenko and Euler-Bernoulli theories are identical for BCC, diamond, and truncated cube structures and both theories have good agreement with the numerical results. However, for the hexagonal packing, truncated octahedron, and rhombicuboctahedron structures, the Euler-Bernoulli and Timoshenko yield strength results are quite different (particularly for the hexagonal packing case), and the Timoshenko analytical results show much better overlapping with numerical results as compared to that for the Euler-Bernoulli analytical curve. The reason why the Timoshenko and Euler-Bernoulli yield strengths are identical for some geometries and different for some other is described extensively in [40], but it is explained here briefly. As introduced above, = and =

Why the New Relationships Give Much Better Accuracy?
In this paper, new relationships based on Timoshenko beam theory have been derived for elastic modulus, Poisson's ratio, and yield strength of several topologies. In addition to implementing Timoshenko beam theory, some adjustments in deriving analytical relationships based on Euler-Bernoulli beam theory have been implemented for hexagonal packing unit cell (see Appendix C). Timoshenko beam theory takes into account shear deformation and rotational bending effects, making it suitable for describing the behaviour of thick beams. As a result, the new relationships based on Timoshenko beam theory give much better accuracy even at high relative densities. This improvement is significant for analytical/experimental agreement and exceptional for analytical/numerical agreement. To give a more physically tangible understanding, it must be noted that taking into account the shear deformation effect increases the flexibility of the beam, which effectively leads to larger deflections of the struts (and therefore the whole lattice structure) under an imposed load. This leads to respectively lower and higher elastic modulus and Poisson's ratio of the structure for Timoshenko theory as compared to Euler-Bernoulli theory. It is worth noting that in the bending-dominated unit cells, the discrepancy between the Euler-Bernoulli and Timoshenko results is much greater as compared to that for stretch-dominated unit cells ( Figure 6). In other words, the inability of the Euler-Bernoulli beam theory to predict the effective elastic moduli of open-cell structures based on BCC, diamond, and truncated octahedron unit cells can be attributed to their main mode of deformation (bending-domination) and the importance of considering the shear deformation effect.

Yield Strength
As mentioned in the Results section, the normalized yield stress relationships based on Timoshenko and Euler-Bernoulli theories are identical for BCC, diamond, and truncated cube structures and both theories have good agreement with the numerical results. However, for the hexagonal packing, truncated octahedron, and rhombicuboctahedron structures, the Euler-Bernoulli and Timoshenko yield strength results are quite different (particularly for the hexagonal packing case), and the Timoshenko analytical results show much better overlapping with numerical results as compared to that for the Euler-Bernoulli analytical curve. The reason why the Timoshenko and Euler-Bernoulli yield strengths are identical for some geometries and different for some other is described extensively in [40], but it is explained here briefly. As introduced above, T = 12E s I l 3 and V = 6E s I l 2 for Euler-Bernoulli beam theory, and T = the theories give identical relationships in these two unit cells. As for the truncated cube case, the normalized yield stress in the critical strut could be obtained from the equilibrium of forces and moments, and the failure is caused merely due to axial stress. Therefore, the Euler-Bernoulli and Timoshenko theories give the same results again. Nevertheless, for the three other unit cells, first, the displacement caused by bending needs to be taken into account in the calculations, and second, the formulas in the final stages of derivation are not a mere function of V T , and the V and T terms are present independently and not in a fraction form in such a way that one is a factor of the other. Therefore, the results differ for Euler-Bernoulli and Timoshenko theories for the three other unit cells.

Some Points Regarding Experimental Data Points
Both analytical and numerical techniques over-predict the experimental elastic modulus of the diamond and truncated cube structures, while they under-predict their experimental yield stress. The random irregularities and imperfections created during the AM process diminish the elastic modulus of the lattice structures. These defects create weak links in the structure that lower the mechanical properties of structures. However, in the case of yield stress, when the initial regions of the lattice structures in the critical points (which experience the highest stress levels in the whole lattice structure) are yielded, their local yield stress increases due to strain-hardening phenomenon, and the structure is still able to remain almost intact as the strain-hardening strengthens the structure at the initial damage points. However, when the external load is increased to higher values, the second (and possibly the third) set of critical points are damaged. After the plasticity of the second type (and possibly third type) of failure points following the failure of initial points, the structure is unable to keep its integrity as it was before, as now damage and softening has propagated throughout the whole structure. That is why, in practice, the yielding in the structure usually occurs under external load levels, which are between the external load levels that theoretically cause the first and second set of critical points become locally yielded. This was shown in [36].

Application to Biomedical Implants
There are several applications for cellular materials including heat exchangers, filters, load bearing components, and biomedical implants. Lattice structures can improve the implants' performance significantly, from both mechanical and biological points of view. Obtaining the accurate characteristics and mechanical properties of the lattice structures is necessary to facilitate their use in orthopaedic implants. The mechanical properties of lattice structures depend mainly on the following three parameters: the material they are made of, the cell topology, and the relative density. Obtaining accurate analytical relationships for different cell topologies is a crucial factor in developing computational tools for designing patient-specific implants with non-uniform mechanical property distribution. In this paper, we presented a new and convenient method to transform the analytical relationships of a lattice structure from Euler-Bernoulli theory into Timoshenko theory. Furthermore, using the transformation relationships presented in this paper (Table 1), developing new analytical relationships based on Timoshenko beam theory becomes easier than how it has been before.
For implant design, it is preferable to have a high yield strength and yet access to a wide range of stiffness. Moreover, having similar mechanical properties in three main directions is another factor that should be considered when choosing a unit cell for an implant. This will avoid unwanted deformations generated in the implant when the implant is placed inside the complex geometrical void it is designed for. The results of this study (Figures 6 and 8) show that among the six unit cells considered, truncated cube followed by diamond and rhombicuboctahedron can best satisfy the above-mentioned characteristics.

Limitations of the Present Approach
Although the conversion methodology and the analytical solutions based on Timoshenko beam theory presented in this study could give very accurate results for elastic mechanical properties of open-cell unit cells and lattice structure, they are still based on linear elastic deformation of structures. Therefore, the analytical relationships are valid only for small deformations (i.e. in the elastic range). The analytical relationships presented in this study could be beneficial for several applications where the deformations remain in the elastic regime, such as in biomedical implants. In some applications of lattice structures and porous materials such as energy absorption and actuation, the matrix material undergoes plastic or hyperelastic deformation, and hence the lattice structures behave non-linearly. Therefore, the range of strains and whether or not the lattice structure's behaviour is linear is one of the most important factors that should be considered before utilization of the analytical relationships presented in this study.

Conclusions
In this paper, a new methodology to conveniently convert analytical relationships based on Euler-Bernoulli to equivalent Timoshenko ones was presented. Six unit cells for which Euler-Bernoulli analytical relationships could be found in the literature, but Timoshenko theories could not be found were considered: BCC, hexagonal packing, rhombicuboctahedron, diamond, truncated cube, and truncated octahedron [26,[30][31][32]35,39].
The results of this study demonstrated that converting analytical relationships based on Euler-Bernoulli to equivalent Timoshenko ones can decrease the discrepancy between the analytical and numerical values by one order of magnitude, which is a significant improvement in the accuracy of the analytical formulas. The highest improvement in the analytical relationships was observed in the bending-dominated structures such as BCC, diamond, and truncated octahedron topologies and especially at higher relative densities, where the shear deformation effect becomes significant. While the conversion methodology introduced in this study had a significant effect on improvement of elastic modulus and Poisson's ratio relationships, its effect on yield stress relationships was much less, and in some cases (such as BCC, diamond, and truncated cube), it was negligible. The methodology presented in this study is not only beneficial for improving the already-existing analytical relationships, but it also facilitates the derivation of analytical relationships for other, yet unexplored, unit cell types.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to their large size and unavailability of storage.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Derivation of Forces and Moments for a Single Strut
Appendix A. 1

. Cantilever Beam with Displacement without Rotation in the End
The governing equations of Timoshenko beam theory are: Sequential integration from the first line of Equation (A1) gives: And integration from the second line of Equation (A1) gives: There are four boundary conditions at the root and at the end of the beam: which in turn gives By differentiation from Equation (A9), the point load and at the end of the cantilever beam can be obtained:

A.2. Cantilever Beam with Rotation without Displacement in the End
The beam governing equations are the same (Equation (A1-A3)). Applying the relevant boundary condition for this beam:

Appendix C. New Analytical Relationships for Hexagonal Packing Geometry
In this appendix, the derivation of new analytical relationships for the hexagonal packing geometry is presented.

Appendix C. New Analytical Relationships for Hexagonal Packing Geometry
In this appendix, the derivation of new analytical relationships for the hexagonal packing geometry is presented. Hexagonal closed-packing (HCP) is known as the mostly efficient way a space can be filled by an arrangement of spheres ( Figure A1a). Sphere constructing an HCP structure fill up 74% of space. For an open-cell lattice structure from the concept of HCP, the spheres can be visualized to inflate from all direction until they create a second lattice structure made of polyhedrons (demonstrated by black solid lines in Figure A1a). To find the unit cell of the new polyhedral, first we select two rows of spheres on top of each other and the corresponding polyhedral lattice structure surrounding them ( Figure A1b). After removal of the spheres we are left with six polyhedral cells, in which we can find a unit cell which after tessellation in space can create the HCP-based polyhedral lattice structure (the red unit cell highlighted in Figure A1c). Dimensions and loads acting on a hexagonal packing unit cell are shown in Figure A1d. The initial spheres demonstrated in Figure A1a  Hexagonal closed-packing (HCP) is known as the mostly efficient way a space can be filled by an arrangement of spheres ( Figure A1a). Sphere constructing an HCP structure fill up 74% of space. For an open-cell lattice structure from the concept of HCP, the spheres can be visualized to inflate from all direction until they create a second lattice structure made of polyhedrons (demonstrated by black solid lines in Figure A1a). To find the unit cell of the new polyhedral, first we select two rows of spheres on top of each other and the corresponding polyhedral lattice structure surrounding them ( Figure A1b). After removal of the spheres we are left with six polyhedral cells, in which we can find a unit cell which after tessellation in space can create the HCP-based polyhedral lattice structure (the red unit cell highlighted in Figure A1c). Dimensions and loads acting on a hexagonal packing unit cell are shown in Figure A1d. The initial spheres demonstrated in Figure A1a create an angle θ for which we have cos θ = 1 3 and sin θ = 2 √ 2 3 . The cross-sectional area of the hexagonal packing unit cell is therefore A UC = 3 √ 3 2 l 2 sin 2 θ. This structure is composed of two types of vertical struts OA (type A) and BC (type C) and one type of inclined strut OB (type B). In each unit cell, there are two type A struts and 6 × 1/3 = 2 type B struts. The term 1/3 is due to the fact that each type B strut is shared by three adjacent unit cells. Displacement of each of the three types of struts are as follows: (i) Type-A strut: Under external load P, the change in the length of strut OA is (Equation (2) in [35]) (ii) Type-C strut: Under external load P , the change in the length of strut BC is (Equation (iii) Type-B strut OB: Displacements of B with respect to O in directions respectively parallel and perpendicular to strut OB are (see Figure 2a,c in the main paper) Total displacement of C with respect to O in the z and x directions are respectively: Therefore, the strains in the z and x directions are: On the other hand, we know that due to continuity of the material and symmetry of each unit cell with respect to the neighbouring cell: δ A,z = ∆ C,z . This means that Poisson's ratio can be simply found by And normalized elastic modulus can be found by Timoshenko beam theory can be obtained (see Table 1 in the main paper).

Appendix D. Stress and Strain Contours in the Lattice Structures
Appl. Sci. 2021, 11, x FOR PEER REVIEW 32 of 38

Appendix E. Effect of Considering Shear Deformation on the Forces/Moments of a Single Strut
To show the shear deformation effect on the final results of analytical relationships of a single strut, we defined different ratios each being defined as the ratio of the resultant forces and moments (required to create a displacement without rotation or rotation without displacement) in the free end of a cantilever beam based on Euler-Bernoulli to the resultant forces and moments required based on Timoshenko beam theory. The ratios , and respectively represent the ratio of the following parameters calculated based on the Euler-Bernoulli and Timoshenko beam theories: the force required for lateral displacement without rotation, , the moment required for lateral displacement without rotation (and the force required to create rotation without displacement), , and the moment required to create rotation without displacement, . The relationships for ratios are derived below:

Appendix E. Effect of Considering Shear Deformation on the Forces/Moments of a Single Strut
To show the shear deformation effect on the final results of analytical relationships of a single strut, we defined different α i ratios each being defined as the ratio of the resultant forces and moments (required to create a displacement without rotation or rotation without displacement) in the free end of a cantilever beam based on Euler-Bernoulli to the resultant forces and moments required based on Timoshenko beam theory. The ratios α 1 , α 2 and α 3 respectively represent the ratio of the following parameters calculated based on the Euler-Bernoulli and Timoshenko beam theories: the force required for lateral displacement without rotation, T, the moment required for lateral displacement without rotation (and the force required to create rotation without displacement), V, and the moment required to create rotation without displacement, U. The relationships for α i ratios are derived below: The variation of ratios α 1 , α 2 and α 3 versus the parameter r/l is presented in the Figure A6. The first important result revealed from the figure is that α 1 and α 2 are equal for all values of r/l. Therefore, considering the shear deformation effect has the same effect on the force required to create lateral displacement without rotation, T, on the one hand and moment required to create lateral displacement without rotation (and the force required to create rotation without displacement), V, on the other hand. On the other hand, the shear deformation has smaller effect on α 3 in comparison with α 1 and α 2 . Although the r/l is not the ultimate parameter for evaluating the effect of shear deformation in lattice structures for different relative densities, but it is a good measure to predict the difference between the resultant force and moments obtained based on Timoshenko and Euler-Bernoulli beam theories. According to Figure A6, without considering the shear deformation effect in the beam theory, the forces and moments required to create a particular deformation could be predicted by 15-20% higher for r/l as large as 0.15. hand, the shear deformation has smaller effect on in comparison with and . Although the / is not the ultimate parameter for evaluating the effect of shear deformation in lattice structures for different relative densities, but it is a good measure to predict the difference between the resultant force and moments obtained based on Timoshenko and Euler-Bernoulli beam theories. According to Figure A6, without considering the shear deformation effect in the beam theory, the forces and moments required to create a particular deformation could be predicted by 15-20% higher for / as large as 0.15.