Effective Mechanical Properties of Periodic Cellular Solids with Generic Bravais Lattice Symmetry via Asymptotic Homogenization

In this paper, the scope of discrete asymptotic homogenization employing voxel (cartesian) mesh discretization is expanded to estimate high fidelity effective properties of any periodic heterogeneous media with arbitrary Bravais’s lattice symmetry, including those with non-orthogonal periodic bases. A framework was developed in Python with a proposed fast–nearest neighbour algorithm to accurately estimate the periodic boundary conditions of the discretized representative volume element of the lattice unit cell. Convergence studies are performed, and numerical errors caused by both voxel meshing and periodic boundary condition approximation processes are discussed in detail. It is found that the numerical error in periodicity approximation is cyclically dependent on the number of divisions performed during the meshing process and, thus, is minimized with a refined voxel mesh. Validation studies are performed by comparing the elastic properties of 2D hexagon lattices with orthogonal and non-orthogonal bases. The developed methodology was also applied to derive the effective properties of several lattice topologies, and variation of their anisotropic macroscopic properties with relative densities is presented as material selection charts.


Introduction
Periodic cellular solids, also known as lattice materials, are formed by tessellating a unit cell in an infinite periodicity to fill a design space.The unit cell, the representative volume element (RVE), is the smallest depiction of the lattice that has the most accurate statistical representation of its physical macroscale properties.The dimensions of the unit cell are at least an order of magnitude smaller than the characteristic length of the macroscopic structure.
Multiscale numerical simulations to evaluate lattice material characteristics considering structural attributes at all length scales are normally challenging and computationally expensive.Thus, homogenization methods are developed to analyze the effective qualities of a multiscale material such that its macroscale properties may be represented statistically.
Some notable works [35][36][37] have provided closed-form representations of the effective mechanical characteristics of lattice materials.Their methods, however, were only applicable to simple topology with straightforward arrangement of lattice cells.
The characteristics of planar lattice materials have also been homogenized using matrixbased methods based on Bloch's theorem and the Cauchy-Born hypothesis [24,25].In order to derive the material's macroscopic stiffness parameters, Hutchinson and Fleck [24] first converted the microscopic nodal deformations of a lattice in terms of the macroscopic strain field.A method was developed to describe cell topologies, such as the Kagome lattice and the Triangular-Triangular lattice, which have a particular degree of symmetry.This approach was expanded by ElSayed and Pasini [25] and ElSayed [38] to handle planar topologies that can have any arbitrary cell geometry.A more generic matrix-based method for the analysis of arbitrary bi-dimensional and tri-dimensional cell topologies with open and closed cells was described by Vigliotti and Pasini [23,26,39].
The characterization of cellular materials has also been effectively accomplished using discrete homogenization approaches [27][28][29].These methods simulate the lattice cell walls using discrete elements like beam or rod elements.The homogenized characteristics are produced by converting the discrete sum of equilibrium equations into a continuous relation of stress and strain.
The asymptotic homogenization (AH) theory has been effectively used among other numerical methods for predicting the effective mechanical characteristics of periodic lattice materials [30][31][32][33][34][40][41][42][43][44].Their results have been validated by experimental tests, demonstrating the effectiveness of the AH method [45][46][47][48].The main discrepancy between the experimental and the theoretical analyses was mainly due to defects and deviations present in the manufactured part.The defects invalidated the periodic boundary condition and caused deviations in the measured experimental values.A standout benefit of AH in comparison to other homogenization systems is its ability to precisely identify the stress distribution in the lattice unit cell, which can then be utilized for an in-depth analysis of the strength and damage of heterogeneous periodic materials [29,49].
The double-scale AH method, considering the use of 2D iso-parametric plate elements, was employed by Andreassen and Andreasen [33] to determine the elastic properties, thermal expansion, thermal conductivity, and fluid permeability of 2D periodic cellular lattices and composite materials.Andreassen and Andreasen [33] were able to analyze RVE's cell envelope with monoclinic, orthorhombic, tetragonal, and hexagonal periodicities by changing the shape of the 2D iso-parametric quadrilateral plate element.
Considering the use of solid elements for discrete homogenization techniques, the homogenized elastic tensor could be obtained by applying unit strains and performing volume averaging of the stresses and strains in the discretized elements.This approach has been performed using popular commercial finite element packages, such as Abaqus [50] and ANSYS Material Designer [51].The main limitation for these commercial applications adopting this approach is that they are only developed for RVEs with orthorhombic, tetragonal, or cubic periodicity [50][51][52] mainly due to the limitations imposed by the meshing process, leading to an incompatibility while applying the periodic boundary conditions.When discretizing the RVE, the locations and the number of nodes along a periodic boundary face would not necessarily match while applying periodic boundary conditions.Dong et al. [34] expanded the work of Andreassen and Andreasen [33] by considering the use of 3D solid elements (iso-parametric hexahedral element shown in Figure 1, also hereafter referred to as 3D voxel).The code developed by Dong et al. [34] was limited to analyzing the periodic cell envelope with the orthogonal periodic basis.hereafter referred to as 3D voxel).The code developed by Dong et al. [34] was limited to analyzing the periodic cell envelope with the orthogonal periodic basis.This paper expands on the work performed by Dong et al. [34].The prior work only considered cell envelopes with orthogonal periodicity, whereas the current work expands upon this by (a) proposing methods to discretize unit cells with non-orthogonal bases and (b) applying approximated periodic boundary conditions for non-orthogonal cell envelopes to evaluate the homogenized elastic properties.
The paper is organized into five sections.After the introduction, in Section 2, the methodology for discretizing the geometry and approximating the periodic boundary condition, along with the numerical homogenization process for determining high fidelity homogenized elastic properties, are presented.In Section 3, the validation for the numerical analysis is detailed.The developed methodology is then applied to triclinic and monoclinic Bravais grid lattice topologies, and the results are documented in Section 4. The paper is concluded in Section 5.
The developed methodology was also used to derive the effective properties of 35 lattice topologies, and variation of their anisotropic macroscopic properties with relative densities which is presented in the appendices included in the Supplementary File.

Methodology
The periodic cellular solids can have an open or a closed cell construction.The former can be modelled as a micro-truss-like structure, while the latter is commonly represented with shells and plates.In this paper, an open-cell micro-truss with circular cross-section elements is considered.
The numerical homogenization performed in this paper is based on the asymptotic double-scale homogenization [30][31][32][33]53,54] with a workflow written in Python [55].Here, the periodic cellular lattices are categorized, based on the shape of the cell envelope, as either orthogonal or non-orthogonal.Cell envelopes with orthogonal periodic bases could This paper expands on the work performed by Dong et al. [34].The prior work only considered cell envelopes with orthogonal periodicity, whereas the current work expands upon this by (a) proposing methods to discretize unit cells with non-orthogonal bases and (b) applying approximated periodic boundary conditions for non-orthogonal cell envelopes to evaluate the homogenized elastic properties.
The paper is organized into five sections.After the introduction, in Section 2, the methodology for discretizing the geometry and approximating the periodic boundary condition, along with the numerical homogenization process for determining high fidelity homogenized elastic properties, are presented.In Section 3, the validation for the numerical analysis is detailed.The developed methodology is then applied to triclinic and monoclinic Bravais grid lattice topologies, and the results are documented in Section 4. The paper is concluded in Section 5.
The developed methodology was also used to derive the effective properties of 35 lattice topologies, and variation of their anisotropic macroscopic properties with relative densities which is presented in the appendices included in the Supplementary File.

