## Abstract

Computational first-order homogenization theory is used for the elastic analysis of generally anisotropic lattice materials within classical continuum mechanics. The computational model is tailored for structural one-dimensional (1D) elements, which considerably reduces the computational cost comparing to previously developed models based on solid elements. The effective elastic behavior of lattice materials is derived consistently with several homogenization approaches including strain- and stress-based methods together with volume and surface averaging. Comparing the homogenization based on the Hill–Mandel Lemma and constitutive approach, a shear correction factor is also introduced. In contrast to prior studies that are usually limited to a specific class of lattice materials such as lattices with cubic symmetry or similarly situated joints, this computational tool is applicable for the analysis of any planar or spatial stretching- and bending-dominated lattices with arbitrary topology and anisotropy. Having derived the elasticity of the lattice, the homogenization is then complemented by the symmetry identification based on the monoclinic distance function. This step is essential for lattices with non-apparent symmetry. Using the computational model, nine different spatial anisotropic lattices are studied among which four are fully characterized for the first time, i.e., non-regular tetrahedron (with trigonal symmetry), rhombicuboctahedron type a (with cubic symmetry), rhombicuboctahedron type b (with transverse isotropy), and double-pyramid dodecahedron (with tetragonal symmetry).

## 1 Introduction

Lattice materials are periodic space-filling structures formed by tessellating an open unit-cell in two- or three-dimensional (3D) space. The unit-cell is the building block of such materials and is made of meso/micro scale struts. The mechanical properties of lattice materials highly depend on their unit-cell architecture. Depending on the topology of the unit-cell, the main deformation mode is dominated by either stretching or bending of the struts. On such basis, these structures are classified as bending-dominated lattice materials (BDLM) and stretching-dominated lattice materials (SDLM).

Lattice materials have great weight-saving potential and are exceptional candidates for multifunctional applications, e.g., impact energy absorption, load carrying, acoustic damping, enabling fluid flow and heat exchanging [1]. Previously, the manufacturing of lattice materials was challenging and often expensive, especially for complicated architectures. However, with the aid of additive manufacturing, it is now possible to build lattices with almost any arbitrary topology, and thus, they have gained great attention in a wide range of industries, e.g., see Ref. [2].

The earlier works in the literature focused on the understanding of the behavior of such materials. Ashby [3] studied the mechanical behavior of cellular structures and introduced their performance graphs. Following that, extensive studies addressed the derivation of effective properties of lattice materials. According to the literature, there exist two main approaches toward the derivation of effective properties, i.e., analytical approaches and homogenization methods.

In the analytical approach, a lattice unit-cell is treated as a truss/frame structure and is analyzed according to the classical beam theories using Neumann (uniform traction) boundary condition. Nonetheless, the analytical approach does not have a unified framework for different lattices [4]. In other words, the implementation of this approach differs for each topology and is not straightforward when it comes to complex geometries with arbitrary anisotropy. In fact, most of the analytical studies are focused on only one lattice architecture and not a family of lattices, e.g., octet lattice in Ref. [5], honeycombs in Ref. [6], Kelvin in Ref. [7], rhombicuboctahedron in Ref. [8], dodecahedron in Ref. [9], and double-pyramid dodecahedron in Ref. [10]. To generalize the existing analytical models, Abdelhamid and Czekanski [11] extended the previously developed expressions to include the elastic properties of octet truss with various lattice angles, but their study is also limited to octet topology.

In addition to topology-specific cases, some studies presented more generic analytical methods and targeted special classes of lattice materials. Gurtner and Durand [12] analytically studied stiff elastic networks (SDLM) and Norris [13] developed a method to derive the effective properties of lattice materials with similarly situated central nodes. However, still many lattice materials do not belong to this category. Moreover, the majority of the literature is devoted to the geometries with cubic symmetry, e.g., Refs. [5,7,8] and Ref. [13], while those with higher levels of anisotropy have not been fully characterized. For example, the shear properties of dodecahedron were not reported in Ref. [9], or both shear and Poisson’s ratio were excluded for double-pyramid dodecahedron in Ref. [10], and only the elastic modulus in certain directions was reported. Even in studies aiming at lattices with higher levels of anisotropy, e.g., Refs. [9] and [10], the unit-cells possess a few apparent symmetry planes, which made their analysis possible for certain loadings and/or in certain directions. However, the derivation of effective properties with the aid of analytical methods for lattice materials with arbitrary anisotropy (with no apparent symmetry plane) is a very daunting task.

To overcome the challenges and narrow applicability range of analytical methods, more recent studies tend to employ homogenization theories for the analysis of lattice materials. In homogenization, a set of mathematical theories are implemented to replace a discrete unit-cell with its equivalent homogenous medium. Yvonnet [14] employed a computational homogenization model based on averaging theory and solid elements and kinematically uniform boundary conditions to obtain the effective properties of diamond and octet lattices. Zhu et al. [15] used a different computational homogenization model, based on variational asymptotic method, with solid elements and obtained the elastic properties of Kelvin cell. Vigliotti and Pasini [16] developed a multiscale computational method to study the stiffness and strength properties of tridimensional lattices. Nevertheless, the analysis in Ref. [16] is limited to lattices with cubic symmetry, and higher anisotropies were excluded. In another work, Arabnejad and Pasini [17] utilized asymptotic homogenization to derive the effective properties of a number of two-dimensional (2D) lattices with specific symmetries. Dong et al. [18] presented a computational scheme (based on asymptotic homogenization and solid elements) to assess the elasticity of multi-material lattices with various 3D topologies. Karathanasopoulos et al. [19] have recently developed a discrete mechanics code for elastic analysis of so-called 2D metamaterials. Tancogne-Dejean et al. [20] established structure-elastic property relations for 2D hexachiral lattices using finite element simulations.

Similar to averaging and asymptotic homogenization, elastodynamic homogenization based on Bloch-wave theorem is an effective tool for the homogenization of lattice materials. Hutchinson and Fleck [21] applied the Bloch-wave theorem to analyze 2D Kagome and triangular–triangular lattices. Later, Elsayed and Pasini [22] utilized the same theory for the analysis of 2D stretching-dominated lattices with square and /or hexagonal symmetry. More recently, Patil and Matlack [23] implemented the Bloch-wave homogenization to derive the static properties of several 3D lattices with different symmetries. Bloch-wave theorem suits for addressing simultaneously the statistics and dynamics of lattice materials. This increases the complexity of the problem and leads to higher computational cost in comparison to static homogenization. Accordingly, we employ the elastostatic homogenization in this study to address generally anisotropic lattice materials.

