Anisotropy and Asymmetry of the Elastic Tensor of Lattice Materials

Abstract

Lattice materials formed by hinged springs or linear elastic bonds may exhibit diverse anisotropy and asymmetry features of the overall elastic behavior depending on their unit cell configuration. The recently developed singum model transfers the force-displacement relationship of the springs in the lattice to the stress-strain relationship in the continuum particle, and provides the analytical form of tangential elasticity. When a pre-stress exists in the lattice, the stiffness tensor significantly changes due to the effect of the configurational stress; existing methods like the lattice spring method, relying on a scalar energy equivalence, are insufficient in such situations. Instead, a tensorial homogenization method with the new definition of singum stress and strain, should be preferred. Different lattice structures lead to different symmetries of the stiffness tensors, which are demonstrated by five lattices. When all bonds exhibit the same length, regular hexagonal, honeycomb, and auxetic lattices demonstrate that the stiffness changes from an isotropic to anisotropic, from symmetric to asymmetric tensor. When the central symmetry of the unit cell is not satisfied, the primitive cell will contain more than one singums and the Cauchy–Born rule fails by the loss of equilibrium of the single singum. A secondary stress is induced to balance the singums. Displacement gradient $d_{ij}=u_{j,i}$ is proposed to replace strain in the constitutive law for the general case because $d_{12}$ and $d_{21}$ can produce different stress states. Although the hexagonal and honeycomb lattices may exhibit isotropic behavior, for general auxetic lattices, an anisotropic and asymmetric elastic tensor is obtained with the loss of both minor and major symmetry, which is also demonstrated in a square lattice with unbalanced central symmetry and a chiral lattice. The modeling procedure and results can be generalized to three dimensions and other lattices with the anisotropic and asymmetric stiffness.

Appendices

Appendix A: Determination of the Force Transfer Matrix $\textbf{R}$

Consider a truss system of N bars, which connect to one singum node at one end with another end fixed by hinges as the cutting point. All bars are assumed to exhibit the same length and elastic constants for simplicity, and the number of the bars is large enough to make the truss system stable or indeterminate. When all the cutting points are fixed, given a force $P_{i}$ on the singum node, the force transferred to the $I$th bar can be written as

$$ T_{j}^{I} = P_{i} R^{I}_{ij}, $$

(66)

which is along the bar with the direction from the singum node to the cutting point $\textbf{n}^{I}$. Here, the conventional summation of double index notation is applied to the subscript with lower case symbols only.

Consider the force in each bar is given at $\textbf{f}^{I} = \frac{V_{,\lambda}n_{i}}{2l_{p}^{0}}$. Without the loss of any generality, set the origin $\textbf{0}$ at the singum node and the cutting points are fixed at $\textbf{x}^{I}$ ($I=1,2,\ldots , N$) with $l_{p}=|\textbf{x}^{I}| =\lambda l_{p}^{0}$, so $\textbf{n}^{I}=\textbf{x}^{I}/l_{p}$.

When all cutting points are fixed, a small variation of the singum node, $d\textbf{x}$, will change the bond vectors into $\textbf{x}^{I}-d\textbf{x}$, which leads to a length change

$$ dl^{I} = -n_{i}^{I} dx_{i} $$

(67)

and an orientational change

$$ dn_{i}^{I} = \frac{n^{I}_{i}n^{I}_{j}-\delta _{ij}}{\lambda l_{p}^{0}} dx_{j} $$

(68)

for the $I$th bonds. Note that $d\lambda ^{I} = dl^{I}/l_{p}^{0}$.

The force variation for each bar can be obtained as

$$ df_{i}^{I} = \frac{V_{,\lambda \lambda}n_{i}^{I} d\lambda ^{I} + V_{,\lambda}dn_{i}^{I}}{2l_{p}^{0}} = \left [ \frac{-V_{,\lambda \lambda}n_{i}^{I} n_{j}^{I}}{2{l_{p}^{0}}^{2} } + \frac{V_{,\lambda}(n_{i}^{I} n_{j}^{I}-\delta _{ij})}{2\lambda{l_{p}^{0}}^{2} } \right ]dx_{j}, $$

(69)

which includes two parts: the first is caused by the length change and the second orientational change. The resultant force variation on the singum is

$$ df_{i} = \sum _{I=1}^{N} df_{i}^{I}. $$

(70)

Then the force transfer matrix $R_{i}^{I}$ can be determined by the classical displacement method with the following procedure in 3D case, which can be reduced to 2D case straightforwardly:

• 1. Given a unit displacement $\textbf{d}_{1}=(1,0,0)$ on the singum node, the resultant force on the singum is calculated by Eq. (70), which is a sum of N vectors in 3D, namely $\textbf{P}_{1} = \sum _{I=1}^{N} \textbf{f}^{I1}$.

• 2. Similarly, given a unit displacement $\textbf{d}_{2}=(0,1,0)$ or $\textbf{d}_{3}=(0,0,1)$ on the singum node, the resultant force on the singum can also be calculated as $\textbf{P}_{2} = \sum _{I=1}^{N} \textbf{f}^{I2}$ and $\textbf{P}_{3}= \sum _{I=1}^{N} \textbf{f}^{I3}$.

• 3. Given any displacement $\textbf{d}=(a,b,c)$, the resultant force will be written as $\textbf{P}= a\textbf{P}_{1} +b\textbf{P}_{2} + c\textbf{P}_{3}$.