Methodology
The periodic cellular solids can have an open or a closed cell construction.The former can be modelled as a micro-truss-like structure, while the latter is commonly represented with shells and plates.In this paper, an open-cell micro-truss with circular cross-section elements is considered.
The numerical homogenization performed in this paper is based on the asymptotic double-scale homogenization [30][31][32][33]53,54] with a workflow written in Python [55].Here, the periodic cellular lattices are categorized, based on the shape of the cell envelope, as either orthogonal or non-orthogonal.Cell envelopes with orthogonal periodic bases could be part of an orthorhombic, tetragonal, or cubic Bravais lattice system, as shown in Figure 2. On the other hand, non-orthogonal cell envelopes have non-orthogonal periodic bases and could be part of triclinic, monoclinic, and hexagonal systems.
be part of an orthorhombic, tetragonal, or cubic Bravais lattice system, as shown in Figure 2. On the other hand, non-orthogonal cell envelopes have non-orthogonal periodic bases and could be part of triclinic, monoclinic, and hexagonal systems.

Representative Volume Element
The RVE of a lattice is its unit cell represented by a cell envelope and the lattice structure, as shown in Figure 3.This figure shows the honeycomb lattice topology with different possibilities of cell envelope representation, including both orthogonal and non-orthogonal bases.Figure 3b shows a discretized 2D hexagon geometry with a non-orthogonal cell envelope.Three different types of voxels are shown in this figure.The green voxels represent the voids within the RVE's cell envelope, and they have a non-zero volume and zero material property.The green voxels are only considered when calculating total volume but ignored when finding periodic node pairs and calculating the stiffness matrix.The red voxels represent the discretized geometry, with both non-zero material property and non-zero volume.The blue voxels are assigned a zero volume and zero material property because they are outside the RVE's cell envelope.The cell envelope is defined as a face with a normal pointing outside the cell envelope, as shown in Figure 4. Based on the cell envelope's plane definition, any voxel elements that are outside the cell envelope (same side as the normal vector) are removed during the analysis.

Representative Volume Element
The RVE of a lattice is its unit cell represented by a cell envelope and the lattice structure, as shown in Figure 3.This figure shows the honeycomb lattice topology with different possibilities of cell envelope representation, including both orthogonal and non-orthogonal bases.Figure 3b shows a discretized 2D hexagon geometry with a non-orthogonal cell envelope.Three different types of voxels are shown in this figure.The green voxels represent the voids within the RVE's cell envelope, and they have a non-zero volume and zero material property.The green voxels are only considered when calculating total volume but ignored when finding periodic node pairs and calculating the stiffness matrix.The red voxels represent the discretized geometry, with both non-zero material property and non-zero volume.The blue voxels are assigned a zero volume and zero material property because they are outside the RVE's cell envelope.The cell envelope is defined as a face with a normal pointing outside the cell envelope, as shown in Figure 4. Based on the cell envelope's plane definition, any voxel elements that are outside the cell envelope (same side as the normal vector) are removed during the analysis.The parent material property for the individual voxel is defined using the two L parameters, namely,  and , which are expressed as: where  and  are the Young's modulus and the Poisson's ratio, respectively.Th ysis presented in this paper employs isotropic steel with  = 2.0 × 10 Pa and 0.3, which correspond to lame parameters of  = 1.153 × 10 Pa and  = 7.692 Pa.
The voxel element is formulated as an iso-parametric brick element whose ed orthogonal with lengths of  ,  and  aligned in the global cartesian coordina tem.The individual voxels with non-zero volume, which are inside the cell envelo treated as having an isotropic material property, and the element's stiffness mater trix ( ( ) ) is formulated as: The parent material property for the individual voxel is defined using the two Lame's parameters, namely, λ and µ, which are expressed as: where E and ν are the Young's modulus and the Poisson's ratio, respectively.The analysis presented in this paper employs isotropic steel with E iso = 2.0 × 10 11 Pa and ν iso = 0.3, which correspond to lame parameters of λ = 1.153 × 10 11 Pa and µ = 7.692 × 10 11 Pa.
The voxel element is formulated as an iso-parametric brick element whose edges are orthogonal with lengths of l x , l y and l z aligned in the global cartesian coordinate system.The individual voxels with non-zero volume, which are inside the cell envelope, are treated as having an isotropic material property, and the element's stiffness material matrix (C (e) ) is formulated as: 2 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1

Discretization and Approximation of the Periodic Boundary Condition
The process to determine the translational periodicity of the node pairs is to translate the coordinates of a node using the periodic basis vector and then find the closest node within a prescribed search radius.
A periodic physical quantity of a structure can be represented as: where x = [x 1 , x 2 , x 3 ] T is the position vector of a point where the physical quantity F is evaluated, N = n 1 , n 2 , n 3 is a 3 × 3 diagonal matrix which consists of arbitrary integer values, and Y = [Y 1 , Y 2 , Y 3 ] T is a vector that denotes the period of the structure, where this value could be a scalar, vector, or tensor function of x [30].
The static response of a periodic structure can then be represented as: where C ijkl is the stiffness tensor, ε kl is the strain tensor, and σ ij is the stress tensor.By definition of x, N, and Y, the stiffness tensor C ijkl can be expanded as The voxel nodes have been discretized such that the nodal coordinates (x i 1 , x i 2 , x i 3 ) of the ith node are defined as: where, A 1 , A 2 , and A 3 are integers and {O x , O y , O z } is the origin.The i superscript indicates node numbers; l x , l y , and l z are the lengths of the voxel along the global x, y, and z axis, respectively; and the x divs , y divs and z divs are integers denoting the number of voxel discretization along the global x, y, and z axes, respectively.Typically, l x , l y , and l z are the same for all the voxels and can be calculated as: The periodicity vector Y can be decomposed such that: where B 1 , B 2 , and B 3 are integers and G 1 , G 2 , and G 3 are the residual or the remainder which are smaller than the lengths of the voxel l x , l y , and l z , respectively.
For each node of the voxel defined by Equation ( 6) that is part of the RVE's cell envelope, there exists a periodic node pair that is also part of the discretized RVE's cell envelope.
If any of the residuals G 1 , G 2 , or G 3 are non-zero, then it would not be possible to find the exact node that is part of the voxels, as shown in Figure 5, (which leads to approximating the periodic node pairs, causing some numerical errors which are discussed in Sections 3.1 and 3.4).Approximating the location of the periodic node pair involves ignoring the residuals (G i ).A perfect voxel mesh can be obtained by minimizing the residual term.For an orthogonal periodic basis, aligned with the global cartesian coordinate system, the residuals are automatically zero.But for an RVE with a non-orthogonal periodic basis, the x divs , y divs , and z divs must be altered to minimize the residuals.
For a coarse mesh (where  ,  ,  values are similar or larger relative to the truss radius), it is suggested to use angular constraints, where the angle  , shown in Figure 5, is more than 5°, rather than the suggested radius constraint in Equation (9).The angle  , as shown in Figure 5, can be calculated by performing the dot product operation between the two normalized vectors,  ⃗ and   ⃗ .The misalignment of the periodic basis with the voxel's element natural coordinate system introduces slight numerical errors, which are discussed in Section 3.
Figure 5. Determination of a periodic node pair for node  located at ( ,  ) with its periodic pair ′ at   ⃗ and the search radius  ⃗ (green arrow) required to search the approximate periodic node pair due to the difference between the voxel element basis and RVE's periodic basis.This search radius can be used with a KD-Tree algorithm [56] to query the closest point.
To determine the periodicity of the lattice, the nodal coordinates of all the voxels are used to create a k-d tree from [56] to implement a fast-nearest neighbour algorithm.Then, Figure 5. Determination of a periodic node pair for node i located at (x i 1 , x i 2 ) with its periodic pair i at N i → Y i and the search radius → R (green arrow) required to search the approximate periodic node pair due to the difference between the voxel element basis and RVE's periodic basis.This search radius can be used with a KD-Tree algorithm [56] to query the closest point.
Thus, to find the periodic node pair for an RVE with a non-orthogonal periodic basis, a local search must be conducted within a radius of R, such that at least one point lies within the defined search radius, as shown in Figure 5.If the periodic node is not found for the node in the cell boundary, or the node pair is incorrectly identified, it can introduce small numerical errors.A suggested search radius is proposed as follows: R = max 0.51 max l x , l y , l z 0.51 For a coarse mesh (where l x , l y , l z values are similar or larger relative to the truss radius), it is suggested to use angular constraints, where the angle θ p , shown in Figure 5, is more than 5 • , rather than the suggested radius constraint in Equation (9).The angle θ p , as shown in Figure 5, can be calculated by performing the dot product operation between the two normalized vectors, → Y and → x i p.The misalignment of the periodic basis with the voxel's element natural coordinate system introduces slight numerical errors, which are discussed in Section 3.
To determine the periodicity of the lattice, the nodal coordinates of all the voxels are used to create a k-d tree from [56] to implement a fast-nearest neighbour algorithm.Then, similar to Equation ( 5), an offset from Equation ( 8) is applied to all of the voxel coordinates based on all of the unique combinations of the basis, and the closest point within a certain minimum distance R and a minimum angle θ is chosen as its periodic nodal pair.Unit cells with the orthogonal periodic basis that align with the voxel's natural coordinate system do not need the minimum distance and the angle restrictions.But for a lattice geometry with a periodic basis that is not an integer multiple of the voxel's natural coordinate system, the process of finding the periodic nodal pair requires an additional step because the offset coordinate is not coincident with the second nodal pair.Thus, a local search must be conducted to find the closest point that is subject to minimum distance and angle restrictions described by (9).
pq is defined as: Based on displacement fields χ kl , which are found by solving the elasticity equations with prescribed macroscopic strains: where ν is the virtual displacement field.Homogenization is normally performed numerically using discretization and finite element method to solve Equation ( 12).