The discrete lattice, once possessing the necessary scale separation, can be homogenized toward a continuum. Depending on the scale separation, first-order or higher-order homogenization would lead to a classical continuum or generalized continua, respectively. In this study, the first-order homogenization is employed to address a generally anisotropic lattice within classical continuum mechanics. For studies on generalized continua, see for example, Refs. [24–28] and Ref. [29].

The topic of elastic analysis of lattice materials has been studied extensively, but the development of an efficient generic model, capable of addressing complex 3D lattices with an arbitrary level of anisotropy, is still an ongoing research, see for example Refs. [14,15,23,30,31] and Ref. [32]. The current study intends to contribute to the field by developing a computational first-order homogenization model for the elastostatic analysis of low-density and generally anisotropic lattice materials. In comparison with the extensive literature on lattice materials, the novel features of this study include the following:

*Structural one-dimensional (1D) elements for 3D elastic analysis*: Previously developed computational models are based on solid elements, and this often leads to considerable computational time. For instance, Zhu et al. [15] used 245,006 solid elements for the elastic analysis of Kelvin cell with 6% relative density, while, at the given density, this can be done with only 24 beam elements (see Sec. 5.6). Indeed, lattice materials, as a group of lightweight materials, are often built with low densities, and a customized model based on 1D elements can drastically reduce the computational cost. However, the analysis of tridimensional periodic structures with 1D elements brings several challenges, which will be treated in this paper.*General anisotropy*: In contrast to many previous studies, the proposed method in this study has no limitation on the geometrical complexity and anisotropy level of the unit-cell. Additionally, it addresses both SDLM and BDLM. For lattices with non-apparent symmetry, the homogenization is complemented with an anisotropy identification algorithm based on the elasticity distance function.*Characterization of new lattice materials*: Using the proposed computational model, several lattice materials are fully characterized for the first time including non-regular tetrahedron, rhombicuboctahedron type a, rhombicuboctahedron type b, and double-pyramid dodecahedron possessing, respectively, trigonal symmetry, cubic symmetry, transverse isotropy, and tetragonal symmetry.

The paper is organized as follows. In Sec. 2, a general homogenization model is developed for 3D elastic analysis of lattices, i.e., both SDLM and BDLM, using structural 1D elements. The practical implementation through the finite element method with periodic boundary conditions is also described. In Sec. 3, a simplified efficient model based on bar elements is suggested for the analysis of SDLM. This is complemented by the introduction of a new criterion for the determinacy analysis of periodic lattice materials to distinguish SDLM from BDLM. Following that, in Sec. 4, an extension to the simplified model is presented including stress-based homogenization and strain-based homogenization using surface averaging. In Sec. 5, the general homogenization model is employed to characterize the elastic properties of nine distinct 3D lattice materials with different symmetries including octet, tetrakaidekahedron or Gurtner–Durand [12], non-regular tetrahedron, cube, diamond, rhombicuboctahedron type a, rhombicuboctahedron type b, Kelvin and double-pyramid dodecahedron. The directional behavior of the lattices are discussed in Sec. 6. In Sec. 7, the simplified homogenization model is applied for five 2D and three 3D SDLM. Finally, the study is concluded in Sec. 8.

## 2 Homogenization of Lattices

The computational model in this study is developed based on a homogenization technique known as averaging method, e.g., see Ref. [14], and is implemented using the finite element method with structural 1D (generalized beam) elements. Homogenization can be implemented by applying either a uniform macroscopic strain field (strain-based) or a uniform macroscopic stress field (stress-based). Additionally, the averaging itself can be done both on the volume and on the surface of the representative volume element (RVE). Therefore, several computational models have been developed to include these approaches and compare their results.

Defining a proper RVE is the very first step in the homogenization of a heterogenous medium toward a specific continuum framework. Considering the phenomena to be addressed (such as elasticity, wave propagation, and buckling), an RVE should effectively represent all the microstructural features of the heterogeneous material. For the elastostatic analysis of anisotropic lattice materials, one should be able to build the infinite macroscale lattice structure by repeating the RVE along its periodicity directions (here, it is assumed to be along the perpendicular coordinate axes). Therefore, the RVE for the elastostatic homogenization of a lattice material might be different from its unit-cell.

In this section, a computational scheme is developed for the elastic analysis of generally anisotropic lattice materials, i.e., including both SDLM and BDLM. For this purpose, the effective elasticity is obtained based on the method of homogenization. Afterward, the implementation of the finite element methodology is elaborated.

### 2.1 Effective Elasticity.

*Ω*of the space, subjected to a uniform macroscopic strain field. For the elastostatic analysis of the material, the six macroscopic strain states are prescribed as follows [14]:

*e*

_{i}(

*i*= 1, 2, 3) are the unit basis vectors.

*ɛ*can be related to an arbitrary macroscopic strain field $\epsilon \xaf$, using a localization matrix [

*M*]. This can be written in the matrix form (Voigt notation) as

Here [*B*] is the derivative of the shape function (also referred to as gradient matrix) and [*d*^{e}] is the element nodal deformation matrix in local coordinate system, corresponding to the six possible strain states, see Eq. (1). The matrix [*d*^{e}] can be obtained using finite element method (FEM) and periodic boundary conditions (PBC), in which the lattice RVE is discretized into general beam elements. In the absence of body forces, each lattice strut is simulated with a two-node general beam element with rigid joints. Moreover, the formulation of the element bending is based on Euler–Bernoulli (E–B) beam theory.

*u*

_{i},

*v*

_{i}, and

*w*

_{i}are the nodal displacements along element’s local

*x*,

*y*, and

*z*axes, and

*φ*

_{ix},

*φ*

_{iy}, and

*φ*

_{iz}are the nodal rotations around element’s local

*x*,

*y*, and

*z*axes. Here

*i*is the node number, i.e.,

*i*= 1, 2, and

*x*is the axial beam direction. The matrix of nodal deformation of a two-node general beam element is given by

