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. [2428] 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.

Consider the RVE of a two-phase heterogeneous material (a porous media in this case) defined in region 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]:
(1)
where ei (i = 1, 2, 3) are the unit basis vectors.
In a linear elastic problem, the microscopic strain field ɛ can be related to an arbitrary macroscopic strain field ε¯, using a localization matrix [M]. This can be written in the matrix form (Voigt notation) as
(2)
where
(3)

Here [B] is the derivative of the shape function (also referred to as gradient matrix) and [de] is the element nodal deformation matrix in local coordinate system, corresponding to the six possible strain states, see Eq. (1). The matrix [de] 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.

Assume ui, vi, and wi 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
(4)
For each element, the [d] matrix should be calculated for all six macroscopic strain states and stored as the column matrices to form [de]. Hence, [de] is a 12 × 6 matrix, in which each column corresponds to one macroscopic strain state. That is
(5)
where [d](ij) is the matrix of nodal deformation corresponding to macroscopic state ε¯(ij), see relation (1). Since the nodal deformations obtained by FEM are expressed in the global coordinate system, i.e., [De], the nodal deformation in the local coordinate system reads
(6)
where [Tbeam] 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.,
(7)
Following Voigt notation and within local Cartesian coordinates of the element, these gradient matrices given by
(8)
and
(9)
and
(10)
where, r=y2+z2, cos(θ)=y/y2+z2 and sin(θ)=z/y2+z2. Here, [Bshaft], 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.
Following the definition of a general beam element, the element’s strain vector can be written by superposing the strain vectors as
(11)
Consequently, the stress field can be obtained as
(12)
where [σ] is the microscopic stress vector of the element and [Ce] is the element elasticity matrix defined in its local coordinate system. For a general beam element, all components of the [Ce] are zero except the (1,1) which is equal the elastic modulus (Es) as well as (5,5) and (6,6) components which are equal to the shear modulus (Gs) of the solid material.
The effective elasticity matrix can be derived using Hill–Mandel Lemma (energy approach) given by
(13)
where [σ¯] is the macroscopic stress vector. It is pointed out that the stress vector defined in Eq. (12) does not include the stresses resulting from shear deformation. This is due to the simplifying assumptions in E–B theory, i.e., neglecting shear deformations (γzx = γxy = 0). Yet, in this method, the effect of shear forces (reflected in beam bending) is captured through bending deformation.
Using Eqs. (12) and (13), the effective macroscopic elasticity matrix is obtained by
(14)
The superscript 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
(15)
where VRVE is the volume of the RVE and Ne 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 ≤ xL, r2z2y+r2z2 and −rz ≤ + 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.

There is an alternative approach for deriving the effective elasticity matrix, which is referred to as constitutive approach in this paper. In contrast to energy method that is based on the deformation energy, the constitutive approach is only based on the stress field. However, as was mentioned earlier, the stress field obtained by E–B theory does not contain the stress terms resulting from shear forces and cannot capture shear effects on its own, i.e., without considering the strain field (deformations). This can be treated by separately adding the shear stress matrix of the element resulting from shear forces, [τs], along with a shear correction factor, ks. The element’s modified stress vector reads as
(16)
Defining the following shear stress matrix
(17)
we have
(18)
where, fy(ij) and fz(ij) are the shear forces along element’s local y and z axes corresponding to macroscopic state ε¯(ij), and A is the area of the element cross section.
By averaging the Eq. (16) over the RVE volume and using constitutive equation, one can derive the effective stiffness matrix as
(19)

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 (ks) turned out to be ks = 2.

Considering the discrete nature of the problem, the above integral can be discretized as
(20)

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 [Ts] (see Appendix  B). Such transformation is not needed in Eq. (15), where the energy is used for homogenization.

The nodal deformation matrix [de] appeared in Eqs. (15) and (20) is to obtained using the finite element methodology. This procedure is explained in the following section.

2.2 Finite Element Analysis.