Validation
In this section, the analysis of a voxelized RVE with a non-orthogonal base is validated against results of a commercial finite element software for the same lattice topology but represented by a different RVE with orthogonal bases, as reported by Gibson and Ashby [36], Vigliotti and Pasini [26], and other homogenization codes by Andreassen and Andreasen [33], as well as Dong et al. [34].For low relative density lattice, the results of the voxelized non-orthogonal hexagon are compared against similar non-orthogonal hexagons but formulated using beam theory [27][28][29].Grid convergence studies have been performed in this section, and the numerical errors caused by both the voxel meshing and the periodic boundary condition approximation processes are discussed in detail.The material used in the validation study is isotropic steel with E iso = 2.0 × 10 11 Pa and ν iso = 0.3.

Grid Convergence Study
For the grid convergence study, the lattice materials' truss radius is held constant, and the number of divisions along the global axis is varied.For the first mesh sensitivity study, a 3D primitive cubic lattice with an orthogonal cell envelope was analyzed, and the results are shown in Figure 6.The primitive cubic lattice, which had a unit length of 1 with a truss radius of 0.157, was targeted, which converged to a volume fraction of 0.18.The convergence of elastic properties was satisfied around a discretization of 60 voxels in each direction.Then, a 2D closed hexagon lattice was analyzed, and a volume fraction of 0.25 was targeted using a radius of 0.117 units.The length of the truss was 1 unit, and the cell angle was 60°.The elastic constants of the stiffness matrix and the elastic properties are presented in Figure 7 It was observed that increasing the number of divisions affects the volume fraction of the discretized geometry; this is due to the voxel meshing process, where at some critical point, the volume fraction increases due to the voxel center coordinate suddenly being included as part of the truss radius (refer to Figure 8).In the plots below, the other quantities generally follow a similar trend as the volume fraction (black dashed lines).However, the observed discrepancy is generally due to the errors in approximating periodicity, which generally converges while refining the mesh.Then, a 2D closed hexagon lattice was analyzed, and a volume fraction of 0.25 was targeted using a radius of 0.117 units.The length of the truss was 1 unit, and the cell angle was 60 • .The elastic constants of the stiffness matrix and the elastic properties are presented in Figure 7 It was observed that increasing the number of divisions affects the volume fraction of the discretized geometry; this is due to the voxel meshing process, where at some critical point, the volume fraction increases due to the voxel center coordinate suddenly being included as part of the truss radius (refer to Figure 8).In the plots below, the other quantities generally follow a similar trend as the volume fraction (black dashed lines).However, the observed discrepancy is generally due to the errors in approximating periodicity, which generally converges while refining the mesh.Then, a 2D closed hexagon lattice was analyzed, and a volume fraction of 0.25 was targeted using a radius of 0.117 units.The length of the truss was 1 unit, and the cell angle was 60°.The elastic constants of the stiffness matrix and the elastic properties are presented in Figure 7 It was observed that increasing the number of divisions affects the volume fraction of the discretized geometry; this is due to the voxel meshing process, where at some critical point, the volume fraction increases due to the voxel center coordinate suddenly being included as part of the truss radius (refer to Figure 8).In the plots below, the other quantities generally follow a similar trend as the volume fraction (black dashed lines).However, the observed discrepancy is generally due to the errors in approximating periodicity, which generally converges while refining the mesh.The theoretical zero terms presented in an orthotropic stiffness matrix are plotted in Figure 7 above.It can be observed that the supposed zero terms fluctuate between a high and low value.The higher values are at least two orders of magnitude smaller than the main diagonals of the matrix, whereas the lower error terms are at least eight orders of magnitude smaller than the main diagonals.The main cause for this is because the periodicity is discretized where the remainder term ( ) mentioned in Equation ( 8) is minimized; this means that the periodicity basis vector could be resolved as an integer multiple of the voxel lengths.
The effect of reduced remainder term leads to a case where the periodic nodes are matched in an orderly fashion, as shown in Figures 8 and 7b.When comparing the periodicity lines from Figure 7a's blue circle to the magenta circle, the periodicity lines are not parallel, as shown in Figure 7b.Furthermore, comparing the periodicity from the blue The theoretical zero terms presented in an orthotropic stiffness matrix are plotted in Figure 7 above.It can be observed that the supposed zero terms fluctuate between a high and low value.The higher values are at least two orders of magnitude smaller than the main diagonals of the matrix, whereas the lower error terms are at least eight orders of magnitude smaller than the main diagonals.The main cause for this is because the periodicity is discretized where the remainder term (G i ) mentioned in Equation ( 8) is minimized; this means that the periodicity basis vector could be resolved as an integer multiple of the voxel lengths.
The effect of reduced remainder term leads to a case where the periodic nodes are matched in an orderly fashion, as shown in Figures 8 and 7b.When comparing the periodicity lines from Figure 7a's blue circle to the magenta circle, the periodicity lines are not parallel, as shown in Figure 7b.Furthermore, comparing the periodicity from the blue circle to the red circle, a single periodic node is shared with two other nodes; this can lead to the periodicity being mismatched.Several of these scenarios are visualized in Figure 9. Based on the visualization, the planes of symmetry are disrupted when periodicity is not matched properly, which causes the orthotropic cellular lattices to have smaller magnitude non-zero terms like a monoclinic or a triclinic material.An in-depth analysis is presented in Section 3.4.circle to the red circle, a single periodic node is shared with two other nodes; this can lead to the periodicity being mismatched.Several of these scenarios are visualized in Figure 9.
Based on the visualization, the planes of symmetry are disrupted when periodicity is not matched properly, which causes the orthotropic cellular lattices to have smaller magnitude non-zero terms like a monoclinic or a triclinic material.An in-depth analysis is presented in Section 3.4.The periodic node pairs are connected using randomly colored lines.

Comparison with ANSYS
The results obtained from the voxelization process are compared with the results from ANSYS Material Designer, which can perform homogenization of elastic and thermal properties for RVE with an orthogonal periodic basis.The results obtained from AN-SYS are comparable with the results obtained with the voxelized method, with differences of less than 2%.The anisotropy plot for the primitive cubic lattice and hexagon lattice is shown in Figure 10.The primitive cubic lattice in ANSYS was meshed using tetrahedrons,