*d*] matrix should be calculated for all six macroscopic strain states and stored as the column matrices to form [

*d*

^{e}]. Hence, [

*d*

^{e}] is a 12 × 6 matrix, in which each column corresponds to one macroscopic strain state. That is

*d*]

^{(ij)}is the matrix of nodal deformation corresponding to macroscopic state $\epsilon \xaf(ij)$, see relation (1). Since the nodal deformations obtained by FEM are expressed in the global coordinate system, i.e., [

*D*

^{e}], the nodal deformation in the local coordinate system reads

*T*

_{beam}] is the transformation matrix for the element, e.g., see Ref. [33]. The general beam element is a combination of bar, shaft, and beam elements for axial, torsional, and bending analyses, respectively. Thus, its gradient matrix [

*B*] is composed of three different gradient matrices to capture all three modes of deformation including stretching, bending, and twisting (torsion), i.e.,

*B*

_{shaft}], which is commonly expressed in cylindrical coordinate system, is transformed into the Cartesian coordinate system (see Appendix A). The formulation is given for E–B theory for struts with circular cross section of radius

*r*. In this case, the analytical solution of torsion as well as symmetric bending simplify the formulation. If a different cross section or bending theory is of interest, the formulation should be adopted accordingly.

*σ*] is the microscopic stress vector of the element and [

*C*

^{e}] is the element elasticity matrix defined in its local coordinate system. For a general beam element, all components of the [

*C*

^{e}] are zero except the (1,1) which is equal the elastic modulus (

*E*

_{s}) as well as (5,5) and (6,6) components which are equal to the shear modulus (

*G*

_{s}) of the solid material.

*γ*

_{zx}=

*γ*

_{xy}= 0). Yet, in this method, the effect of shear forces (reflected in beam bending) is captured through bending deformation.

*HML*denotes that the elasticity matrix is derived using Hill–Mandel Lemma. Since the lattice unit-cell is composed of discrete elements, the above integral can be written in the following discretized form

*V*

_{RVE}is the volume of the RVE and

*N*

_{e}is the number of elements.

It is noted that the integral bounds in Eq. (15) are written for a circular cross section, i.e., 0 ≤ *x* ≤ *L*, $\u2212r2\u2212z2\u2264y\u2264+r2\u2212z2$ and −*r* ≤ *z* ≤ + *r*. These bounds are written in general form, i.e., for the complete struts (internal members). However, if the struts are located on the boundary of the RVE, the integral bounds on *y* and *z* are halved or quartered, depending on the geometry.

*τ*

_{s}], along with a shear correction factor,

*k*

_{s}. The element’s modified stress vector reads as

*y*and

*z*axes corresponding to macroscopic state $\epsilon \xaf(ij)$, and

*A*is the area of the element cross section.

The superscript *Const* indicates that the above elasticity matrix is obtained using the constitutive approach. Moreover, to achieve consistent results with the energy approach, the shear correction factor (*k*_{s}) turned out to be *k*_{s} = 2.

The stress vector in Eq. (16) was obtained in element’s local coordinate system. Thus, for each element, it should be transformed back into the global coordinate system using the stress transformation matrix [*T*_{s}] (see Appendix B). Such transformation is not needed in Eq. (15), where the energy is used for homogenization.

### 2.2 Finite Element Analysis.

*A*and

*B*being on two opposite sides of RVE, the periodicity condition reads [14]

*u*is the displacement field and

*X*is the position of the nodes. For a lattice unit-cell composed of general beam elements, with both translational (displacement) and rotational degrees-of-freedom (DOF), the periodicity condition can be formulated in the matrix form as

*D*] is the vector of unit-cell nodal deformations and [

*P*] is a matrix relating the deformation of paired nodes to the whole set of nodal deformations. For each set of paired nodes, e.g.,

*A*and

*B*, [

*R*] is a column matrix with six rows. Its first three components satisfy the condition on periodicity of the displacement field, obtained by Eq. (21), while its last three components are zero to enforce the continuity of the rotations. The latter condition is based on the classical continuum mechanics theory in which the martial particles do not have rotational DOF.

Based on the concept of minimum potential energy, constrained minimization method using Lagrange multipliers is implemented to find the unit-cell nodal displacement vector [*D*]. The [*D*] vector can be then used in Eq. (6) to obtain [*d*^{e}] to be inserted in Eqs. (15) and (20) to derive the effective elasticity matrix of the material.

*K*] is the unit-cell’s global stiffness matrix, and the matrix [

*Q*] removes the unit-cell rigid-body motion by constraining translational degrees of a chosen node. It is noted that the rotational degrees of the unit-cell are automatically constrained by PBC. The Lagrange multipliers [Λ

_{1}] and [Λ

_{2}] are used for periodicity and rigid-body-translation constraints, respectively. Also,

*N*

_{n}is the number of nodes and

*N*

_{p}is the number of paired sets. Upon solving the system of Eqs. (23), the unit-cell nodal deformation matrix (in the global coordinate) is found.

The Lagrange multipliers give rise to these external forces on the boundary nodes of the unit-cell to enforce the periodicity and rigid-body-translation constraints. This information is required once the homogenization is carried out by surface averaging (see Sec. 4.2.).

The above-mentioned process should be repeated for six different strain states to obtain their corresponding nodal force and deformation matrices.

## 3 Homogenization of SDLM

The introduced computational model presented in Sec. 2 is based on the rigid-node assumption to consider all possible DOF and include the analysis of lattice materials in general, i.e., both BDLM and SDLM. However, the model can be simplified for SDLM by using the pin-joint assumption and excluding nodal rotations, moments, and torsions. Since in SDLM the members deform only axially, they can be simulated as a pin-joint network. The simplified model can considerably reduce the computational time when it comes to the analysis of SDLM. The procedure is similar to what was described in Sec. 2.1, with some differences which are discussed here.

*d*], reduces to the nodal displacement matrix, [

*u*], defined as $[u]=[u1v1w1u2v2w2]T$. Similar to [

*d*

^{e}] in Eq. (5), the nodal displacement matrix is calculated for the six macroscopic strain states and is stored as the column vectors to form a 6 × 6 matrix [

*u*

^{e}]. Additionally, the gradient matrix described in Eq. (7) reduces to [

*B*] = [

*B*

_{bar}], and the element’s strain vector reads

*u*