• 4. Solve $a\textbf{P}_{1} +b\textbf{P}_{2} + c\textbf{P}_{3} = (-1,0,0)$ with three equations for $a, b, c$. Using $d \textbf{x} = (a,b,c)$, one can calculate force variation of each bar by Eq. (69), which defines $R^{I}_{1j} = df_{j}^{I}$.

• 5. Similarly, solve $a\textbf{P}_{1} +b\textbf{P}_{2} + c\textbf{P}_{3} = (0, -1,0)$ or $(0,0,-1)$ with three equations for $a, b, c$. One can obtain $R^{I}_{2j}$ or $R^{I}_{3j}$, respectively.

Therefore, $R^{I}_{ij}$, which shows the force for member $I$ due to a unit force in $x_{i}$, can be obtained.

Appendix B: A Case Study of a Honeycomb Lattice Truss System

To demonstrate the accuracy of the singum model, a case study is presented for a honeycomb lattice with harmonic potential or linear spring bonds between neighboring nodes. A MATLAB program is developed to simulate the elastic behavior of the lattice for verification of the formulation in Sect. 3, which can be straightforwardly extended to other lattices. In the program, an array of nodes are automatically generated based on the given lattice with the number of unit cells in $X$ and $Y$ directions. The $X-Y$ coordinate origin is set up on the center of the lattice. When the lattice is undergone a deformation, the new coordinate $x-y$ shares the same origin but the coordinates of the nodes change. A list of neighbors for each node is detected with the corresponding bonds and saved for the force computation step. The mid-point of each bond is also important as the potential cutting point of singum surface.

Figure 6 schematically illustrates the honeycomb lattices. The boundary points can be selected by two cases. Figure 6(a) uses the mid-points of the bonds for the boundary; whereas Fig. 6(b) exhibits the boundary on the nodes. If the lattice approaches the infinite domain, the boundary selection produces negligible effects to the effective mechanical behavior. However, when finite unit cells are used, they may produce big differences with different convergence rates to the solution. Due to the periodicity of the singum, Fig. 6(a) provides a better performance and will be used in the study. The algorithm is structured and implemented as follows:

• 1. Initialize the simulation box by periodically extending the singum along $x$ and $y$ direction with $N_{x}$ and $N_{y}$ replications. The 4 sides of the box are made of loading bars as a boundary layer. The node on the boundary is connected to the loading bars by hinges. The initial bond length and force are at $r=2l_{p}^{0}$ and $\textbf{F}^{I}=0$, respectively, so $\lambda =1$.

• 2. Given a testing mode, such as tension or shearing, apply the corresponding uniform singum strain variation $\delta d_{ij}=10^{-6}$ to the simulation box according to the Cauchy–Born (CB) rule. The new positions of all mid-points are updated through an affine transformation with $\delta d_{ij}$.

• 3. Calculate the coordinate of each node from the mid-points by the equilibrium of the node. The length change of each fiber and $\lambda$, and use the potential function $V(\lambda )$ to calculate the equilibrium bond forces. For each loading bar, collect all bond forces and calculate the effective stress vector on each bar in the deformed configuration. Using the effective stress vector on the 4 edges, one can obtain the stress variation caused by $\delta d_{ij}=10^{-6}$ at $\lambda =1$, and thus calculate the elastic tensor.

• 4. For any normalized fiber length or stretch ratio, namely $\lambda ^{i}$, the coordinate of each node and the force in each bond can be calculated with the harmonic potential $V$ in Eq. (5). Repeat Steps 2 and 3 to calculate the elastic tensor at $\lambda ^{i}$. Therefore, the relation of elastic tensor and $\lambda$ can be calculated.

Fig. 6

The honeycomb lattice truss system with boundary: (a) at the mid-points of the bonds and (b) at the nodes

Although the above calculation can guarantee the equilibrium of the internal nodes immediately with the periodic microstructure during the deformation, the result of the elastic tensor is affected by the cut-off for the boundary length but will quickly converge with the number of unit cell $N_{x}$ and $N_{y}$.

On the contrary, if Fig. 6(b) is used directly to generate the new positions of nodes following the C-B rule, the equilibrium positions of the internal nodes will not be periodic anymore, and a very large lattice is required to reach the convergent solution, which is demonstrated in Fig. 7(a).

Fig. 7

Numerical simulation results of elastic tensors changing with different simulation scales and different boundary construction types(a) the comparison of the numerical simulation and the prediction of elastic tensors with different $\lambda$ (b)

As a numerical example, consider a lattice with $k=1000N/m$, $l_{p}^{0}=1mm$, and $\lambda =1$. Equation (33) provides $C_{1111}=C_{2211}=288.6751$N/m and $C_{1212}=0$N/m respectively. 7(a) shows the prediction of $C_{1111}$ and $C_{1122}$ changing with the increase of $N_{x}$ or $N_{y}$ (here we use $N_{x}=N_{y}$). Obviously, for type A lattice in Fig. 6(a), the simulation converges quickly with a few nodes. For instance, with only 184 nodes in total, the simulation results gave $C_{1111}=288.6520$N/m, compared to the prediction value from the singum model, the difference was only about 0.008%. However, for type B lattice in Fig. 6(b), the simulation converges much slower and requires more nodes to dilute the boundary effects as the displacement of nodes on top or bottom boundary does not follow the periodic distribution.

Using the type A lattice, we change the stretch ratio $\lambda$ between 1.0 to 1.2, and show the comparison of the numerical simulation and the prediction of Eq. (33) in 7(b), which exhibit excellent agreement. Indeed, the singum model provides the exact solution for the lattice with the potential between short-range bonds.