Comparison with ANSYS
The results obtained from the voxelization process are compared with the results from ANSYS Material Designer, which can perform homogenization of elastic and thermal properties for RVE with an orthogonal periodic basis.The results obtained from ANSYS are comparable with the results obtained with the voxelized method, with differences of less than 2%.The anisotropy plot for the primitive cubic lattice and hexagon lattice is shown in Figure 10.The primitive cubic lattice in ANSYS was meshed using tetrahedrons, and the meshing strategy made sure that the meshes in all of the faces were similar, such that the periodicity of the lattice could be found easily.The hexagon lattice in ANSYS follows a similar meshing strategy, but the Honeycomb was meshed with a mix of hexahedral (brick) and tetrahedral elements.However, the cell envelope's faces are meshed using brick elements, which ensures that the cell envelope's nodes can be matched easily to apply the periodic boundary conditions.ANSYS also provides an option to calculate the homogenized properties by applying symmetric boundary conditions for some of the simpler lattices.Three planes of symmetry exist in the x, y, and z-planes for both the primitive cubic lattice in 3D and the hexagon lattice in 2D (with periodic conditions imposed in the outof-plane).Thus, the homogenized elastic properties for these two lattices can be considered as an orthotropic.Theoretically, the homogenized elasticity tensor in the Voigt notation for an orthotropic lattice cell takes the form of: Due to some small numerical residuals, the zero terms shown in Equation ( 13) were not zeros when the hexagon and the cubic lattice were homogenized.A comparison of the terms for  for the hexagon and cubic lattice is shown in Figure 11.The data exported Three planes of symmetry exist in the x, y, and z-planes for both the primitive cubic lattice in 3D and the hexagon lattice in 2D (with periodic conditions imposed in the out-ofplane).Thus, the homogenized elastic properties for these two lattices can be considered as an orthotropic.Theoretically, the homogenized elasticity tensor in the Voigt notation for an orthotropic lattice cell takes the form of: Due to some small numerical residuals, the zero terms shown in Equation ( 13) were not zeros when the hexagon and the cubic lattice were homogenized.A comparison of the terms for C H 41 for the hexagon and cubic lattice is shown in Figure 11.The data exported from ANSYS, which is marked as "D [1,1]" and "D [1,4]" in the figures, are from the grid and the hexagon lattice cell, respectively.Based on the data presented in Figure 11, it can be observed that results from ANSYS have a lower  term for the hexagon lattice, and a higher error for the grid lattice, in comparison to the voxelized homogenization process developed in this paper.This small numerical error is caused by the meshing scheme, as the grid lattice was meshed using tetrahedrons, while the hexagons were meshed mostly using brick elements, as shown in Figure 10c,d, respectively.In the next section, numerical errors caused by approximating periodicity are discussed.

RVE Rotation
The anisotropic behaviour of the elastic properties of the lattice cell is investigated by rotating the microstructure of the unit cell with respect to the applied macroscopic strain field; this is carried out by applying the rotation tensor as follows: where  is the orthogonal tensor, which corresponds to an orthogonal transformation from  to  basis.The    in the Voigt notation is expanded to    , and an Einstein summation [57] was performed, where   =   =   =   = .The 2D anisotropy plot produced in this paper (shown in Figure 10) corresponds to a rotation in the  the axis of the global coordinate system, where: The 3D anisotropy plot produced in this paper corresponds to a rotation in the  axis and the  axis, where: cos () −sin () 0 cos () 0 sin () Based on the data presented in Figure 11, it can be observed that results from ANSYS have a lower C H 41 term for the hexagon lattice, and a higher error for the grid lattice, in comparison to the voxelized homogenization process developed in this paper.This small numerical error is caused by the meshing scheme, as the grid lattice was meshed using tetrahedrons, while the hexagons were meshed mostly using brick elements, as shown in Figure 10c,d, respectively.In the next section, numerical errors caused by approximating periodicity are discussed.

RVE Rotation
The anisotropic behaviour of the elastic properties of the lattice cell is investigated by rotating the microstructure of the unit cell with respect to the applied macroscopic strain field; this is carried out by applying the rotation tensor as follows: where Q is the orthogonal tensor, which corresponds to an orthogonal transformation from x i to x i basis.The C H ij in the Voigt notation is expanded to E H ijkl , and an Einstein summation [57] was performed, where The 2D anisotropy plot produced in this paper (shown in Figure 10) corresponds to a rotation in the x 3 the axis of the global coordinate system, where: The 3D anisotropy plot produced in this paper corresponds to a rotation in the x 3 axis and the x 2 axis, where: Using the compliance matrix by finding the inverse of the elasticity matrix, the following elastic properties can be determined: where S ij is the compliance matrix, and E ij is the elastic stiffness matrix.The sensitivity of the primitive cubic lattice's elastic property to the cell rotation was performed and presented in Figure 12.The radius for each plot was adjusted such that a volume fraction of 0.3 could be achieved.The RVE was rotated along the z-axis by an angle of α, and the E 11 was plotted in the anisotropy plot, which was plotted for α = 0 • .For a volume fraction of 0.3, it was noted that 50 voxel discretization was sufficient to prove that applying the rotation matrix to the homogenized elastic tensor is equivalent to rotating the RVE.
terials 2023, 16, x FOR PEER REVIEW 14 of 33 where  is the compliance matrix, and  is the elastic stiffness matrix.The sensitivity of the primitive cubic lattice's elastic property to the cell rotation was performed and presented in Figure 12.The radius for each plot was adjusted such that a volume fraction of 0.3 could be achieved.The RVE was rotated along the z-axis by an angle of , and the  was plotted in the anisotropy plot, which was plotted for  = 0°.For a volume fraction of 0.3, it was noted that 50 voxel discretization was sufficient to prove that applying the rotation matrix to the homogenized elastic tensor is equivalent to rotating the RVE.

Numerical Errors Due to Approximating Periodicity
In this section, a hexagon was meshed with an equal number of divisions along the global x and y axes and one element along the z-axis.The following term ( ) has been visualized in Figure 14:

Numerical Errors Due to Approximating Periodicity
In this section, a hexagon was meshed with an equal number of divisions along the global x and y axes and one element along the z-axis.The following term (τ e ) has been visualized in Figure 14: where χ and χ (i) e are the element's macroscopic and microscopic displacement due to unit strain i, and τ e represents the contribution of element e to the homogenized elastic tensor C H ij normalized by the total C H ij .By visualizing τ e for each element, it is possible to observe the elements that could cause the numerical errors shown in Figure 11.The colour scheme used in Figure 11 is such that the maximum absolute value dictates the positive (red) and negative (blue) limits of the plot.Values close to zero are given a green colour.
where  ( ) and  ( ) are the element's macroscopic and microscopic displacement due to unit strain , and  represents the contribution of element  to the homogenized elastic tensor  normalized by the total  .By visualizing  for each element, it is possible to observe the elements that could cause the numerical errors shown in Figure 11.The colour scheme used in Figure 11 is such that the maximum absolute value dictates the positive (red) and negative (blue) limits of the plot.Values close to zero are given a green colour.For the hexagon closed topology with a volume fraction of 0.29 using isotropic steel with  = 2.0 × 10 and  = 0.3, which has been discretized as 30 × 30 × 1 voxels ( / = 1.155), the homogenized elastic tensor was computed to be: After the mesh resolution was increased to 150 × 150 × 1 voxels ( / = 1.155), with a volume fraction of 0.30, the refined, homogenized elastic tensor was computed to For the hexagon closed topology with a volume fraction of 0.29 using isotropic steel with E iso = 2.0 × 10 11 and ν iso = 0.3, which has been discretized as 30 × 30 × 1 voxels (l x /l y = 1.155), the homogenized elastic tensor was computed to be: After the mesh resolution was increased to 150 × 150 × 1 voxels (l x /l y = 1.155), with a volume fraction of 0.30, the refined, homogenized elastic tensor was computed to be: The contributions of the voxels that contribute the small non-zero terms for the 30 × 30 × 1 and the 150 × 150 × 1 are shown in Figure 14a b, respectively.The comparison of the eigen values with and without the zero terms is tabulated in Table 1.Comparing the hexagon corresponding to C H 51 in Figure 14a with Figure 15a, it can be seen that the voxels that have the highest contributions belong to the group whose periodicity has not been matched properly due to discretizing the periodicity basis vector; this can be seen in Figure 15a-the line that marks the periodic node pair for the nodes near the marked area of a 1 to a 5 , is not parallel with the periodic basis (bold green line), and a single node from a 1 is periodic, with two nodes from a 5 zone.This behaviour causes over-stiffening of the nodal degree of freedom because the stiffness corresponding to another node would be accumulated into one node during the global stiffness formulation.Similarly, if the nodes on the cell envelope are not periodically matched with another node on the envelope, it would cause the stiffness of the nodal degree of freedom to be smaller than the ones that have been matched.The over-stiffening and the under-stiffening impact the microscopic displacements, which contribute to values that are above or below the actual threshold.The mismatch in the periodic degrees of freedom would remove planes of symmetry and would cause the RVE with hexagonal symmetry to be represented as a monoclinic or a triclinic material, as shown in Figure 9.    Furthermore, it can be observed from Figures 14 and 16 that increasing the voxel discretization reduces the small non-zero components in the C H ij tensor because each voxel contributes a smaller value due to the smaller voxel volume.In Figure 14, for the hexagon with a larger voxel size, the residual terms are contributed by the voxels along the edges caused by imperfect periodic boundary conditions; this is also the case for the refined hexagon lattice, but the ratio of voxels in the edge is smaller than the hexagon with larger voxel size, so the contribution of the residual error is decreased.