^{e}] = [

*T*

_{bar}][

*U*

^{e}]. Here, [

*U*

^{e}] is the nodal displacement matrix in the global coordinate system obtained by FEM, and [

*T*

_{bar}] is the transformation matrix for the bar element [33].

*C*

^{e}] is equal to the elastic modulus (

*E*

_{s}) of the solid material. Finally, following the same procedure described in Sec. 2.1, the elasticity matrix based on the constitutive approach reads

*B*

_{bar}] matrix, Eqs. (26) and (27) can be written in the form of a summation series as

*V*

^{e}is the volume of element.

Due to the simplicity and computational efficiency, the simplified model is preferred for the analysis of SDLM. However, this requires a criterion to distinguish SDLM from BDLM. The discussion of the determinacy involves both the unit-cell (finite structure) and the lattice material (infinite structure). Considering the differences of finite and infinite structures, a criterion will be introduced in the following based on the computational homogenization scheme.

### 3.1 Determinacy Analysis.

Consider the lattice unit-cell as a finite pin-joint truss structure. For *kinematically indeterminate structures*, under an external load, the structure may collapse by the rotation of its struts around joints creating neither stress nor strain. Such structures act as *mechanism*. On the other hand, the structure may support external loads by extension or compression of its members without collapsing. In this case, if there is no redundant member (a member whose removal does not make the structure a mechanism), the structure is said to be *statically determinate*. However, if the structure possesses redundant members, it is *statically indeterminate* and have *states of self-stress* equal to the number of redundant members.

To assess the rigidity of a finite structure, the generalized Maxwell rule was introduced as *M* = *s* − *m*, where *s* and *m* and are number of self-stress states and number of mechanisms, respectively [34]. The structure is statically and kinematically determinate, if both *s* and *m* are zero. The values for *s* and *m* can be determined using the equilibrium equation [*A*_{eq}][*t*_{strut}] = [*f*], where [*A*_{eq}] is the equilibrium matrix containing the cosine vector of the members with respect to the global coordinate system, [*t*_{strut}] is the matrix of strut tensions, and *f* is the matrix of nodal forces. The constrained nodes are not included when forming this equation [34]. Using the rank of equilibrium matrix ($rAeq$), the value for *s* can be obtained by $s=b\u2212rAeq$, where *b* is the number of struts. Finally, *m* can be found using the generalized Maxwell’s rule. In the literature, this approach is often referred to as determinacy analysis.

Here, for a range of 2D (see Fig. 1) and 3D unit-cell topologies (see Fig. 2), the determinacy analysis is carried out using an in-house matlab code. In addition to the equilibrium matrix, [*A*_{eq}], the stiffness matrix governed by classical FEM, i.e., [*K*], is also used for the determinacy analysis. It turned out that both equilibrium and stiffness matrices provide the same results, i.e., $rAeq=rk$. This can be simply proved by rewriting the equilibrium equation in terms of nodal displacements.

The results for the determinacy analysis of the unit-cells are presented in Tables 1 and 2. These tables also include the behavior of the periodic lattices (i.e., stretch- or bending-dominated). It is noted that the analysis is performed on the unit-cell basis, and consequently, *s* and *m* values belong to the lattice unit-cell and not the periodic lattice material. Besides, the selection of lattice unit-cell would affect the rank of equilibrium and stiffness matrices (and consequently the values for *m* and *s*). The determinacy of a unit-cell is a function of the number of nodes, numbers of struts, and spatial arrangement of struts (topology). Therefore, the reported result for *m* and *s* may change if the unit-cell is selected differently.

Topology | s | m | M | Periodic behavior | |
---|---|---|---|---|---|

1 | Triangular | 0 | 0 | 0 | SD |

2 | Kagome | 0 | 3 | −3 | SD |

3 | Mix-square type a | 1 | 0 | 1 | SD |

4 | Mix-square type b | 1 | 0 | 1 | SD |

5 | Mix-square type c | 0 | 0 | 0 | SD |

6 | Diamond | 0 | 0 | 0 | SD |

7 | Square web | 2 | 1 | 1 | BD |

8 | Square | 0 | 1 | −1 | BD |

9 | Hexagonal | 0 | 3 | −3 | BD |

Topology | s | m | M | Periodic behavior | |
---|---|---|---|---|---|

1 | Triangular | 0 | 0 | 0 | SD |

2 | Kagome | 0 | 3 | −3 | SD |

3 | Mix-square type a | 1 | 0 | 1 | SD |

4 | Mix-square type b | 1 | 0 | 1 | SD |

5 | Mix-square type c | 0 | 0 | 0 | SD |

6 | Diamond | 0 | 0 | 0 | SD |

7 | Square web | 2 | 1 | 1 | BD |

8 | Square | 0 | 1 | −1 | BD |

9 | Hexagonal | 0 | 3 | −3 | BD |

Topology | s | m | M | Periodic behavior | |
---|---|---|---|---|---|

1 | Octet | 0 | 0 | 0 | SD |

2 | Tetrakaidekahedral | 1 | 14 | −13 | SD |

3 | Tetrahedron | 5 | 0 | 5 | SD |

4 | Cube | 0 | 6 | −6 | BD |

Topology | s | m | M | Periodic behavior | |
---|---|---|---|---|---|

1 | Octet | 0 | 0 | 0 | SD |

2 | Tetrakaidekahedral | 1 | 14 | −13 | SD |

3 | Tetrahedron | 5 | 0 | 5 | SD |

4 | Cube | 0 | 6 | −6 | BD |

Interpreting Tables 1 and 2 highlights that being statically and kinetically determinate, i.e., *s* = *m* = 0, at the unit-cell level guarantees the stretching-dominated behavior of the periodic structure but the violation of this criterion does not necessarily mean the periodic structure is not stretching-dominated. In fact, the results of determinacy analysis at the unit-cell level may not be informative enough for the assessment of the periodic structure, see also Refs. [21] and [35]. To overcome this issue, a criterion (based on Bloch’s theorem) is suggested by Hutchinson and Fleck for determinacy analysis of periodic structures [21].

In this study, using the concept of computational homogenization, a new criterion is derived for the analysis of deformation behavior of periodic networks. The proposed method involves less complexity (in comparison with Bloch’s theorem [21] for example) and has no limitations on the geometry and can accurately distinguish stretching-dominant (SD) from bending-dominant (BD) periodic networks. The last columns in Tables 1 and 2 are reported based on this criterion.