The matrices of nodal deformations and nodal forces are calculated by FEM. It is possible to solve the FE problem using different types of boundary conditions. However, in order to effectively consider the macroscopic behavior of an infinite lattice material, periodic boundary conditions (PBC) are used. For a set of paired nodes, e.g., A and B being on two opposite sides of RVE, the periodicity condition reads [14]
(21)
where 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
(22)
where [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 [de] to be inserted in Eqs. (15) and (20) to derive the effective elasticity matrix of the material.

The system of equations takes the form
(23)
where [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, Nn is the number of nodes and Np 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 matrix of external nodal forces reads
(24)

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.

For SDLM, the nodal deformation matrix, [d], reduces to the nodal displacement matrix, [u], defined as [u]=[u1v1w1u2v2w2]T. Similar to [de] 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 [ue]. Additionally, the gradient matrix described in Eq. (7) reduces to [B] = [Bbar], and the element’s strain vector reads
(25)
where the nodal displacement in local coordinates is given by [ue] = [Tbar][Ue]. Here, [Ue] is the nodal displacement matrix in the global coordinate system obtained by FEM, and [Tbar] is the transformation matrix for the bar element [33].
Since the struts are only loaded axially, there is no shear effects in SDLM and thus, the stress vector is simply defined as [σe]=[Ce][Bbar][ue][ε¯]. For bar elements, [Ce] is equal to the elastic modulus (Es) of the solid material. Finally, following the same procedure described in Sec. 2.1, the elasticity matrix based on the constitutive approach reads
(26)
According to Hill–Mandel Lemma, the effective elasticity matrix alternatively reads
(27)
Since the lattice unit-cell is composed of discrete elements with constant [Bbar] matrix, Eqs. (26) and (27) can be written in the form of a summation series as
(28)
and
(29)
where Ve 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 = sm, 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 [Aeq][tstrut] = [f], where [Aeq] is the equilibrium matrix containing the cosine vector of the members with respect to the global coordinate system, [tstrut] 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=brAeq, 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, [Aeq], 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.

Fig. 1
Schematic of 2D unit-cells with x and y axes along the edges of the unit-cell: (a) triangular, (b) Kagome, (c) mix-square type a, (d) mix-square type b, (e) mix-square type c, (f) diamond, (g) square web, (h) square, and (i) hexagonal
Fig. 1
Schematic of 2D unit-cells with x and y axes along the edges of the unit-cell: (a) triangular, (b) Kagome, (c) mix-square type a, (d) mix-square type b, (e) mix-square type c, (f) diamond, (g) square web, (h) square, and (i) hexagonal
Close modal
Fig. 2
Schematic of the selected lattice RVEs: (a) Octet, (b) tetrakaidekahedron, (c) irregular tetrahedron, (d) cube, (e) diamond, (f) Kelvin, (g) rhombicuboctahedron, and (h) double-pyramid dodecahedron. (Color version online.)
Fig. 2
Schematic of the selected lattice RVEs: (a) Octet, (b) tetrakaidekahedron, (c) irregular tetrahedron, (d) cube, (e) diamond, (f) Kelvin, (g) rhombicuboctahedron, and (h) double-pyramid dodecahedron. (Color version online.)
Close modal

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.

Table 1

Determinacy analysis of selected 2D topologies

TopologysmMPeriodic behavior
1Triangular000SD
2Kagome03−3SD
3Mix-square type a101SD
4Mix-square type b101SD
5Mix-square type c000SD
6Diamond000SD
7Square web211BD
8Square01−1BD
9Hexagonal03−3BD
TopologysmMPeriodic behavior
1Triangular000SD
2Kagome03−3SD
3Mix-square type a101SD
4Mix-square type b101SD
5Mix-square type c000SD
6Diamond000SD
7Square web211BD
8Square01−1BD
9Hexagonal03−3BD
Table 2

Determinacy analysis of selected 3D topologies

TopologysmMPeriodic behavior
1Octet000SD
2Tetrakaidekahedral114−13SD
3Tetrahedron505SD
4Cube06−6BD
TopologysmMPeriodic behavior
1Octet000SD
2Tetrakaidekahedral114−13SD
3Tetrahedron505SD
4Cube06−6BD

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.

The criterion is constructed by supplementing the stiffness matrix with periodicity and rigid-body-translation constraints. This approach is an alternative to the one based on eliminating the DOF of the dependent joints, which leads to reduced equilibrium equations [21]. Recall the Eq. (23) in Sec. 2.2. This relation can be adopted for bar elements (pin-joint assumption), by replacing the matrix of nodal deformations, [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
(30)
We refer to this relation as “periodicity equilibrium” equation which governs the behavior of a lattice material as a periodic medium and not as a finite unit-cell. The square coefficient matrix appeared on the left-hand side is defined as periodicity equilibrium matrix [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.,
(31)
where 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.

Consider the RVE described in Sec. 2.1. In a linear elastic problem, the microscopic stress field σ can be related to an arbitrary macroscopic strain field σ¯, using a localization matrix [N] (see Appendix  C for different macroscopic stress states). This can be written in the matrix form as [14]
(32)
By multiplying the above equation to the compliance matrix, the strain vector reads
(33)
By averaging over RVE domain, we have
(34)
where, 〈…〉 denotes averaging over the RVE volume. Considering the macroscopic constitutive law, the effective elastic tensor is given by
(35)
In the stress-based approach, statically uniform boundary condition (SUBC) (also referred to as Neumann boundary condition) should be satisfied, i.e.,
(36)
where n is the normal unit vector of the RVE surfaces.
Lattice materials are two-phase materials consisting of a solid phase and voids. Since the strain is not defined within voids, it is not possible to spatially average the strain over the RVE volume in Eq. (34). However, using the divergence theorem, it is possible to transform volume integration into surface integration [14] according to
(37)
where 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, ε¯. These matrices are written in vector form and stored as the columns of [ε¯total] as
(38)
where [ε¯](ij) are the macroscopic strain vectors obtained from Eq. (37) (corresponding to the applied [σ¯](ij)). By rewriting the six macroscopic stress matrices in vector from and arranging them in the same manner, the matrix [σ¯total] is obtained as
(39)
Ultimately, considering Hook’s law, the effective elasticity matrix reads
(40)

4.2 Strain Approach With Surface Averaging.

In Sec. 3, the effective properties were obtained by taking the average of properties over the RVE volume under macroscopic strain fields. However, it is possible to use the same methodology, i.e., strain approach, but derive the effective properties by averaging over the RVE surface. Using divergence theorem, it is possible to transform volume integration into surface integration [14] according to
(41)
where 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 σ¯ 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 [σ¯total] can be formed. The macroscopic strain matrices are also written in vector form to give [ε¯total]. 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=(FTFI)/2=const. This ultimately results in a linear displacement field (called affine displacement field), i.e., U=EX, in which X is the position of the material points.

The deformation field is homogenous (affine) if and only if the undeformed configuration (position) of the material points can be uniquely determined. This is the case for SDLM, as their deformation is characterized by pure stretching, and thus, the undeformed configuration can be uniquely determined leading to an affine microscopic displacement field, which coincides with the macroscopic displacement field. Consequently, the six RVE planes in SDLM remain plane after deformation, and by considering the linear displacement field in the RVE planes, the integrals in Eqs. (37) and (41) can, respectively, be simplified to
(42)
and
(43)
where A is the area of the RVE faces (identical for cubic RVE), ni is the outward unit normal vector of the face i, uci is the displacement of the centroid of the face i, Fi 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 ρ¯0.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.

The regular octet lattice is a classical 3D lattice with high stiffness-to-weight ratio and is an interesting candidate for load bearing applications. It is a well-established SDLM in the literature, e.g., see Refs. [5,13], and Ref. [16], and is a good candidate for validation purpose. Assuming all the struts have length (l) and cross-section area (A), the relative density of the lattice reads
(44)
The results of the computational study showed that the octet lattice has cubic symmetry, with three independent elastic constants. The lattice engineering constants are obtained as
(45)
where Es 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 A1 and A2 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.

The relative density of the lattice can be characterized by the length of the unit-cell edge (l) as
(46)
The results of the computational study showed that the G–D lattice is isotropic with only two engineering constants. The lattice effective properties can be normalized as
(47)

The results agree with those reported by Gurtner and Durand [12] and Norris [13].

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.

Assuming all the struts have the same cross-sectional area (A), the relative density of the lattice can be characterized by the length of the unit-cell edge (l) as
(48)
The tetrahedron lattice is computationally homogenized using the algorithm described in the previous sections (for ρ¯=0.16). The resulting stiffness matrix in the original coordinate system (defined on the unit-cell in Fig. 2) is given by
(49)
As the original coordinate system corresponding with the RVE edges is not the material natural coordinate system, the material symmetry of the lattice is not obvious in this stiffness matrix. To identify the symmetry class of this lattice, the monoclinic distance function is constructed using the methodology described in Ref. [39] as
(50)

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.

Fig. 3
Monoclinic distance function (MPa2) for tetrahedron lattice
Fig. 3
Monoclinic distance function (MPa2) for tetrahedron lattice
Close modal
Fig. 4
Contours of the monoclinic distance function for tetrahedron lattice
Fig. 4
Contours of the monoclinic distance function for tetrahedron lattice
Close modal
Studying the patterns in Figs. 3 and 4 and locating the roots of the monoclinic distance function, three mirror planes are obtained for the material. The unit normal vectors of these mirror planes, with respect to the original coordinate system, are given by
(51)
These normal vectors are three coplanar vectors equally spaced with 60deg angles and characterize the trigonal symmetry class. To find the natural coordinate system based on the material principal direction, the 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
(52)
This matrix highlights the symmetry of the lattice, i.e., trigonal symmetry with six independent constants. To find the engineering constants, the corresponding compliance matrix is characterized as
(53)
where ηi,jk and μij,kl are the effective interaction and Chentsov coefficients, respectively. From the above matrix, the effective engineering constants are obtained as
(54)
with ν13/E1 = ν31/E3, ν23/E2 = ν32/E3, G12 = 1/2((1/E1) + (ν21/E2))−1, and μ23,12/G12 = −2(η1,13/G13). To the authors’ knowledge, this lattice is characterized for the first time in the literature.

5.4 Regular Cube.

Regular cubic lattice is a BDLM with high stiffness in tensile loading and very low stiffness in shear loading. This lattice has been extensively studied in the literature, e.g., see Refs. [13] and [16], and is a good choice for the assessment of the developed computational model for BDLM. Once all the struts have length (l) and cross-section area (A), the relative density of the lattice is given by
(55)

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 ρ¯. However, the elastic properties of BDLM may scale dissimilarly under various macroscopic strains states ε¯(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 ρ¯2.

The results of the computational study showed that the regular cube lattice possess cubic symmetry, with three independent elastic constants. The regular cube deforms axially under first three strains states and its elastic modulus scales with ρ¯ (as for SDLM), while it bends in share states and its shear modulus scales with ρ¯2(as for BDLM). Therefore, the lattice engineering constants are characterized as
(56)

The above equations are consistent with the findings of Refs. [13] and [16]. Additionally, from the expression for shear modulus, it can be concluded that the low-density cubic lattice behaves similar to Penta mode materials, in a sense that its shear modulus at low densities is negligible.

5.5 Diamond.

Diamond lattice is another example of BDLM. The diamond lattice has been analyzed in several studies, e.g., Refs. [13] and [14], and is brought here for validation. Assuming all the struts have length (l) and cross-section area (A), the relative density of the lattice is
(57)
The results of the computational study showed that the diamond lattice has cubic symmetry, with three independent elastic constants. Unlike the regular cube in which the deformation modes were separable, the diamond lattice has mixed mode deformation, i.e., both stretching and bending, in all strain states. For this, the characterization of elastic properties with simple equations is not possible. Yet, using the results of the computational study along with the analytical expressions presented in Ref. [13], the effective engineering constants of diamond lattice are characterized as below:
(58)

For very low densities, i.e., ρ¯0.1, the above expressions can be simplified as E(1/3)Esρ¯2, ν ≃ 0.5 and G(1/3)Esρ¯2. 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.

Kelvin cell, also referred to as truncated octahedron or tetrakaidekahedron, is a class of Archimedean solids which can fill the space with no gaps. This cell nearly satisfies the minimum surface energy and can reasonably represent the unit-cell of an open-cell foam [40]. Its elastic properties has been studied earlier, e.g., see Refs. [7,15,40] and Ref. [41]. Here, the selected RVE, see Fig. 2, for elastostatic analysis of Kelvin lattice is different from its unit-cell. Assuming all the struts have length (l) and cross-section area (A), the relative density of the lattice is
(59)

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.

The results of the computational study showed that the Kelvin lattice has cubic symmetry, with three independent elastic constants. The equivalent engineering constants can be characterized as
(60)
where cf 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 cf = (8 EsI + GsJ)/(5 EsI + GsJ), where Gs 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, ρ¯0.1, the above elastic properties can be approximated as E0.6Esρ¯2, ν ≃ 0.5, and G0.6Esρ¯2/cf. The obtained elastic properties in this computational study agreed with those reported in Refs. [7] and [40].

5.7 Rhombicuboctahedron (Type a).

Rhombicuboctahedron is another class of Archimedean solids, with promising application in biomedical industry, especially for bone replacement [8]. The schematic of its RVE, which is identical to its unit-cell, is shown in Fig. 2. The relative density of the lattice reads
(61)
where 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.
Fig. 5
Spatial configuration of rhombicuboctahedron type a and type b. For both types, the colored unit-cells are obtained by tessellation of the blue (original) unit-cell along their periodicity directions.
Fig. 5
Spatial configuration of rhombicuboctahedron type a and type b. For both types, the colored unit-cells are obtained by tessellation of the blue (original) unit-cell along their periodicity directions.
Close modal

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.

The results of the computational study revealed that the rhombicuboctahedron (type a) is a BDLM lattice with cubic symmetry, with three independent elastic constants. The elastic modulus predicted by the computational model is very close to the results of [8], and thus can be characterized as
(62)

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.

Fig. 6
Comparison of the effective Poisson’s ratio of rhombicuboctahedron lattice obtained by computational study with the closed-from analytical expressions provided in Ref. [8]
Fig. 6
Comparison of the effective Poisson’s ratio of rhombicuboctahedron lattice obtained by computational study with the closed-from analytical expressions provided in Ref. [8]
Close modal
No analytical solution is available for the shear modulus of rhombicuboctahedron (type a). In computational analysis, the effective properties are calculated in an automated process using an in-house matlab code. The results are eventually obtained numerically and therefore as expected for a numerical method, it is not possible to present them as closed-form analytical expressions, if no prior knowledge exist on the type/form of the characterization function. For that, the shear modulus of the lattice is characterized by curve fitting according to the results of the computational study. This gives the shear modulus as
(63)

The above expression was derived for low-density (ρ¯0.2) metallic lattice materials. The computational model was normalized with respect to the elastic modulus of the solid material (Es), 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).

The rhombicuboctahedron (type b) is different from (type a) in the sense that its periodicity directions have rotated by 45deg around the 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
(64)

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 (ρ¯0.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 ρ¯2) is negligible for the third and sixth strain states. In other words, the lattice out-of-plane elastic modulus (E3) and in-plane shear modulus (G12) can be approximated with linear functions (as for SDLM). However, this is not the case for the shear moduli G23 and G13, and elastic moduli E1 and E2, 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.

The schematic of the RVE of the lattice is illustrated in Fig. 2. The edges of the hexagonal have length 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
(65)
The results of the computational study showed that the double-pyramid dodecahedron is a BDLM with transverse isotropy and thus have five independent elastic constants. By fitting curves over data points obtained using computational model, the effective engineering constants read
(66)
for low-density lattices (ρ¯0.2) made of an isotropic solid material with ν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 (ρ¯2) 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 G23 and G13, 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.

Fig. 7
Deformed shape of double-pyramid dodecahedron under various strain states. The solid lines indicate the undeformed unit-cell, while dashed lines represent the deformed lattice. (a) state ɛ11, (b) state ɛ22, (c) state ɛ33, (d) state ɛ23, (e) state ɛ13, and (f) state ɛ12.
Fig. 7
Deformed shape of double-pyramid dodecahedron under various strain states. The solid lines indicate the undeformed unit-cell, while dashed lines represent the deformed lattice. (a) state ɛ11, (b) state ɛ22, (c) state ɛ33, (d) state ɛ23, (e) state ɛ13, and (f) state ɛ12.
Close modal

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.

Fig. 8
Directional behavior of the effective elastic modulus (E11) for the selected lattices with ρ¯≃11.5%, Es = 210 GPa and νs = 0.3 (measured in GPa): (a) Octet (cubic symmetry), (b) non-regular tetrahedron (trigonal symmetry), (c) cube (cubic symmetry), (d) diamond (cubic symmetry), (e) Kelvin (cubic symmetry), (f) rhombicuboctahedron A (cubic symmetry), (g) rhombicuboctahedron B (tetragonal symmetry), and (h) Double-pyramid dodecahedron (transverse isotropy)
Fig. 8
Directional behavior of the effective elastic modulus (E11) for the selected lattices with ρ¯≃11.5%, Es = 210 GPa and νs = 0.3 (measured in GPa): (a) Octet (cubic symmetry), (b) non-regular tetrahedron (trigonal symmetry), (c) cube (cubic symmetry), (d) diamond (cubic symmetry), (e) Kelvin (cubic symmetry), (f) rhombicuboctahedron A (cubic symmetry), (g) rhombicuboctahedron B (tetragonal symmetry), and (h) Double-pyramid dodecahedron (transverse isotropy)
Close modal

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.

In addition to the selected 3D topologies in Sec. 7.1, the simplified computational model is tested here for a range of 2D SDLM. The methodology is identical to the 3D case, except that the macroscopic strain states reduce to the following three in-plane states
(67)

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 ρ¯=5t/(3l), and its elastic properties are characterized as E1=0.2(Esρ¯), E2=0.35(Esρ¯), ν12 = 0.33, ν21 = 0.6, and G=0.15(Esρ¯). More specifically, the expressions obtained for E1, ν12, and G agree with those reported in Ref. [6] but those for E2 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.

Fig. 9
Directional elastic behavior of (a) mix-square type a, (b) mix-square type b, and (c) diamond lattice; for ρ¯=0.2 and Es = 210 GPa (the effective properties are measured in GPa). For colors, see the electronic version of the article: (a) mix-square type a, (b) mix-square type b, and (c) diamond.
Fig. 9
Directional elastic behavior of (a) mix-square type a, (b) mix-square type b, and (c) diamond lattice; for ρ¯=0.2 and Es = 210 GPa (the effective properties are measured in GPa). For colors, see the electronic version of the article: (a) mix-square type a, (b) mix-square type b, and (c) diamond.
Close modal
Table 3

Stiffest/weakest direction for selected lattice topologies

LatticeEffective propertiesStrongest direction (θdeg)Weakest direction (θdeg)
Mix-square AE1 = E245deg,135deg,225deg and 315 deg0deg,90deg,180deg and 270 deg
G0deg,90deg,180deg and 270 deg45deg,135deg,225deg and 315 deg
Mix-square BE1 = E20deg,90deg,180deg and 270 deg45deg,135deg,225deg and 315 deg
G45deg,135deg,225deg and 315 deg0deg,90deg,180deg and 270 deg
DiamondE160deg,120deg,240deg and 300 deg0deg and 180 deg
E230deg,150deg,210deg and 330 deg90deg and 270 deg
G0deg,90deg,180deg and 270 deg45deg,135deg,225deg and 315 deg
LatticeEffective propertiesStrongest direction (θdeg)Weakest direction (θdeg)
Mix-square AE1 = E245deg,135deg,225deg and 315 deg0deg,90deg,180deg and 270 deg
G0deg,90deg,180deg and 270 deg45deg,135deg,225deg and 315 deg
Mix-square BE1 = E20deg,90deg,180deg and 270 deg45deg,135deg,225deg and 315 deg
G45deg,135deg,225deg and 315 deg0deg,90deg,180deg and 270 deg
DiamondE160deg,120deg,240deg and 300 deg0deg and 180 deg
E230deg,150deg,210deg and 330 deg90deg and 270 deg
G0deg,90deg,180deg and 270 deg45deg,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

Consider a cylindrical and a Cartesian coordinate system characterized with (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
(A1)
If a torsion moment is applied around the element axis (i.e., x axis), the strain vector in cylindrical coordinate is given as
(A2)
in which, the only non-zero component is γθx. Rewriting [ɛ]polar in the matrix form yields
(A3)
The strain matrix [ɛ]polar can be transformed into the Cartesian coordinate system using the transformation R matrix as according to
(A4)
The strain matrix in Cartesian coordinate, [ɛ]cart, is now rewritten in the vector from to obtain the strain vector of the shaft element in the Cartesian coordinate as
(A5)
where γ = r((φx2φx1)/L).

Appendix B: Stress Transformation From Local to Global Coordinate System

Consider two sets of Cartesian coordinate system which have a mutual origin, corresponding to local and global coordinate systems characterized with (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
(B1)
This matrix can be used when the stress or strains are given in tensor form. However, it is a common practice to write the constitutive equation in vector from for computational purpose. In that case, a transformation matrix [Ts] can be defined as
(B2)
where Rij 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
(B3)

Appendix C: Macroscopic Stress States

In stress approach, the macroscopic stress fields, σ¯(ij), are prescribed as
(C1)

References

1.
Queheillalt
,
D. T.
, and
Wadley
,
H. N. G.
,
2005
, “
Cellular Metal Lattices With Hollow Trusses
,”
Acta Mater.
,
53
(
2
), pp.
303
313
. 10.1016/j.actamat.2004.09.024
2.
Molavitabrizi
,
D.
, and
Laliberte
,
J.
,
2020
, “
Methodology for Multiscale Design and Optimization of Lattice Core Sandwich Structures for Lightweight Hopper Railcars
,”
Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci.
,
234
(
21
), pp.
4224
4238
. 10.1177/0954406220920694
3.
Ashby
,
M. F.
,
2006
, “
The Properties of Foams and Lattices
,”
Philos. Trans. R. Soc. Math. Phys. Eng. Sci.
,
364
(
1838
), pp.
15
30
. 10.1098/rsta.2005.1678
4.
Yan
,
J.
,
Cheng
,
G.
,
Liu
,
S.
, and
Liu
,
L.
,
2006
, “
Comparison of Prediction on Effective Elastic Property and Shape Optimization of Truss Material With Periodic Microstructure
,”
Int. J. Mech. Sci.
,
48
(
4
), pp.
400
413
. 10.1016/j.ijmecsci.2005.11.003
5.
Deshpande
,
V. S.
,
Fleck
,
N. A.
, and
Ashby
,
M. F.
,
2001
, “
Effective Properties of the Octet-Truss Lattice Material
,”
J. Mech. Phys. Solids
,
49
(
8
),, p.
23
.
6.
Wang
,
A.-J.
, and
McDowell
,
D. L.
,
2004
, “
In-Plane Stiffness and Yield Strength of Periodic Metal Honeycombs
,”
ASME J. Eng. Mater. Technol.
,
126
(
2
), pp.
137
156
. 10.1115/1.1646165
7.
Warren
,
W. E.
, and
Kraynik
,
A. M.
,
1997
, “
Linear Elastic Behavior of a Low-Density Kelvin Foam With Open Cells
,”
ASME J. Appl. Mech.
,
64
(
4
), pp.
787
794
. 10.1115/1.2788983
8.
Hedayati
,
R.
,
Sadighi
,
M.
,
Mohammadi-Aghdam
,
M.
, and
Zadpoor
,
A. A.
,
2016
, “
Mechanics of Additively Manufactured Porous Biomaterials Based on the Rhombicuboctahedron Unit Cell
,”
J. Mech. Behav. Biomed. Mater.
,
53
, pp.
272
294
. 10.1016/j.jmbbm.2015.07.013
9.
Babaee
,
S.
,
Jahromi
,
B. H.
,
Ajdari
,
A.
,
Nayeb-Hashemi
,
H.
, and
Vaziri
,
A.
,
2012
, “
Mechanical Properties of Open-Cell Rhombic Dodecahedron Cellular Structures
,”
Acta Mater.
,
60
(
6
), pp.
2873
2885
. 10.1016/j.actamat.2012.01.052
10.
Mahbod
,
M.
, and
Asgari
,
M.
,
2019
, “
Elastic and Plastic Characterization of a New Developed Additively Manufactured Functionally Graded Porous Lattice Structure: Analytical and Numerical Models
,”
Int. J. Mech. Sci.
,
155
, pp.
248
266
. 10.1016/j.ijmecsci.2019.02.041
11.
Abdelhamid
,
M.
, and
Czekanski
,
A.
,
2018
, “
Impact of the Lattice Angle on the Effective Properties of the Octet-Truss Lattice Structure
,”
ASME J. Eng. Mater. Technol.
,
140
(
4
), p.
041010
. 10.1115/1.4040409
12.
Gurtner
,
G.
, and
Durand
,
M.
,
2014
, “
Stiffest Elastic Networks
,”
Proc. R. Soc. Math. Phys. Eng. Sci.
,
470
(
2164
), p.
20130611
. 10.1098/rspa.2013.0611
13.
Norris
,
A. N.
,
Dec. 2014
, “
Mechanics of Elastic Networks
,”
Proc. R. Soc. Math. Phys. Eng. Sci.
,
470
(
2172
), p.
20140522
. 10.1098/rspa.2014.0522
14.
Yvonnet
,
J.
,
2019
,
Computational Homogenization of Heterogeneous Materials With Finite Elements
,
Springer International Publishing
.
15.
Zhu
,
W.
,
Blal
,
N.
,
Cunsolo
,
S.
, and
Baillis
,
D.
,
2017
, “
Micromechanical Modeling of Effective Elastic Properties of Open-Cell Foam
,”
Int. J. Solids Struct.
,
115–116
, pp.
61
72
. 10.1016/j.ijsolstr.2017.02.031
16.
Vigliotti
,
A.
, and
Pasini
,
D.
,
2012
, “
Stiffness and Strength of Tridimensional Periodic Lattices
,”
Comput. Methods Appl. Mech. Eng.
,
229–232
, pp.
27
43
. 10.1016/j.cma.2012.03.018
17.
Arabnejad
,
S.
, and
Pasini
,
D.
,
2013
, “
Mechanical Properties of Lattice Materials via Asymptotic Homogenization and Comparison With Alternative Homogenization Methods
,”
Int. J. Mech. Sci.
,
77
, pp.
249
262
. 10.1016/j.ijmecsci.2013.10.003
18.
Dong
,
G.
,
Tang
,
Y.
, and
Zhao
,
Y. F.
,
2019
, “
A 149 Line Homogenization Code for Three-Dimensional Cellular Materials Written in Matlab
,”
ASME J. Eng. Mater. Technol.
,
141
(
1
), p.
011005
. 10.1115/1.4040555
19.
Nikolaos
,
K.
,
Francisco
,
D. R.
,
Hadjidoukas
,
P.
, and
Ganghoffer
,
J.-F.
,
2020
, “
LatticeMech: A Discrete Mechanics Code to Compute the Effective Static Properties of 2D Metamaterial Structures
,”
SoftwareX
,
11
, p.
100446
. 10.1016/j.softx.2020.100446
20.
Tancogne-Dejean
,
T.
,
Karathanasopoulos
,
N.
, and
Mohr
,
D.
,
2019
, “
Stiffness and Strength of Hexachiral Honeycomb-Like Metamaterials
,”
ASME J. Appl. Mech.
,
86
(
11
), p.
111010
. 10.1115/1.4044494
21.
Hutchinson
,
R. G.
, and
Fleck
,
N. A.
,
2006
, “
The Structural Performance of the Periodic Truss
,”
J. Mech. Phys. Solids
,
54
(
4
), pp.
756
782
. 10.1016/j.jmps.2005.10.008
22.
Elsayed
,
M. S. A.
, and
Pasini
,
D.
,
2010
, “
Analysis of the Elastostatic Specific Stiffness of 2D Stretching-Dominated Lattice Materials
,”
Mech. Mater.
,
42
(
7
), pp.
709
725
. 10.1016/j.mechmat.2010.05.003
23.
Patil
,
G. U.
, and
Matlack
,
K. H.
,
2019
, “
Effective Property Evaluation and Analysis of Three-Dimensional Periodic Lattices and Composites Through Bloch-Wave Homogenization
,”
J. Acoust. Soc. Am.
,
145
(
3
), pp.
1259
1269
. 10.1121/1.5091690
24.
Khakalo
,
S.
, and
Niiranen
,
J.
,
2019
, “
Lattice Structures as Thermoelastic Strain Gradient Metamaterials: Evidence From Full-Field Simulations and Applications to Functionally Step-Wise-Graded Beams
,”
Compos. Part B Eng.
,
177
, p.
107224
. 10.1016/j.compositesb.2019.107224
25.
Rokoš
,
O.
,
Ameen
,
M. M.
,
Peerlings
,
R. H. J.
, and
Geers
,
M. G. D.
,
2019
, “
Micromorphic Computational Homogenization for Mechanical Metamaterials With Patterning Fluctuation Fields
,”
J. Mech. Phys. Solids
,
123
, pp.
119
137
. 10.1016/j.jmps.2018.08.019
26.
Auffray
,
N.
,
Bouchet
,
R.
, and
Bréchet
,
Y.
,
2010
, “
Strain Gradient Elastic Homogenization of Bidimensional Cellular Media
,”
Int. J. Solids Struct.
,
47
(
13
), pp.
1698
1710
. 10.1016/j.ijsolstr.2010.03.011
27.
Nikolaos
,
K.
,
Francisco
,
D. R.
,
Reda
,
H.
, and
Ganghoffer
,
J.-F.
,
2018
, “
Computing the Effective Bulk and Normal to Shear Properties of Common Two-Dimensional Architectured Materials
,”
Comput. Mater. Sci.
,
154
, pp.
284
294
. 10.1016/j.commatsci.2018.07.044
28.
Nikolaos
,
K.
,
Francisco
,
D. R.
,
Diamantopoulou
,
M.
, and
Ganghoffer
,
J.-F.
,
2020
, “
Mechanics of Beams Made From Chiral Metamaterials: Tuning Deflections Through Normal-Shear Strain Couplings
,”
Mater. Des.
,
189
, p.
108520
. 10.1016/j.matdes.2020.108520
29.
Dos Reis
,
F.
, and
Ganghoffer
,
J. F.
,
2012
, “
Construction of Micropolar Continua From the Asymptotic Homogenization of Beam Lattices
,”
Comput. Struct.
,
112–113
, pp.
354
363
. 10.1016/j.compstruc.2012.08.006
30.
Munford
,
M.
,
Hossain
,
U.
,
Ghouse
,
S.
, and
Jeffers
,
J. R. T.
,
2020
, “
Prediction of Anisotropic Mechanical Properties for Lattice Structures
,”
Addit. Manuf.
,
32
, p.
101041
. 10.1016/j.addma.2020.101041
31.
Lohmuller
,
P.
,
Favre
,
J.
,
Kenzari
,
S.
,
Piotrowski
,
B.
,
Peltier
,
L.
, and
Laheurte
,
P.
,
2019
, “
Architectural Effect on 3D Elastic Properties and Anisotropy of Cubic Lattice Structures
,”
Mater. Des.
,
182
, p.
108059
. 10.1016/j.matdes.2019.108059
32.
Xu
,
S.
,
Shen
,
J.
,
Zhou
,
S.
,
Huang
,
X.
, and
Xie
,
Y. M.
,
2016
, “
Design of Lattice Structures With Controlled Anisotropy
,”
Mater. Des.
,
93
, pp.
443
447
. 10.1016/j.matdes.2016.01.007
33.
Logan
,
D. L.
,
2016
,
A First Course in Finite Element Method: Sixth Edition
, Si Edition,
Cengage Learning
,
Mason, OH
.
34.
Pellegrino
,
S.
, and
Calladine
,
C. R.
,
1986
, “
Matrix Analysis of Statically and Kinematically Indeterminate Frameworks
,”
Int. J. Solids Struct.
,
22
(
4
), pp.
409
428
. 10.1016/0020-7683(86)90014-4
35.
Guest
,
S. D.
, and
Hutchinson
,
J. W.
,
2003
, “
On the Determinacy of Repetitive Structures
,”
J. Mech. Phys. Solids
,
51
(
3
), pp.
383
391
. 10.1016/S0022-5096(02)00107-2
36.
Landau
,
L. D.
,
1989
,
Theory of Elasticity
, Ed 3 Vol
7
,
Pergamon Press
,
Oxford, UK
.
37.
Goldberg
,
M.
,
1974
, “
Three Infinite Families of Tetrahedral Space-Fillers
,”
J. Comb. Theory Ser. A
,
16
(
3
), pp.
348
354
. 10.1016/0097-3165(74)90058-2
38.
Arabnejad
,
S.
,
Burnett Johnston
,
R.
,
Pura
,
J. A.
,
Singh
,
B.
,
Tanzer
,
M.
, and
Pasini
,
D.
,
2016
, “
High-strength Porous Biomaterials for Bone Replacement: A Strategy to Assess the Interplay Between Cell Morphology, Mechanical Properties, Bone Ingrowth and Manufacturing Constraints
,”
Acta Biomater.
,
30
, pp.
345
356
. 10.1016/j.actbio.2015.10.048
39.
Diner
,
Ç
,
Kochetov
,
M.
, and
Slawinski
,
M. A.
,
2011
, “
Identifying Symmetry Classes of Elasticity Tensors Using Monoclinic Distance Function
,”
J. Elasticity
,
102
(
2
), pp.
175
190
. 10.1007/s10659-010-9272-7
40.
Zhu
,
H. X.
,
Knott
,
J. F.
, and
Mills
,
N. J.
,
1997
, “
Analysis of the Elastic Properties of Open-Cell Foams With Tetrakaidecahedral Cells
,”
J. Mech. Phys. Solids
,
45
(
3
), pp.
319
343
. 10.1016/S0022-5096(96)00090-7
41.
Jang
,
W.-Y.
,
Kyriakides
,
S.
, and
Kraynik
,
A. M.
,
2010
, “
On the Compressive Strength of Open-Cell Metal Foams With Kelvin and Random Cell Structures
,”
Int. J. Solids Struct.
,
47
(
21
), pp.
2872
2883
. 10.1016/j.ijsolstr.2010.06.014
42.
Reis
,
F. D.
, and
Ganghoffer
,
J. F.
,
2010
, “
Discrete Homogenization of Architectured Materials
,”
Tech. Mech. Eur. J. Eng. Mech.
,
30
(
1–3
), pp.
85
109
.