Comparison of Young's Modulus with Volume Fraction
The  elastic property for the hexagon lattice that is obtained from the voxelization code and ANSYS are compared with different homogenization procedures obtained from Gibson and Ashby [36] and Vigliotti and Pasini [26].In Figure 17, the  , properties calculated from the voxelized closed cell hexagon, open cell hexagon, and hexagon cell envelope (all shown in Figure 18) with orthogonal basis are compared.No significant deviations were observed when comparing it to the results from ANSYS, which were obtained from a hexagon cell envelope with an orthogonal basis.

Comparison of Young's Modulus with Volume Fraction
The E 11 elastic property for the hexagon lattice that is obtained from the voxelization code and ANSYS are compared with different homogenization procedures obtained from Gibson and Ashby [36] and Vigliotti and Pasini [26].In Figure 17, the E 11 , properties calculated from the voxelized closed cell hexagon, open cell hexagon, and hexagon cell envelope (all shown in Figure 18) with orthogonal basis are compared.No significant deviations were observed when comparing it to the results from ANSYS, which were obtained from a hexagon cell envelope with an orthogonal basis.

Comparison of Young's Modulus with Volume Fraction
The  elastic property for the hexagon lattice that is obtained from the voxelization code and ANSYS are compared with different homogenization procedures obtained from Gibson and Ashby [36] and Vigliotti and Pasini [26].In Figure 17, the  , properties calculated from the voxelized closed cell hexagon, open cell hexagon, and hexagon cell envelope (all shown in Figure 18) with orthogonal basis are compared.No significant deviations were observed when comparing it to the results from ANSYS, which were obtained from a hexagon cell envelope with an orthogonal basis.

Comparison of Young's Modulus with Volume Fraction
The  elastic property for the hexagon lattice that is obtained from the voxelization code and ANSYS are compared with different homogenization procedures obtained from Gibson and Ashby [36] and Vigliotti and Pasini [26].In Figure 17, the  , properties calculated from the voxelized closed cell hexagon, open cell hexagon, and hexagon cell envelope (all shown in Figure 18) with orthogonal basis are compared.No significant deviations were observed when comparing it to the results from ANSYS, which were obtained from a hexagon cell envelope with an orthogonal basis.A hexagon lattice made from a circular truss (i.e., without considering periodicity in the out-of-plane direction) and a hexagon with a plate wall (3D voxel with a single layer and periodicity in the z-direction) are considered below.Gibson and Ashby's model [36] starts to overpredict the elasticity, whereas Vigliotti and Pasini's [26] model, which uses a beam model, underpredicts the elasticity at higher relative density.

Comparison of Elastic Properties for Hexagon Lattice with Varying Cell Angle
In this section, the cell angle is changed, and Young's modulus and shear modulus are plotted for a hexagon with a relative density of 0.2.The results obtained from the voxel homogenization process match very well, with some slight discrepancies; this is caused by the geometry discretization process, where a volume fraction of 0.2 was not achieved due to the linear interpolation used with the optimizer to define the relative density as a function of the radius of the truss elements.The variation of the homogenized elastic properties for the hexagon lattice w.r.t cell angle is presented in Figure 19.A hexagon lattice made from a circular truss (i.e., without considering periodicity in the out-of-plane direction) and a hexagon with a plate wall (3D voxel with a single layer and periodicity in the z-direction) are considered below.Gibson and Ashby's model [36] starts to overpredict the elasticity, whereas Vigliotti and Pasini's [26] model, which uses a beam model, underpredicts the elasticity at higher relative density.

Comparison of Elastic Properties for Hexagon Lattice with Varying Cell Angle
In this section, the cell angle is changed, and Young's modulus and shear modulus are plotted for a hexagon with a relative density of 0.2.The results obtained from the voxel homogenization process match very well, with some slight discrepancies; this is caused by the geometry discretization process, where a volume fraction of 0.2 was not achieved due to the linear interpolation used with the optimizer to define the relative density as a function of the radius of the truss elements.The variation of the homogenized elastic properties for the hexagon lattice w.r.t cell angle is presented in Figure 19.

Comparison of Elastic Constants for a 2D Monoclinic RVE Lattice with Varying Cell Angle
In this section, a square-like 2D lattice is considered; the angle shown varies from 45° to 90°.The results are compared with code published by Andreassen and Andreasen [33] as shown in Figure 20.Andreassen's work is considered more accurate because the 2D voxel elements change shape as the cell angle changes, and there are no approximations performed when applying periodicity boundary conditions to the lattice.The voxel code underestimates the  component of the stiffness matrix because of the periodicity approximations and the voxel meshing algorithm.The main cause is the voxel meshing, where the truss elements are jagged, as shown in Figure 8.
Homogenized elastic properties of a square lattice with varying cell angles.

Comparison of Elastic Constants for a 2D Monoclinic RVE Lattice with Varying Cell Angle
In this section, a square-like 2D lattice is considered; the angle shown varies from 45 • to 90 • .The results are compared with code published by Andreassen and Andreasen [33] as shown in Figure 20.Andreassen's work is considered more accurate because the 2D voxel elements change shape as the cell angle changes, and there are no approximations performed when applying periodicity boundary conditions to the lattice.The voxel code underestimates the C 22 component of the stiffness matrix because of the periodicity approximations and the voxel meshing algorithm.The main cause is the voxel meshing, where the truss elements are jagged, as shown in Figure 8.

Comparison of Elastic Constants for a 2D Monoclinic RVE Lattice with Varying Cell Angle
In this section, a square-like 2D lattice is considered; the angle shown varies from 45° to 90°.The results are compared with code published by Andreassen and Andreasen [33] as shown in Figure 20.Andreassen's work is considered more accurate because the 2D voxel elements change shape as the cell angle changes, and there are no approximations performed when applying periodicity boundary conditions to the lattice.The voxel code underestimates the  component of the stiffness matrix because of the periodicity approximations and the voxel meshing algorithm.The main cause is the voxel meshing, where the truss elements are jagged, as shown in Figure 8.
Homogenized elastic properties of a square lattice with varying cell angles.

Application
In this section, application of the developed voxel code is applied to different lattice topologies with non-orthogonal bases.The material used in this study is isotropic steel with  = 2.0 × 10 Pa and  = 0.3.

Application
In this section, application of the developed voxel code is applied to different lattice topologies with non-orthogonal bases.The material used in this study is isotropic steel with E iso = 2.0 × 10 11 Pa and ν iso = 0.3.