*D*], with the matrix of nodal displacements, [

*U*] (described in Sec. 3), and adjusting the dimensions of other matrices accordingly. This results in the following system of equations

*T*]. Equation (30) is solvable ([

*T*] is non-singular) for SDLM, while [

*T*] is singular for BDLM. This provides a novel criterion for the identification and classification of lattice materials. A lattice is stretching-dominated if its periodicity equilibrium matrix [

*T*] is of a full-rank, i.e.,

*r*

_{(T)}is the rank of periodicity equilibrium matrix [

*T*],

*d*is the dimension of the problem, and

*c*is the number of motions constrains. A lattice material is bending-dominated, if it does not satisfy the above equation. Generally, due to the removal of rigid-body rotations by PBC, the value of

*c*to prevent rigid-body motion is 2 and 3 for 2- and 3D problems, respectively. However, in special cases where RVE faces cut some struts,

*c*may take some other values (more than 2 or 3) to treat the virtual nodes which are introduced in Sec. 7.

## 4 An Extension to the Simplified Homogenization

In Sec. 3, the homogenization scheme using strain-based volume averaging is elaborated. Nevertheless, homogenization via averaging can also be done by applying macroscopic stress fields (stress-based) and/or averaging over the RVE boundaries (surface averaging). In this section, these methods are presented as extensions to the computational scheme.

### 4.1 Stress Approach.

*σ*can be related to an arbitrary macroscopic strain field $\sigma \xaf$, using a localization matrix [

*N*] (see Appendix C for different macroscopic stress states). This can be written in the matrix form as [14]

*n*is the outward unit normal vector of the RVE surface,

*u*is the nodal displacement, and Γ is the RVE boundary. The nodal displacements of the unit-cell are again obtained by FEM. The finite element problem solved for six different macroscopic stress fields results in six different nodal displacements. Using Eq. (37), the six sets of nodal displacements yield six different macroscopic strain matrices, $\epsilon \xaf$. These matrices are written in vector form and stored as the columns of $[\epsilon \xaftotal]$ as

### 4.2 Strain Approach With Surface Averaging.

*t*

_{(x)}is the traction vector acting on the RVE faces and

*X*is the position of the boundary nodes. The external forces derived in FEM (see Eq. (24)) is used to determine these traction vectors. The macroscopic stress $\sigma \xaf$ is given in Eq. (41) in the tensor form, which should be converted into the vector form according to the Voigt notation for current computational scheme. By repeating the calculations for the six macroscopic strain states and storing the results as column vectors, the matrix $[\sigma \xaftotal]$ can be formed. The macroscopic strain matrices are also written in vector form to give $[\epsilon \xaftotal]$. Then, the effective elasticity matrix can be obtained using Eq. (40).

### 4.3 Affine Deformation in SDLM.

The surface averaging in Eqs. (37) and (41) can be simplified for the SDLM. Once a homogenous medium is loaded uniformly, it undergoes a homogeneous deformation [36]. Having an affine deformation field implies that the deformation gradient is position independent and constant at a given instance of time, i.e., $F=const$, leading to a constant strain field, i.e., $E=(FTF\u2212I)/2=const$. This ultimately results in a linear displacement field (called affine displacement field), i.e., $U=E\u22c5X$, in which *X* is the position of the material points.

*A*is the area of the RVE faces (identical for cubic RVE),

*n*

_{i}is the outward unit normal vector of the face

*i*,

*u*

_{ci}is the displacement of the centroid of the face

*i*,

*F*

_{i}is the external uniformly distributed force acting on the face

*i*, and $Xci$ is the position of the of the centroid of face

*i*.

## 5 Elasticity of Selected Anisotropic Lattices

An in-house matlb code was developed to implement the computational homogenization scheme described in Sec. 2. The code is used to analyze nine different 3D lattices with different levels of elastic anisotropy. The selected 3D topologies include octet (cubic symmetry), tetrakaidekahedron or Gurtner–Durand (isotropic), non-regular tetrahedron (trigonal symmetry), cube (cubic symmetry), diamond (cubic symmetry), rhombicuboctahedron type a (cubic symmetry), rhombicuboctahedron type b (tetragonal symmetry), Kelvin (cubic symmetry), and double-pyramid dodecahedron (transverse isotropy). The RVEs of the above-mentioned lattice materials are shown in Fig. 2.

From the list of candidate materials, four are being fully characterized for the first time, i.e., non-regular tetrahedron, rhombicuboctahedron type a, rhombicuboctahedron type b, and double-pyramid dodecahedron. The model is applicable to any arbitrary low-density lattice material with any level of anisotropy and symmetry. The computational model is fully automated and determines the effective stiffness and compliance matrix of lattice materials (both SDLM and BDLM), from which the engineering constants and the material symmetry class can be obtained. For each topology, the effective elasticity matrix is obtained using HML method.