Three-Dimensional Triclinic and a Monoclinic Bravais Grid Lattice
In this section, a triclinic and a monoclinic Bravais grid lattice are analyzed by applying the approximated periodicity boundary condition.The geometry was discretized using voxels, as shown in Figure 21.The lattice definition used for the current analysis to create the triclinic and monoclinic Bravais grid lattice is shown in Figure 22.The lattice is defined such that the red trusses are in the global x-y plane.The trusses T 14 and T 23 are aligned parallel to the global x-axis.The other red trusses T 12 and T 43 , are rotated from the global y-axis using angle gamma (γ) in the x-y plane.The green trusses are created by offsetting the red trusses using the blue trusses.The blue trusses are defined using a spherical coordinate system with φ and θ rotation angles, as shown in Figure 22.The periodic basis for the lattice was defined using the trusses T 12 , T 14 , and T 15 .The angular restrictions for φ, θ, and γ for the Bravais lattice system are tabulated in Table 2.

Phi (φ) Theta (θ) Gamma (γ)
The key difference between the triclinic and the monoclinic Bravais grid lattice used for the current analysis is the slight change in the truss angles and their corresponding periodicity bases.It can be noted that for the monoclinic lattice, the trusses that make up the red and the green frames are orthogonal to each other (γ = 0 • ).Furthermore, the point set {1,4,5,8} and {2,3,6,7} are co-planar with the z-x plane.
For the monoclinic Bravais grid lattice shown in Figure 23b, which had a volume fraction of 0.29 and was discretized as 30 × 30 × 1 voxels, the homogenized elastic tensor was computed to be: The elastic properties of the monoclinic and the triclinic lattice are tabulated in Table 3.The evolution of the elastic properties visualized as 2D and 3D anisotropic plots is shown in Figures 24 and 25 for the triclinic and monoclinic Bravais grid lattice, respectively.The key difference between the 3D anisotropy plot between the triclinic and monoclinic Bravais grid lattice is how the peaks of the anisotropy plots are warped.Similarly, for the triclinic Bravais grid lattice shown in Figures 23a and 23, which had a volume fraction of 0.29 and was discretized as 30 × 30 × 1 voxels, the homogenized elastic tensor was computed to be The elastic properties of the monoclinic and the triclinic lattice are tabulated in Table 3.The evolution of the elastic properties visualized as 2D and 3D anisotropic plots is shown in Figures 24 and 25 for the triclinic and monoclinic Bravais grid lattice, respectively.The key difference between the 3D anisotropy plot between the triclinic and monoclinic Bravais grid lattice is how the peaks of the anisotropy plots are warped.

Two-Dimensional Non-Orthogonal Lattice
In this section, the elastic properties of the 3.4.6.4 2D lattice are plotted.The homogenized elastic tensor of the 3.4.6.4 lattices shown in Figure 26 is tabulated in Table 3.The tabulated tensor is based on voxels that are generated as a single layer with an additional periodicity added out of a plane (z-axis).Considering this periodicity means that the single layer is equivalent to a fully tessellated geometry that has infinite depth in the out-of-plane direction.The tabulated tensors are obtained by fitting a third-order polynomial for each of the tensor entries, where the volume fraction varies from 0.1 to 0.9.

Two-Dimensional Non-Orthogonal Lattice
In this section, the elastic properties of the 3.4.6.4 2D lattice are plotted.The homogenized elastic tensor of the 3.4.6.4 lattices shown in Figure 26 is tabulated in Table 3.The tabulated tensor is based on voxels that are generated as a single layer with an additional periodicity added out of a plane (z-axis).Considering this periodicity means that the single layer is equivalent to a fully tessellated geometry that has infinite depth in the out-ofplane direction.The tabulated tensors are obtained by fitting a third-order polynomial for each of the tensor entries, where the volume fraction varies from 0.1 to 0.9. ≅ 0.05  ≅ 0.28  ≅ 0.56  ≅ 0.86   4.

Three-Dimensional Sandwich Panel Lattice
In this section, the elastic properties of the sandwiched X (body-centred cubic) lattice are plotted.For the sandwich panel analysis, periodicity in the z-direction (sandwich plate normal) is assumed; thus, the elastic properties that are plotted in this section are for sandwich panels that are assumed to be stacked.For all the sandwich panel geometries analyzed in this paper, the thickness of the sandwich panel is held constant at a unit of 0.05 of the cell length.The homogenized elastic tensor of the Sandwich X lattice is shown and tabulated in Table 3.The 2D and 3D anisotropy plots for the Sandwich X lattice are shown in Figure 28.The coefficients for the homogenized elastic tensor as a function of relative density is tabulated in Table 5.

Three-Dimensional Sandwich Panel Lattice
In this section, the elastic properties of the sandwiched X (body-centred cubic) lattice are plotted.For the sandwich panel analysis, periodicity in the z-direction (sandwich plate normal) is assumed; thus, the elastic properties that are plotted in this section are for sandwich panels that are assumed to be stacked.For all the sandwich panel geometries analyzed in this paper, the thickness of the sandwich panel is held constant at a unit of 0.05 of the cell length.The homogenized elastic tensor of the Sandwich X lattice is shown and tabulated in Table 3.The 2D and 3D anisotropy plots for the Sandwich X lattice are shown in Figure 28.The coefficients for the homogenized elastic tensor as a function of relative density is tabulated in Table 5.   4.  Table 5. Homogenized elastic tensor coefficients (based on Equation ( 13)) and correlation coefficient (R 2 ) for the polynomial fit for the Sandwich X lattice.

Three-Dimensional Sandwich Panel Lattice
In this section, the elastic properties of the sandwiched X (body-centred cubic) lattice are plotted.For the sandwich panel analysis, periodicity in the z-direction (sandwich plate normal) is assumed; thus, the elastic properties that are plotted in this section are for sandwich panels that are assumed to be stacked.For all the sandwich panel geometries analyzed in this paper, the thickness of the sandwich panel is held constant at a unit of 0.05 of the cell length.The homogenized elastic tensor of the Sandwich X lattice is shown and tabulated in Table 3.The 2D and 3D anisotropy plots for the Sandwich X lattice are shown in Figure 28.The coefficients for the homogenized elastic tensor as a function of relative density is tabulated in Table 5.

Conclusions
The homogenized elastic properties of several 2D and 3D lattices can be performed using the voxel mesh approach regardless of whether the periodicity is orthogonal or nonorthogonal.Approximating the periodicity on an imperfect voxel mesh misaligns the periodic node pairs; this alters or removes the planes of symmetry within the RVE and introduces errors in the stiffness matrix.These errors are at least two orders of magnitude smaller than the main diagonals for the imperfect voxel mesh.The mesh is considered perfectly voxelized if the translated nodes of the cell envelope, using the periodicity basis, align themselves with another node.For this perfectly voxelized mesh, the numerical errors are at least eight orders of magnitude smaller than the main diagonals; this can be achieved by altering the number of discretizations performed along the global coordinate system.Furthermore, the numerical error caused by approximating the non-orthogonal periodicity decreases as the voxel size is reduced.We have shown that it is possible to evaluate the elastic properties of periodic cellular materials whose periodic basis is of any Bravais lattice system.

Figure 1 .
Figure1.Iso-parametric hexahedral voxel element, with the global coordinate system (x, y, z), local coordinate system (ξ, ζ, η), and voxel edges aligned with a global cartesian coordinate system with lengths (l x , l y , l z ), respectively.

Figure 3 .
Figure 3. (a) Honeycomb lattice RVE with multiple cell envelope definitions.Periodic bases are shown with red and green arrows for the three proposed envelopes.(b) Visualization of the 2D Open Hexagon RVE, RVE's envelope, voxels, and periodic basis.Voxels' shape is not to scale.

Figure 3 .
Figure 3. (a) Honeycomb lattice RVE with multiple cell envelope definitions.Periodic bases are shown with red and green arrows for the three proposed envelopes.(b) Visualization of the 2D Open Hexagon RVE, RVE's envelope, voxels, and periodic basis.Voxels' shape is not to scale.