The main assumption here is that the unit-cell members can be simulated with E–B beam elements. This requires the members to be slender, i.e., *r*/*l* ≤ 0.1, and will ultimately lead to low-density lattice materials (typically $\rho \xaf\u22640.2$). The model is developed for lattices made of metallic struts with circular cross section but can be adopted for other shapes and/or materials as well.

### 5.1 Regular Octet.

*l*) and cross-section area (

*A*), the relative density of the lattice reads

*E*

_{s}is the elastic modulus of the base material. The obtained results agree with literature, i.e., Refs. [5,13] and Ref. [16]. It is noted that the compliance matrix derived by Deshpande et al. [5] is to be interpreted in the Voigt notation, since they did the analysis based on the engineering shear stain, i.e.,

*γ*

_{ij}= 2

*ɛ*

_{ij}(see Ref. [16] for more details).

### 5.2 Tetrakaidekahedron.

Tetrakaidekahedron lattice is an SDLM which is formed by connecting the center to all faces of the Kelvin cell, as described by Refs. [12] and [13]. This description is more useful when the lattice is to be analyzed from node arrangement perspective, i.e., similarly situated node. Alternatively, this lattice can be constructed using the unit-cell shown in Fig. 1. The latter is used here as the RVE for homogenization. The fabrication of this unit-cell using two different set of struts (with $A2=((33)/4)A1$) yields an optimized isotropic material, known as the stiffest isotropic network, as suggested by Gurtner and Durand [12]. The struts with *A*_{1} and *A*_{2} cross-sectional areas are shown with blue and red colors, respectively, in Fig. 2. This special class of tetrakaidekahedron lattice is referred to as Gurtner–Durand (G–D) lattice in the literature.

*l*) as

### 5.3 Irregular Tetrahedron.

The developed computational model is applicable to a generally anisotropic lattice material. To highlight this feature, tetrahedron lattice is homogenized here. Tetrahedron is an SDLM with interesting structural features but is not space-filling as a single unit-cell. However, it is possible to spatially arrange multiple tetrahedrons and form a space-filling unit-cell (or RVE). Goldberg [37] has introduced different classes of tetrahedral space-fillers. In this section, a specific case of these space-filling tetrahedrons, called Sommerville #3 [37], is studied. The unit-cell is composed of 12 non-regular tetrahedrons as shown in Fig. 2, where the internal members are colored in red and the boundary ones are colored as blue. This lattice is a promising candidate for biomedical applications [38], but its elastostatic behavior has not been reported possibly due to its relatively high level of anisotropy.

*A*), the relative density of the lattice can be characterized by the length of the unit-cell edge (

*l*) as

This function is illustrated on a sphere and in a contour plot in Figs. 3 and 4, respectively. Angles *θ*, *ϕ*, and *ψ* are the intrinsic Euler angles, being the rotation angles around *z*, *x*, and *y* axes, respectively. Setting *θ* = 0, the rotation angles *ϕ* and *ψ* can be used to identify the symmetry class as well as the material principal direction based on the monoclinic distance function.

*z*axis is selected to be normal to the plane of the coplanar normal vectors, and

*x*- and

*y*- axes are selected to have $45deg$ angle with respect to, arbitrarily, one of the three normal vectors. The corresponding stiffness matrix in the natural coordinate system, in which the material symmetry class is apparent, is found to be

*η*

_{i,jk}and

*μ*

_{ij,kl}are the effective interaction and Chentsov coefficients, respectively. From the above matrix, the effective engineering constants are obtained as

*ν*

_{13}/

*E*

_{1}=

*ν*

_{31}/

*E*

_{3},

*ν*

_{23}/

*E*

_{2}=

*ν*

_{32}/

*E*

_{3},

*G*

_{12}= 1/2((1/

*E*

_{1}) + (

*ν*

_{21}/

*E*

_{2}))

^{−1}, and

*μ*

_{23,12}/

*G*

_{12}= −2(

*η*

_{1,13}/

*G*

_{13}). To the authors’ knowledge, this lattice is characterized for the first time in the literature.

### 5.4 Regular Cube.

*l*) and cross-section area (

*A*), the relative density of the lattice is given by

As was observed for octet, G–D and irregular tetrahedron, the effective elastic properties of SDLM have the same order of magnitude and they all scale with $\rho \xaf$. However, the elastic properties of BDLM may scale dissimilarly under various macroscopic strains states $\epsilon \xaf(ij)$, depending on whether the bending, stretching or a mixture of both modes are activated. If the applied macroscopic strain only activates the bending mode, the corresponding elastic properties scale with $\rho \xaf2$.

### 5.5 Diamond.

For very low densities, i.e., $\rho \xaf\u22640.1$, the above expressions can be simplified as $E\u2243(1/3)Es\rho \xaf2$, *ν* ≃ 0.5 and $G\u2243(1/3)Es\rho \xaf2$. The smaller the relative density is, the more accurate are these relations. It is noted that the elasticity matrix derived in computational model, matched with the results of Ref. [13].

### 5.6 Kelvin.

*l*) and cross-section area (

*A*), the relative density of the lattice is

The Kelvin lattice has mixed mode deformation in all strain states and is categorized as a BDLM. Additionally, it has an active twisting mode in shear states. The twisting mode does not affect the elastic modulus, but it has a profound impact on shear modulus. In fact, the studies which have ignored the beam twisting mode, e.g., Ref. [42], cannot properly predict the shear modulus. The reason is that the infinite torsional rigidity assumption in Kelvin cell (as opposed to cube or diamond lattice) leads to an over-estimated shear stiffness.

*c*

_{f}is a torsional correction factor which is equal to 1.52 for lattices made of metallic struts with circular cross section. According to Ref. [40], this factor can generally be obtained for any cross-section shape and material using the expression

*c*

_{f}= (8

*E*

_{s}

*I*+

*G*

_{s}

*J*)/(5

*E*

_{s}

*I*+

*G*

_{s}

*J*), where

*G*

_{s}is the shear modulus of the solid material,

*I*is the second moment of area, and

*J*is the polar moment of area of the beam cross section. For very low densities, $\rho \xaf\u22640.1$, the above elastic properties can be approximated as $E\u22430.6Es\rho \xaf2$,

*ν*≃ 0.5, and $G\u22430.6Es\rho \xaf2/cf$. The obtained elastic properties in this computational study agreed with those reported in Refs. [7] and [40].

### 5.7 Rhombicuboctahedron (Type a).

*l*and

*A*are the length and cross-sectional area of the struts. The rhombicuboctahedron unit-cell can form several lattices with different properties, depending on its tessellation directions. In the current study, two different possible configurations of this unit-cell are investigated. These lattices are differentiated by adding the extensions “type a” and “type b”, see Fig. 5. Both lattices (type a) and (type b) share the same unit-cell but their periodicity directions, and consequently, their spatial arrangements are different.

The current section focuses on (type a), which often only referred to as rhombicuboctahedron. This lattice has not been fully characterized in the literature. The elastic modulus and Poisson’s ratio of the perfect rhombicuboctahedron was analytically derived in Ref. [8], while its shear modulus is not determined.

For the Poisson’s ratio, there is a good agreement between the results (computational model and analytical study [8]) for densities below 0.15, but they tend to deviate at higher densities. This is illustrated in Fig. 6.

The above expression was derived for low-density $(\rho \xaf\u22640.2)$ metallic lattice materials. The computational model was normalized with respect to the elastic modulus of the solid material (*E*_{s}), but it still depends on its Poisson’s ratio (*ν*_{s}). Nonetheless, the *ν*_{s} = 0.3 is often the case for metals, and Eq. (63) is valid for most metallic lattices. In a more general sense, the above expression is valid for any lattice made of an isotropic material with a Poisson’s ratio (*ν*_{s}) close to 0.3. This means that Eq. (63) can also be employed for some classes of polymeric lattice materials. Nonetheless, there exist an inherent error in curve fitting, and if the exact value of the shear modulus is of the interest, one should run the computational code.

### 5.8 Rhombicuboctahedron (Type b).

*z*-axis. This new configuration of the rhombicuboctahedron results in a different symmetry class and a new set of elastic properties, which has not been studied before. The results of the homogenization showed that the rhombicuboctahedron (type b) is a BDLM with tetragonal symmetry, with six independent elastic constants. By fitting curves over data points obtained from the computational model, the effective engineering constants are given by

The above expressions are valid for lattices with arbitrary elastic modulus for solid material. However, these results are specific for solid materials with Poisson’s ratio *ν*_{s} = 0.3. Accordingly, these results are applicable for low-density lattices $(\rho \xaf\u22640.2)$ made of isotropic material with *ν*_{s} = 0.3 (typically metals and some classes of polymers).

Despite the lattice mixed mode deformation in all strain states, the results of the study showed that at low relative densities, the effect of bending (reflected as term with $\rho \xaf2$) is negligible for the third and sixth strain states. In other words, the lattice out-of-plane elastic modulus (*E*_{3}) and in-plane shear modulus (*G*_{12}) can be approximated with linear functions (as for SDLM). However, this is not the case for the shear moduli *G*_{23} and *G*_{13}, and elastic moduli *E*_{1} and *E*_{2}, in which the effect of bending mode is influential.

### 5.9 Double-Pyramid Dodecahedron.

Double-pyramid dodecahedron is a relatively new lattice. The in-plane elastic moduli of this lattice was characterized in Ref. [10]. However, the provided expression was in terms of nodal forces and displacements, which requires the reader to solve the problem analytically or numerically before using the expression. Moreover, the lattice Poisson’s ratio, shear modulus, and out-of-plane elastic modules have not been reported.

*l*, the out-of-plane inclined members have length $2l$, and the members extended from the vertices of the hexagon have the length

*l*/2. Assuming all members have the same cross-sectional area, the relative density of the lattice reads