Figure 3 .
Figure 3. (a) Honeycomb lattice RVE with multiple cell envelope definitions.Periodic ba shown with red and green arrows for the three proposed envelopes.(b) Visualization of Open Hexagon RVE, RVE's envelope, voxels, and periodic basis.Voxels' shape is not to sca

Figure 4 .
Figure 4. Cell envelope definition for filtering voxels with zero and non-zero volume.Re centers are considered to be inside the cell envelope, whereas blue-coloured voxel centers side the cell envelope based on the normal direction of the cell envelope (orange arrow).

Figure 4 .
Figure 4. Cell envelope definition for filtering voxels with zero and non-zero volume.Red voxel centers are considered to be inside the cell envelope, whereas blue-coloured voxel centers are outside the cell envelope based on the normal direction of the cell envelope (orange arrow).

Figure 7 .
Figure 7. Numerical error analysis of a closed hexagon with voxel aspect ratios of unity.(a) The large non-zero elastic constants.For a constant truss radius, the volume fraction changes due to the voxel meshing process, and the volume fraction converges for a smaller voxel mesh.(b) The homogenized effective properties of the hexagon lattice, where the fluctuations are caused by the voxel meshing process and approximating the periodic boundary conditions.(c) The residual terms in a log plot, where the residuals fluctuate between a high and a low value, corresponds to the error caused by approximating the periodic boundary conditions.(compare colored circles in (d,e)) The periodic boundary conditions for the voxelized hexagon lattice are marked in (c).(d) Large error (AR ≈ 1) for 46 by 39 divisions, VF = 0.263.This voxelized hexagon corresponds to the upper red arrow shown in (c).(e) Small error (AR ≈ 1) for 47 by 40 divisions, VF = 0.262.This voxelized hexagon corresponds to the lower red arrow shown in (c).

Figure 7 .
Figure 7. Numerical error analysis of a closed hexagon with voxel aspect ratios of unity.(a) The large non-zero elastic constants.For a constant truss radius, the volume fraction changes due to the voxel meshing process, and the volume fraction converges for a smaller voxel mesh.(b) The homogenized effective properties of the hexagon lattice, where the fluctuations are caused by the voxel meshing process and approximating the periodic boundary conditions.(c) The residual terms in a log plot, where the residuals fluctuate between a high and a low value, corresponds to the error caused by approximating the periodic boundary conditions.(compare colored circles in (d,e)) The periodic boundary conditions for the voxelized hexagon lattice are marked in (c).(d) Large error (AR ≈ 1) for 46 by 39 divisions, VF = 0.263.This voxelized hexagon corresponds to the upper red arrow shown in (c).(e) Small error (AR ≈ 1) for 47 by 40 divisions, VF = 0.262.This voxelized hexagon corresponds to the lower red arrow shown in (c).

Figure 8 .
Figure 8. Variation of volume fraction due to increasing the number of voxel divisions in the ydirection while maintaining an aspect ratio of unity.The voxel's center coordinate has been represented as smaller voxels.The periodic basis is plotted as a bold green line.

Figure 8 .
Figure 8. Variation of volume fraction due to increasing the number of voxel divisions in the y-direction while maintaining an aspect ratio of unity.The voxel's center coordinate has been represented as smaller voxels.The periodic basis is plotted as a bold green line.

Figure 9 .
Figure 9. Visualization of mismatched periodicity (black circles) and its effect on lattice symmetry.The periodic node pairs are connected using randomly colored lines.

Figure 9 .
Figure 9. Visualization of mismatched periodicity (black circles) and its effect on lattice symmetry.The periodic node pairs are connected using randomly colored lines.

Figure 10 .
Figure 10.Comparison of normalized  anisotropy plot at multiple volume fractions for (a) voxelized grid with ANSYS' cubic lattice model, and (b) voxelized hexagon with non-orthogonal periodic basis and ANSYS' Honeycomb lattice cell with orthogonal periodic basis.(c) Voxelized grid and ANSYS' cubic lattice model.(d) Voxelized hexagon with non-orthogonal periodic basis and ANSYS' Honeycomb model with orthogonal periodic basis

Figure 10 .
Figure 10.Comparison of normalized E 11 anisotropy plot at multiple volume fractions for (a) voxelized grid with ANSYS' cubic lattice model, and (b) voxelized hexagon with non-orthogonal periodic basis and ANSYS' Honeycomb lattice cell with orthogonal periodic basis.(c) Voxelized grid and ANSYS' cubic lattice model.(d) Voxelized hexagon with non-orthogonal periodic basis and ANSYS' Honeycomb model with orthogonal periodic basis.

Figure 11 .
Figure 11.Comparison of the [4,1] term of the homogenized elastic tensor for (a) hexagon lattice and (b) grid lattice with the [1,1] term for a material of  = 2× 10 and  = 0.3.The CH and D labels correspond to the data from the voxelization and ANSYS material modeller, respectively.

Figure 11 .
Figure 11.Comparison of the [1,4] term of the homogenized elastic tensor for (a) hexagon lattice and (b) grid lattice with the [1,1] term for a material of E iso = 2 × 10 11 and nu iso = 0.3.The CH and D labels correspond to the data from the voxelization and ANSYS material modeller, respectively.

Figure 12 .
Figure 12.Anisotropy plot of the Youngs Modulus Elastic property of a 3D primitive cubic lattice cell, where the red circles mark the  values at multiple rotation cell angles similar to the rotated cells shown in Figure 13.The dashed  values are overlapped with the  values due to geometrical symmetry.

Figure 13 .
Figure 13.A grid lattice cell discretized by 30 voxels along the x, y, and z axis, where the unit cell

Figure 12 .
Figure 12.Anisotropy plot of the Youngs Modulus Elastic property of a 3D primitive cubic lattice cell, where the red circles mark the E 11 values at multiple rotation cell angles similar to the rotated cells shown in Figure 13.The dashed E 11 values are overlapped with the E 22 values due to geometrical symmetry.

Figure 12 .
Figure 12.Anisotropy plot of the Youngs Modulus Elastic property of a 3D primitive cubic lattice cell, where the red circles mark the  values at multiple rotation cell angles similar to the rotated cells shown in Figure 13.The dashed  values are overlapped with the  values due to geometrical symmetry.

Figure 13 .
Figure 13.A grid lattice cell discretized by 30 voxels along the x, y, and z axis, where the unit cell rotated at 0°, 22°, 45°, and 68° along the z-axis.

Figure 13 .
Figure 13.A grid lattice cell discretized by 30 voxels along the x, y, and z axis, where the unit cell rotated at 0 • , 22 • , 45 • , and 68 • along the z-axis.

Figure 14 .
Figure 14.Visualization of the integral in Equation (18) for individual voxel elements, a hexagon with  / = 1.155 and a voxel shape of (a) 30 × 30 × 1 and (b) 150 × 150 × 1.Each of the individual hexagon visualizations were normalized by the value of the  term.

Figure 14 .
Figure 14.Visualization of the integral in Equation (18) for individual voxel elements, a hexagon with l x /l y = 1.155 and a voxel shape of (a) 30 × 30 × 1 and (b) 150 × 150 × 1.Each of the individual hexagon visualizations were normalized by the value of the C H ij term.

( a )
Approximated periodic boundary condition for a coarsely meshed hexagon lattice.(b) Approximated periodic boundary condition for a finely meshed hexagon lattice.

Figure 15 .
Figure15.Visualization of the periodicity node pairs.Each line of a random colour represents a coupling of the degrees of freedom for that node pair.The periodic basis has been highlighted as thick green lines.The approximated periodic boundary condition can be observed while comparing how the periodic boundary condition is applied at (a) between regions a1 to a5 and a2 to a4.Due to the finer mesh size, the approximation of the periodic boundary condition at (b) is not affected.