*ν*

_{s}= 0.3.

To illustrate the deformation modes of this lattice, the deformed shapes of the lattice are shown in Fig. 7. Although the lattice has mixed mode deformation in all strain states, the results of the study showed that at low relative densities, the effect of bending ($\rho \xaf2$) is negligible for the first three strain states and stretching is the dominant mode of deformation, see Figs. 7(a)–7(c). In other words, the three elastic moduli can be well approximated with linear functions (as for SDLM). However, this is not the case for the shear moduli *G*_{23} and *G*_{13}, in which the effect of bending mode is influential, see Figs. 7(d) and 7(e). Also, it is observed that *ν*_{21} depends on the relative density, while *ν*_{31} is almost constant, i.e., *ν*_{31} ≃ 0.26 at low densities. It is also noted that there is no active twisting mode in shear deformations of this lattice.

## 6 Directional Behavior of Anisotropic Lattices

Having identified the symmetry class and the effective stiffness matrices of the selected lattices, in this section, their elastic modulus in different directions are illustrated in Fig. 8. It is observed that the behavior of lattice materials can vary considerably, even if they are in the same symmetry class. For example, the Kelvin cell has a cubic symmetry, but it has a sphere-like directional graph with very small variations, which implies that the unit-cell have semi-isotropic behavior at low densities and can be regarded as an isotropic material at very low densities. On the other hand, the elastic modulus of both cube and rhombicuboctahedron (type a) unit-cells differ considerably in different directions. In that sense, they can be regarded as the extreme cases of the cubic symmetry class. Finally, octet and diamond lattices both have a typical directional behavior of a material with cubic symmetry with moderate variation of elastic modulus in 3D space.

## 7 Elasticity of Selected SDLM Lattices

In this section, the simplified model (described in Sec. 3) is tested on a range of 2- and 3D SDLM.

### 7.1 Three-Dimensional SDLM.

The simplified computational model, i.e., Eqs. (28) and (29), is tested here for a selected three-dimensional SDLM, i.e., regular octet, Gurtner–Durand, and irregular tetrahedron. The effective elastic properties obtained with this model are consistent with those obtained in Secs 5.1, 5.2, and 5.3. This highlights the absence of bending modes in the SDLM. Yet, it is noted that some SDLM lattices require special treatment when using this simplified model. For instance, to create a space-filling unit-cell for the G–D lattice, some of the struts should be cut by the unit-cell (RVE) faces. These half struts form virtual nodes on the unit-cell boundary that do not exist in the periodic structure. In contrast to regular nodes which acts as pins, the virtual nodes are rigid nodes and should not be modeled as pin-joints in the finite element formulation. Thus, extra constraints are to be defined on such nodes to make them rigid. However, this is not an issue in the generalized model because all nodes are automatically modeled as rigid joints.

Moreover, it is worth mentioning that the extensions to the simplified computational model, i.e., stress approach and surface averaging models, are tested for the octet lattice material. They both yielded consistent results, which are in full agreement with those obtained using the general model (Sec. 5.1) as well as the simplified model.

### 7.2 Two-Dimensional SDLM.

The dimensions of all other matrices should be adopted accordingly to tailor the homogenization scheme for a 2D problem. Here, the computational scheme is implemented on five 2D topologies including equilateral triangle, Kagome, mix-square type a, mix-square type b, and diamond. Similar to 3D lattices, a virtual node is to be considered where a strut is cut by the unit-cell boundary. For the first four cases, the elastic properties obtained by computational study agree with those provided in Ref. [6], but this is not the case for 2D diamond. Consider a 2D diamond lattice (see Fig. 1) for which all the struts have identical length (*l*) and width (*t*). The relative density of the lattice is $\rho \xaf=5t/(3l)$, and its elastic properties are characterized as $E1=0.2(Es\rho \xaf)$, $E2=0.35(Es\rho \xaf)$, *ν*_{12} = 0.33, *ν*_{21} = 0.6, and $G=0.15(Es\rho \xaf)$. More specifically, the expressions obtained for *E*_{1}, *ν*_{12}, and *G* agree with those reported in Ref. [6] but those for *E*_{2} and *ν*_{21} are different. Considering the symmetry of the compliance matrix, the values given in this study are reasonable.

Aside from beam-based analytical studies, comparing the results of the current model with Ref. [17] (which is based on the asymptotic homogenization using 2D elements) reveals that there is a very good agreement between the results for densities below 0.2 for all geometries.

The directional behavior of mix-square type a, mix-square type b, and diamond lattices are plotted in Fig. 9 on polar graphs. The results of these directional graphs are also interpreted in Table 3, where the angle *θ* is a counterclockwise angle representing the rotation of the coordinate system from its initial configuration, i.e., the configuration aligned with the material’s axes of symmetry.

Lattice | Effective properties | Strongest direction $(\theta deg)$ | Weakest direction $(\theta deg)$ |
---|---|---|---|

Mix-square A | E_{1} = E_{2} | $45deg,135deg,225deg$ and 315 deg | $0deg,90deg,180deg$ and 270 deg |

G | $0deg,90deg,180deg$ and 270 deg | $45deg,135deg,225deg$ and 315 deg | |

Mix-square B | E_{1} = E_{2} | $0deg,90deg,180deg$ and 270 deg | $45deg,135deg,225deg$ and 315 deg |

G | $45deg,135deg,225deg$ and 315 deg | $0deg,90deg,180deg$ and 270 deg | |

Diamond | E_{1} | $60deg,120deg,240deg$ and 300 deg | $0deg$ and 180 deg |

E_{2} | $30deg,150deg,210deg$ and 330 deg | $90deg$ and 270 deg | |

G | $0deg,90deg,180deg$ and 270 deg | $45deg,135deg,225deg$ and 315 deg |

Lattice | Effective properties | Strongest direction $(\theta deg)$ | Weakest direction $(\theta deg)$ |
---|---|---|---|

Mix-square A | E_{1} = E_{2} | $45deg,135deg,225deg$ and 315 deg | $0deg,90deg,180deg$ and 270 deg |

G | $0deg,90deg,180deg$ and 270 deg | $45deg,135deg,225deg$ and 315 deg | |

Mix-square B | E_{1} = E_{2} | $0deg,90deg,180deg$ and 270 deg | $45deg,135deg,225deg$ and 315 deg |

G | $45deg,135deg,225deg$ and 315 deg | $0deg,90deg,180deg$ and 270 deg | |

Diamond | E_{1} | $60deg,120deg,240deg$ and 300 deg | $0deg$ and 180 deg |

E_{2} | $30deg,150deg,210deg$ and 330 deg | $90deg$ and 270 deg | |

G | $0deg,90deg,180deg$ and 270 deg | $45deg,135deg,225deg$ and 315 deg |

## 8 Conclusion

Several computational models were developed in this study for the elastostatic analysis of generally anisotropic lattice materials. Lattices from any arbitrary symmetry class, including both SDLM and BDLM, were analyzed and characterized.

In the first part of the paper, a general computational homogenization model was developed for low-density lattice materials. The homogenization scheme is based on rigid nodes and general beam elements with PBC and employs the strain-based averaging technique. The Hill–Mandel Lemma was employed to derive the properties and a shear correction factor was introduced for constitutive approach. This was followed by suggesting a simplified model specifically for the analysis of SDLM, which can further reduce the computational time. Then, a new criterion based on computational homogenization (using periodicity and motion constraints) was introduced for the classification of lattice materials. Finally, two other homogenization techniques, i.e., stress-based averaging with SUBC and strain-based surface averaging with PBC, were introduced as extensions to the simplified computational model.

In the second part, the developed computational models were implemented on a wide range of lattice materials. More specifically, the general model was tested for nine distinct 3D lattices with different symmetries. The obtained results were validated with the previous studies, where applicable. Four of the studied topologies are fully characterized for the first time in the literature. Once needed, monoclinic distance function was considered to identify the anisotropy of the lattice. Next, the simplified model was tested on three 3D topologies and five 2D topologies. While the obtained results for most of the cases were consistent with the available literature, a modification was suggested for the effective properties of the 2D diamond lattice. Additionally, anisotropic behavior of selected 2D topologies was analyzed, and the stiffest/weakest material directions were identified. Eventually, the extensions to the simplified model were tested for the octet lattice topology and provided identical results, which agreed with the literature as well as with the strain-based volume averaging method.

The computational model developed in this study is fully automated and has no limitation on geometrical complexity or symmetry class. In contrast to previous computational models which were based on 3D elements, the current model is developed using 1D elements and is computationally far more efficient when it comes to low-density lattice materials. This study aimed for linear elastic analysis within classical continuum mechanics. The higher-order and higher-grade homogenization of lattice materials toward generalized continua and the corresponding anisotropy can be a topic for future studies.

## Acknowledgment

The authors acknowledge the financial support by the Starting Grant (2018-03636) from the Swedish Research Council (Vetenskapsrådet).

### Appendix A: Strain Transformation from Cylindrical to Cartesian System

*r*,

*θ*,

*x*) and (

*x*,

*y*,

*z*), respectively, while the

*x*axis is along the axis of the 1D element. The transformation matrix from cylindrical to the Cartesian system is defined by

*x*axis), the strain vector in cylindrical coordinate is given as

*γ*

_{θx}. Rewriting [

*ɛ*]

_{polar}in the matrix form yields

*ɛ*]

_{polar}can be transformed into the Cartesian coordinate system using the transformation

*R*matrix as according to

*ɛ*]

_{cart}, is now rewritten in the vector from to obtain the strain vector of the shaft element in the Cartesian coordinate as

*γ*

_{xθ}=

*r*((

*φ*

_{x2}−

*φ*

_{x1})/

*L*).

### Appendix B: Stress Transformation From Local to Global Coordinate System

*x*′,

*y*′,

*z*′) and (

*x*,

*y*,

*z*), respectively. The field variables for elements are often obtained in their local coordinate system, in which the

*x*′ is aligned with the strut axial direction. Thus, to find the properties in global coordinate system, one should use a transformation matrix defined by

*T*

_{s}] can be defined as

*R*

_{ij}corresponds to the components of the [R] matrix. (Please note that the [R] defined here is different from the [R] matrix used in periodicity equations.) The stress and strain vectors are transformed according to