Figure 15 .
Figure15.Visualization of the periodicity node pairs.Each line of a random colour represents a coupling of the degrees of freedom for that node pair.The periodic basis has been highlighted as thick green lines.The approximated periodic boundary condition can be observed while comparing how the periodic boundary condition is applied at (a) between regions a 1 to a 5 and a 2 to a 4 .Due to the finer mesh size, the approximation of the periodic boundary condition at (b) is not affected.

Figure 16 .
Figure 16.Convergence of error terms for the hexagon lattice with non-unity aspect ratio.The datum shown above is a reduced subset of Figure 7c to emphasize the decreasing residual errors (thick red line).

Figure 17 .
Figure 17.Comparison of normalized  for the hexagon lattice at multiple volume fractions across multiple homogenization schemes present in the literature [26,36].

Figure 16 .
Figure 16.Convergence of error terms for the hexagon lattice with non-unity aspect ratio.The datum shown above is a reduced subset of Figure 7c to emphasize the decreasing residual errors (thick red line).

Materials 2023 , 33 Figure 16 .
Figure 16.Convergence of error terms for the hexagon lattice with non-unity aspect ratio.The datum shown above is a reduced subset of Figure 7c to emphasize the decreasing residual errors (thick red line).

Figure 17 .
Figure 17.Comparison of normalized  for the hexagon lattice at multiple volume fractions across multiple homogenization schemes present in the literature [26,36].

Figure 17 . 33 Figure 16 .
Figure 17.Comparison of normalized E 11 for the hexagon lattice at multiple volume fractions across multiple homogenization schemes present in the literature [26,36].

Figure 17 .
Figure 17.Comparison of normalized  for the hexagon lattice at multiple volume fractions across multiple homogenization schemes present in the literature [26,36].

Figure 18 .
Figure 18.Multiple representations of the Honeycomb hexagon lattice with different representative volume elements with their corresponding periodic basis.The representative volume element of the unit cell is shown as red voxels.

Figure 18 .
Figure 18.Multiple representations of the Honeycomb hexagon lattice with different representative volume elements with their corresponding periodic basis.The representative volume element of the unit cell is shown as red voxels.

Figure 19 .
Figure 19.Elastic properties of the 2D hexagon lattice with a relative density of 0.25 made from isotropic material with Young's modulus of 2 × 10 11 Pa and Poisson's ratio of 0.3.Two-dimensional anisotropic diagram of elastic properties for the hexagon with a volume fraction of 0.2 due to the variation of the hexagon cell angle.

Figure 19 .
Figure 19.Elastic properties of the 2D hexagon lattice with a relative density of 0.25 made from isotropic material with Young's modulus of 2 × 10 11 Pa and Poisson's ratio of 0.3.Two-dimensional anisotropic diagram of elastic properties for the hexagon with a volume fraction of 0.2 due to the variation of the hexagon cell angle.

Figure 19 .
Figure 19.Elastic properties of the 2D hexagon lattice with a relative density of 0.25 made from isotropic material with Young's modulus of 2 × 10 11 Pa and Poisson's ratio of 0.3.Two-dimensional anisotropic diagram of elastic properties for the hexagon with a volume fraction of 0.2 due to the variation of the hexagon cell angle.

Figure 20 .
Figure 20.Comparison of homogenized elastic properties of a square lattice with varying cell angles between 2D MATLAB voxel and 3D voxel code.

Figure 20 .
Figure 20.Comparison of homogenized elastic properties of a square lattice with varying cell angles between 2D MATLAB voxel and 3D voxel code.

Materials 2023 ,
16, x FOR PEER REVIEW 21 of 334.1.Three-Dimensional Triclinic and a Monoclinic Bravais Grid LatticeIn this section, a triclinic and a monoclinic Bravais grid lattice are analyzed by applying the approximated periodicity boundary condition.The geometry was discretized using voxels, as shown in Figure21.Visualization of the approximated periodicity.Voxel center (red), voxel nodes (grey), boundary condition (multi-coloured lines), and the cumulative periodic basis (thick green line).

Figure 21 .
Figure 21.Discretized triclinic Bravais grid lattice.The Voxel center is represented as a small sphere for a volume fraction of 0.3.The planes that represent the cell envelope definition and its normal have also been plotted.

Figure 21 .
Figure 21.Discretized triclinic Bravais grid lattice.The Voxel center is represented as a small sphere for a volume fraction of 0.3.The planes that represent the cell envelope definition and its normal have also been plotted.

Figure 23 .Table 3 .
Figure 23.Geometric definitions for the monoclinic and triclinic grid lattice are used in this paper.Table3.Elastic properties of the sample monoclinic and triclinic Bravais lattice.

Figure 23 .
Figure 23.Geometric definitions for the monoclinic and triclinic grid lattice are used in this paper.

Figure 24 .
Figure24.Two-dimensional and three-dimensional anisotropy plots for triclinic Bravais grid lattice (3D).Three-dimensional anisotropy plots are not to scale.Figure24.Two-dimensional and three-dimensional anisotropy plots for triclinic Bravais grid lattice (3D).Three-dimensional anisotropy plots are not to scale.

Figure 25 .
Figure 25.Two-dimensional and three-dimensional anisotropy plots for monoclinic Bravais grid lattice (3D).Three-dimensional anisotropy plots are not to scale.

Figure 26 .
Figure 26.3.4.6.4 lattice cells with a non-orthogonal periodic basis.The representative volume element of the unit cell is shown as red voxels.Evolution of selected 2D latices as the volume fraction increases.

Figure 26
Figure 26 depicts the 3.4.6.4 lattice for multiple volume fractions.This figure shows how the geometry of the lattice changes as the volume fraction increases.This sort of filling behaviour is one of the reasons for the Euler-Bernoulli beam formulation being invalid for larger relative density, as it cannot predict the interactions of the geometry within the lattice.The 2D and 3D anisotropy plots for the 3.4.6.4 lattice are shown in Figure 27.The coefficients for the homogenized elastic tensor as a function of relative density is tabulated in Table4.

Figure 26 .
Figure 26.3.4.6.4 lattice cells with a non-orthogonal periodic basis.The representative volume element of the unit cell is shown as red voxels.Evolution of selected 2D latices as the volume fraction increases.

Figure 26
Figure 26 depicts the 3.4.6.4 lattice for multiple volume fractions.This figure shows how the geometry of the lattice changes as the volume fraction increases.This sort of filling behaviour is one of the reasons for the Euler-Bernoulli beam formulation being invalid for larger relative density, as it cannot predict the interactions of the geometry within the lattice.The 2D and 3D anisotropy plots for the 3.4.6.4 lattice are shown in Figure 27.The coefficients for the homogenized elastic tensor as a function of relative density is tabulated in Table 4. Materials 2023, 16, x FOR PEER REVIEW 28 of 33

23 Figure 26
Figure26depicts the 3.4.6.4 lattice for multiple volume fractions.This figure shows how the geometry of the lattice changes as the volume fraction increases.This sort of filling behaviour is one of the reasons for the Euler-Bernoulli beam formulation being invalid for larger relative density, as it cannot predict the interactions of the geometry within the lattice.The 2D and 3D anisotropy plots for the 3.4.6.4 lattice are shown in Figure27.The coefficients for the homogenized elastic tensor as a function of relative density is tabulated in Table4.

Figure 28 .
Figure 28.Two-dimensional and three-dimensional anisotropy plots for X lattice in a sandwich panel (3D).

Table 1 .
(20)arison of eigenvalues of lower resolution(19)and higher resolution(20)hexagon lattice for an imperfectly voxelized hexagon with a voxel aspect ratio of 1.155.

Table 2 .
Angular constraints for the Bravais lattice are based on the angle definitions used in this paper.

Table 3 .
Elastic properties of the sample monoclinic and triclinic Bravais lattice.

Table 5 .
Homogenized elastic tensor coefficients (based on Equation (13)) and correlation coefficient ( ) for the polynomial fit for the Sandwich X lattice.

Table 5 .
Homogenized elastic tensor coefficients (based on Equation (13)) and correlation coefficient ( ) for the polynomial fit for the Sandwich X lattice.