# Homogenized elasticity and domain of linear elasticity of 2D archituctured materials

Valentin Jeanneau, C. Combescure, M.L.M. François

► To cite this version:

Valentin Jeanneau, C. Combescure, M.L.M. François. Homogenized elasticity and domain of linear elasticity of 2D archituctured materials. International Journal of Solids and Structures, 2023, 269, pp.112185. 10.1016/j.ijsolstr.2023.112185. hal-04173348

HAL Id: hal-04173348

https://hal.science/hal-04173348

Submitted on 28 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Homogenized elasticity and domain of linear elasticity of 2D architecture materials

V. Jeanneau<sup>a*</sup>, C. Combescure<sup>b,c</sup>, M. Francois<sup>a</sup>

## Abstract

This article presents a generic method to determine analytically the equivalent homogeneous linear elastic behavior of any 2D architecture material and to determine its macroscopic linear elasticity domain. The method is applied throughout the paper in the equilateral triangular case, as an example. Additional examples of square, isosceles right-angled triangular and hexagonal architecture materials are presented as a conclusion. The proposed method combines a periodic homogenization for the equivalent homogeneous linear elastic behavior and three criteria to establish the macroscopic linear elasticity domain. The first criterion corresponds to the limit strength of the bars. The second one, which corresponds to the apparition of a periodic buckling, involving a finite number of cells, is studied on a single primitive unit cell using Bloch wave analysis. The buckling modes that appear under equibiaxial compression for the equilateral triangle architecture material are similar to experimental results reported in the literature. The third criterion is the aperiodic buckling, or macroscopic instability, which corresponds to the apparition of shear bands or a macroscopic buckling. The resultant macroscopic linear elasticity domain resulting from these three criteria is presented, for the first time, in the case of the equilateral triangle, isosceles right-angled triangle architecture material and also for two other geometries, square and hexagonal architecture material. It is shown that the theoretical properties of these architecture materials are, in an Ashby chart, in the domain of light but relatively strong materials, a domain in which man-made materials are currently poorly represented.

**keywords:** Architecture materials, periodic homogenization, linear elasticity, stability, limit surface.

<footnote><sup>a</sup> Institut de Recherche en Génie Civil et Mécanique (UMR CNRS 6183), Nantes Université, 2 rue de la Houssinière, 44322 Nantes Cedex 3, France;</footnote>

<footnote><sup>b</sup> Académie militaire de Saint-Cyr Coëtquidan, CReC, 56380 Guer, France;</footnote>

<footnote><sup>c</sup> Univ. Bretagne Sud, UMR CNRS 6027, IRDL, F-56100 Lorient, France;</footnote>

<footnote>* Corresponding author, valentin.jeanneau@univ-nantes.fr</footnote>
Introduction

Architectured structures have been used for a while to create lightweight and high-
performance structures (for example, the Eiffel Tower). With conventional manufacturing
processes, these architectured structures required the assembly of numerous beams and joints.
This complexity limited their application to large civil engineering structures. Nowadays, ad-
ditive manufacturing, possibly associated with older processes such as lost wax casting, allows
to manufacture these structures at millimeter scale. Manufacturing processes of architectured
material at even smaller scale (micrometer or nanometer) are under development [9, 40].

Architectured materials have three different description scales [8] (Figure 1). The first one called macroscopic is at the scale of the overall structure and considers the architectured material as a continuous and homogeneous material. The second one, called microscopic, is at the scale of the architectured material’s constitutive material, meaning the constitutive material of the beams. The last one, called mesoscopic, is an intermediate scale between the macroscopic and the microscopic and describes the geometrical arrangement of the beams constituting the architectured material. The separation of scales between the mesoscopic scale of the elementary cells and the macroscopic scale of the global structure allows us to speak of architectured materials.

<figure><img src="image_1.png" /><figcaption>Figure 1: The three different scales of description of a architectured material.</figcaption></figure>

The complexity of the fabrication is not the only obstacle to the generalization of the
use of architectured materials. The design computation of a structure made of architectured
material is also a problem because the finite element softwares currently used in engineering
are hampered by their geometrical complexity (very large number of finite elements and
important calculation time). Fortunately, the separation of scales allows the use of the
concept of equivalent homogeneous materials.

Chosen by the engineer, the geometry of the cells can lead to very different macroscopic behaviors such as high specific stiffness, energy absorption, auxetic or intelligent materials and
so on. In general, if the beam network behavior mainly involves bending, the architectured material will behave in a way close to a mechanism<fnref n="1" /> and will have great capacities of deformation leading to great energy absorption up to contact. The 2D architectured materials with square or hexagonal cells are of this type. On the contrary, if the beam network behavior mainly involves tension, the architectured material will have low deformation capacities but a high specific stiffness and a high specific strength [19, 11], qualities that are desired for structural materials. The 2D architectured materials with triangular cells are of the latter type.

The behavior architectured materials is described, a minima, by an elasticity tensor and a linear elastic behavior domain in the stress or strain space. In the case of 2D triangular architectured material, the elasticity tensor, as well as the analysis of its symmetry group, was obtained analytically by [16]. However, very few papers present the linear elastic behavior domain of such materials and if they do, they mostly rely on a strong hypothesis concerning the number of cells to be considered when considering buckling [20]. This article will show that this hypothesis is often too strong.

The linear elasticity domain of architectured materials has to be established considering more than just the yield strength of the constitutive beams. Indeed, the yield strength of classical solid materials corresponds to a degradation of the material either by plasticity, damage or cracking and are generally at least as strong in compression as in tension. On the contrary, architectured materials show a high sensitivity to buckling under compressive states. This loss of linearity is due to instabilities at different scales as shown in numerous experiments [30, 31, 32]. Results from [37] underline that these non linear elastic elastic non-linearity effects, associated to buckling, can occur while the constitutive material remains in its linear elastic domain. Some of the linear behavior domains of various architectured materials have already been established, considering (i) periodic buckling [5, 12, 27, 28, 36, 37] (buckling of a finite number of cells), (ii) aperiodic buckling (buckling of an infinite number of cells, i.e. at the macroscopic scale) [5, 12], or (iii) the constitutive material non linearity [1, 41]. The intersection of these three domains defines the domain of linear elasticity of the equivalent macroscopic homogeneous material. This macroscopic linear elasticity domain, together with the elasticity tensor are both required by engineers to design structures (global shape of the object) and to locally optimize the geometry of the architectured material's cells.

In this article, a combined method of determination of these three surfaces and their intersection is presented. This method uses the Bloch-wave analysis which does not rely on

<footnote n="1">¹is called mechanism a system that can deform without additional external energy</footnote>
any hypothesis concerning the number of cells to consider when studying periodic buckling. It is the first time such a combined method is proposed for 2D architectured material. This method is valid for any type of 2D shape of the unit cell and for any cross section and limit stress of the three constitutive beams. This method is also valid, by extension, if the mesostructure varies slowly with respect to the gradients of the stress field at the macroscopic scale.

An example of an equilateral architectured material with three identical beams is given, for clarity, at any steps of the article. In this example the beam half-length is $L = 10$ mm, the beam width is $b = 1$ mm, the thickness is $h = 1$ mm and the constitutive material is supposed to be an isotropic high grade aluminium of Young's modulus $E = 72$ GPa and limit strength $\sigma_e = 500$ MPa in both tension and compression.

Section 1 of this article is devoted to the linear elasticity of this architectural material (at the macroscopic scale), obtained by a periodic homogenization method. In section 2, the limit surface associated with the yield strength of the constituent material is determined. Section 3 is devoted to the determination of the surface corresponding to the appearance of a periodic buckling, on a finite number of cells. The surface associated with the initiation of aperiodic buckling, over an infinite number of cells, is studied in section 4. In section 5, the surface resulting from the intersection of these three domains is presented, i.e. the boundary surface limiting the linear elastic behavior of the structured material. These surfaces are also presented in the case of an isosceles, right-angled triangle, square and hexagonal architectured materials. Finally, Section 6 draws some conclusions and perspectives to this study.

# 1. Homogeneous linear elastic response

The two-dimensional architectured material considered in this study is composed of a triangular cell beam assembly. This structure is considered perfectly periodic and with infinite extent. The triangle is of any kind and the three composing beams can have different cross-sections and mechanical properties. It has two translational symmetries of vectors {$t_1, t_2$} and is represented by a rhombohedral unit cell shown in Figure 2. The unit cell $C$ presented in Figure 2 has been chosen so that beams are cut in half in their lengths. It is composed of six half beams linking the internal node 1 with the boundary nodes [2, 7].

In this section, for the homogenization, we consider Euler-Bernoulli beams undergoing small displacements and small deformations. Given a half beam, we denote the local abscissa by $x \in [0, L]$ where $L$ is the length of the considered half beam. Additionally, the axial displacement is denoted by $u(x)$, the transverse displacement by $v(x)$, the rotation of the section by $\theta(x) = v'(x)$. Assuming a linear interpolation of the axial displacement and a cubic interpolation for the transverse displacement, the standard shape functions of the
<figure><img src="image_2.png" /><figcaption>Figure 2: 2D architecture material (grey lines) and unit cell $\mathcal{C}$ (dashed red lines).</figcaption></figure>

finite element theory are used. These shape functions relate the kinematic of the beam end nodes $q^i = [u^i v^i L\theta^i]$, i.e. the degrees of freedom (DOF), to the displacement of any point at abscissa $x$:

$$ \underbrace{\begin{bmatrix} u(x) \\ v(x) \\ L\theta(x) \end{bmatrix}}_{[U(x)]} = \underbrace{\begin{bmatrix} 1-\xi & 0 & 0 & \xi & 0 & 0 \\ 0 & 1-3\xi^2+2\xi^3 & \xi(1-\xi)^2 & 0 & 3\xi^2-2\xi^3 & \xi^2(\xi-1) \\ 0 & -6\xi(1-\xi) & 1-4\xi+3\xi^2 & 0 & 6\xi(1-\xi) & \xi(3\xi-2) \end{bmatrix}}_{[\mathcal{N}]} \underbrace{\begin{bmatrix} u^i \\ v^i \\ L\theta^i \\ u^j \\ v^j \\ L\theta^j \end{bmatrix}}_{[q^i q^j]^T}, \quad (1) $$

where $\xi = x/L$. Here, the degree of freedom for the rotation of the section is chosen as $L\theta(x)$ in order to ensure a good conditioning of the matrix $\mathcal{N}$. For a beam, the elastic energy is given the formula:

$$ W = \frac{1}{2} \int_{0}^{L} [ESu'(x)^2 + EIv''(x)^2] dx. \quad (2) $$

Using Equation (1), this energy is:
$$W = \frac{1}{2} [q^i q^j] K [q^i q^j]^T, \qquad (3)$$

where K is given by:

$$K = \frac{EI}{L^3} \begin{bmatrix} \Gamma & 0 & 0 & -\Gamma & 0 & 0 \\ 0 & 12 & 6 & 0 & -12 & 6 \\ 0 & 6 & 4 & 0 & -6 & 2 \\ -\Gamma & 0 & 0 & \Gamma & 0 & 0 \\ 0 & -12 & -6 & 0 & 12 & -6 \\ 0 & 6 & 2 & 0 & -6 & 4 \end{bmatrix}, \qquad (4)$$

in which $\Gamma = SL^2/I$, $S$ is the cross section, $I$ is the second moment of area, and $E$ is the Young's modulus.

The global energy $W$ of the unit cell $C$ is the sum of its six half beam elastic energies for which $(E, I, L, S)$ can be different for each beam (2, 5), (3, 6) and (4, 7) as pictured on Figure 2.

At the macro-scale, the position of a current point $X$ is referred to the orthonormal frame $(e_1, e_2)$ (see Figure 2). According to the Cauchy-Born hypothesis (also called Periodic Boundary Conditions) [5, 6, 17], the displacement field at the macro scale $U(X)$ is the sum of the homogenous field $\hat{U}(X)$ and some periodic fluctuation $\tilde{U}(X)$:

$$U(X) = \hat{U}(X) + \tilde{U}(X). \qquad (5)$$

The homogenous field $\hat{U}(X)$ results of some homogenous Green-Lagrange strain tensor $E$ such as $\hat{U} = E \cdot X$. The period of the periodic fluctuation $\tilde{U}(X)$ corresponds to the one of the unit cell $C$ (Figure 3).

<figure><img src="image_3.png" /><figcaption>Figure 3: Schematic composition of a displacement field with Cauchy-Born hypothesis.</figcaption></figure>

The homogenization process requires to compute the unknown periodic displacement field $\tilde{U}$ for all case of imposed homogenous strain tensor $E$. Such field expresses as:
$$ \tilde{U}(\mathbf{X} + p^1\mathbf{t}_1 + p^2\mathbf{t}_2) = \tilde{U}(\mathbf{X}), \quad \forall [p^1, p^2] \in \mathbb{Z}^2, \quad \forall \mathbf{X} \in \mathcal{C}, \qquad (6) $$

where $\mathbf{t}_1$ and $\mathbf{t}_2$ are the so called basis vectors representing the minimal translational symmetry of the unit cell $\mathcal{C}$ (Figure 2). Thanks to Equation (1), the displacement field $U(\mathbf{X})$ is only needed at the nodes $i \in \{1, \dots, 7\}$. As a consequence, the Cauchy-Born hypothesis resumes as:

$$ q^i = \hat{q}^i + \tilde{q}^i, \qquad (7) $$

where $\hat{q}^i$ is the macroscopic homogenous imposed DOF at node $i$ with:

$$ \hat{q}^i = \left[ E_{11}X_1^i + E_{12}X_2^i, \ E_{21}X_1^i + E_{22}X_2^i, \ 0 \right], \qquad (8) $$

where last component $\theta^i = 0$ because a homogenous strain field implies no rotation (the nodal rotations are then only created by the periodic fluctuation term), and $\tilde{q}^i$ is the unknown periodic DOF such as:

$$ \forall [p^1, p^2] \in \mathbb{Z}^2, \quad \mathbf{X}^j = \mathbf{X}^i + p^1\mathbf{t}_1 + p^2\mathbf{t}_2 \quad \Rightarrow \quad \tilde{q}^j = \tilde{q}^i. \qquad (9) $$

The seven displacements are gathered in corresponding vectors:

$$ \begin{aligned} \hat{Q} &= [\hat{q}^1, \dots, \hat{q}^7], \\ \tilde{Q} &= [\tilde{q}^1, \dots, \tilde{q}^7]. \end{aligned} \qquad (10) $$

For the periodic part, the boundary nodes (2, 5), (3, 6) and (4, 7) are pairwise similar. This is expressed by:

$$ \begin{aligned} \tilde{Q} &= \tilde{Q}^* Z_0, \\ Z_0 &= \begin{bmatrix} \mathbb{I} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \mathbb{I} & 0 & 0 & \mathbb{I} & 0 & 0 \\ 0 & 0 & \mathbb{I} & 0 & 0 & \mathbb{I} & 0 \\ 0 & 0 & 0 & \mathbb{I} & 0 & 0 & \mathbb{I} \end{bmatrix}, \end{aligned} \qquad (11) $$

where $\mathbb{I}$ stands for the 3 × 3 identity matrix and $\tilde{Q}^* = [\tilde{q}^1, \dots, \tilde{q}^4]$ are the independent DOF. Finally the total nodal displacement is given by:

$$ Q = \hat{Q} + \tilde{Q}^* Z_0. \qquad (12) $$
Gathered with previous results, by assembling the matrices K in a global matrix K, this equation gives access to the global elastic energy:

$$W = \frac{1}{2} Q K Q^T. \qquad (13)$$

One could be surprised that no forces appear in the system. This is due to the infinity of the whole structure as it occurs the infinite spring system [7]. As a consequence, the total potential energy is reduced to the internal energy W. The equilibrium of the unit cell C can be found by the minimization of W:

$$\frac{\partial W}{\partial Q} = 0, \qquad (14)$$

leading to:

$$K Q^T = 0. \qquad (15)$$

Using Equation(12), the equilibrium equation becomes:

$$Z_0 K Z_0^T \tilde{Q}^{*T} = -Z_0 K \hat{Q}^T, \qquad (16)$$

which resembles to a classical finite element formulation in which the second member plays the role of a force term.

The rank of matrix $Z_0 K Z_0^T$ is smaller than the number of unknown displacements in the system, meaning that the solution is not unique due to the persistence of rigid body translations. Indeed, periodic boundary conditions (11) will prevent rigid body rotations but rigid body translations are also solutions to the problem. To ensure the uniqueness of the solution, these rigid body translations are suppressed by imposing a null mean displacement of the independent nodes:

$$\sum_{i=1}^{4} u^i = 0, \quad \sum_{i=1}^{4} v^i = 0, \qquad (17)$$

which can be written as:

$$ \begin{aligned} H Q^{*T} &= 0, \\ H &= \begin{bmatrix} 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \end{bmatrix}. \end{aligned} \qquad (18) $$
The minimum of the energy W (13), under the condition (18), is obtained by minimising the following Lagrange expression $QKQ^T/2 + \Lambda H Q^{*T}$ which, from Equation (16), resumes as:

$$ \begin{bmatrix} Z_0 K Z_0^T & H^T \\ H & 0 \end{bmatrix} \begin{bmatrix} \tilde{Q}^{*T} \\ \Lambda \end{bmatrix} = \begin{bmatrix} -Z_0 K_t \hat{Q}^T \\ 0 \end{bmatrix}. \qquad (19) $$

Solving this problem gives access to the nodal displacements $Q$.

Thanks to this nodal displacement $Q$ we can access to the homogenized elastic behavior of the macroscopic material. The latter is described by the fourth order stiffness tensor $C$ which relates the imposed strain $E$ to the macroscopic Cauchy stress $\Sigma$.

$$ \Sigma = C : E. \qquad (20) $$

The Hill-Mandel lemma [21] sets the equality between the strain energy at both the macroscopic and mesoscopic scales:

$$ \frac{\mathcal{A}}{2} E : C : E = W, \qquad (21) $$

where $\mathcal{A}$ is the area of the unit cell $C$. As a consequence, using Equation (13), the stiffness tensor is given by:

$$ C = \frac{1}{2\mathcal{A}} \frac{\partial^2}{\partial E^2} QKQ^T, \qquad (22) $$

where $Q$ depends upon the imposed deformation $E$ according to Equation (8).

For the chosen example of equilateral cells with identical beams (see the introduction), this method leads to the following analytical result which is in agreement with [16]:

$$ C = \frac{\sqrt{3}bhE}{32L^3} \begin{bmatrix} 12L^2 + b^2 & 4L^2 - b^2 & 0 \\ 4L^2 - b^2 & 12L^2 + b^2 & 0 \\ 0 & 0 & 2(4L^2 + b^2) \end{bmatrix}, \qquad (23) $$

once expressed in a tensor basis [3] for which the component $\epsilon_{33}$ is $2C_{1212}$. As expected for such $D_6$ symmetry group and according to [15, 16], this tensor is isotropic.

## 2. Material non-linearity limit surface

The first criterion taken into account is the limit of elasticity of the constitutive material of the beams. According to the beam theory and supposing the same strength in tension and compression, this criterion is:

$$ \max_x (|\sigma(x)|) < \sigma_e, \qquad (24) $$
where

$$ \sigma(x) = \frac{N(x)}{S} + \frac{M(x)b}{2I}, \qquad (25) $$

with $b$ is the width of the beam, $S$ the cross-section area and $I$ the second moment of area and where $N$ and $M$ are the normal force and the bending moment. The latter are computed from the nodal displacements:

$$ \begin{gathered} \begin{bmatrix} N(\xi) \\ T(\xi) \\ M(\xi)/L \end{bmatrix} = \frac{EI}{L^3} [\mathcal{N}'(\xi)][q^i \ q^j]^T, \\ [\mathcal{N}'(\xi)] = \begin{bmatrix} -\Gamma & 0 & 0 & \Gamma & 0 & 0 \\ 0 & -12 & -6 & 0 & 12 & -6 \\ 0 & -6 + 12\xi & -4 + 6\xi & 0 & 6 - 12\xi & -2 + 6\xi \end{bmatrix}, \end{gathered} \qquad (26) $$

where, according to the previous section, $q^i$ depends upon the imposed homogenous strain tensor $E$ which can be decomposed in a loading factor $\lambda \in \mathbb{R}^+$ and a normalized tensor $e$ as:

$$ E = \lambda e. \qquad (27) $$

According to [37], $e$ is parametrized by two angles $\varphi \in [-\pi, \pi]$ and $\alpha \in [0, \pi]$:

$$ e = \begin{bmatrix} \sin(\alpha) \cos(\varphi) & \cos(\alpha) \\ \cos(\alpha) & \sin(\alpha) \sin(\varphi) \end{bmatrix}, \qquad (28) $$

and the value of $\lambda$ can be computed for each retained set of angles.

In our example of application, the high slenderness ($L/b \gg 1$ and the triangular shape leads to $|M|b/2I \ll |N|/S$ which is called stretch dominated [11]. This leads to a simple and analytical result of a limit surface composed of 6 planes parallel two by two due to the
symmetry of the unit cell and the symmetry of the constitutive material behavior:

$$
\begin{align}
\lambda_M e_{11} + \frac{\sigma_e}{E} &= 0, \nonumber \\
\frac{\lambda_M}{4} (e_{11} + 2\sqrt{3}e_{12} + 3e_{22}) + \frac{\sigma_e}{E} &= 0, \nonumber \\
\frac{\lambda_M}{4} (e_{11} - 2\sqrt{3}e_{12} + 3e_{22}) + \frac{\sigma_e}{E} &= 0, \nonumber \\
\lambda_M e_{11} - \frac{\sigma_e}{E} &= 0, \nonumber \\
\frac{\lambda_M}{4} (e_{11} + 2\sqrt{3}e_{12} + 3e_{22}) - \frac{\sigma_e}{E} &= 0, \nonumber \\
\frac{\lambda_M}{4} (e_{11} - 2\sqrt{3}e_{12} + 3e_{22}) - \frac{\sigma_e}{E} &= 0, \tag{29}
\end{align}
$$

where $\lambda_M$ is the critical load factor for the Material limit of elasticity criterion. In the strain space $(E_{11}, E_{22}, E_{12})$, these six planes define the domain of constitutive material linear elasticity (Figure 4). As supposed for the constitutive material, this surface is symmetric in

<figure><img src="image_4.png" /><figcaption>Figure 4: Limit surface associated with the limit of elasticity of the constitutive material.</figcaption></figure>

tension and compression. However, the section of the domain (Figure 4b) shows the absence of symmetry by a rotation of $\pi/2$ which corresponds to the permutation $1 \leftrightarrow 2$. This is due to the symmetry group of this particular cell which does not include such symmetry.

### 3. Periodic buckling limit surface

When subjected to compression, a single slender beam may buckle before the elasticity limit of the constitutive material is reached. This effect, which corresponds to an instability,
has been first studied by Euler [14]. In architectured materials, buckling concerns many beams and many cells and we study hereafter the periodic case, the periodicity of which can differ from that of the unit cell C. Such buckling cases have been shown by experiments in basic compressive states [22, 29, 32, 35].

We will study these instabilities spreading on possibly larger periodicity using Bloch waves analysis [4, 18]:

$$ \tilde{U}(\mathbf{X}) = \tilde{V}(\mathbf{X}) \exp (2i\pi\eta \cdot \mathbf{X}), \qquad (30) $$

with $\tilde{V}$ a periodic displacement field on the unit cell C as defined in Equation (6) and $\eta = \eta^1\mathbf{t}_1 + \eta^2\mathbf{t}_2 = \eta_1\mathbf{t}^1 + \eta_2\mathbf{t}^2$ the wave-vector expressed in both the direct basis $\{\mathbf{t}_1, \mathbf{t}_2\}$ and the dual basis $\{\mathbf{t}^1, \mathbf{t}^2\}$, such as $\mathbf{t}_i \cdot \mathbf{t}^j = \delta_i^j$. As a consequence, the former periodicity condition of Equation (6) is replaced by the Bloch periodicity condition:

$$ \begin{aligned} \tilde{U}(\mathbf{X} + \mathbf{t}) &= \tilde{V}(\mathbf{X} + \mathbf{t}) \exp (2i\pi\eta \cdot (\mathbf{X} + \mathbf{t})), \\ \tilde{U}(\mathbf{X} + \mathbf{t}) &= \tilde{U}(\mathbf{X}) \exp (2i\pi\eta \cdot \mathbf{t}) \quad \text{if} \quad (p^1, p^2) \in \mathbb{Z}^2, \\ &= \tilde{U}(\mathbf{X}) \quad \text{if} \quad (\eta \cdot \mathbf{t}, p^1, p^2) \in \mathbb{Z}^3. \end{aligned} \qquad (31) $$

One can thus see that the actual periodicity vector $\mathbf{t}$ of function $\tilde{U}$ is linked to $\eta$ and cannot be smaller than the periodicity of the unit cell C because it is the periodicity of $\tilde{V}$.

Additionally, the complex exponential function appearing in Equation (31) is a periodic function in $\eta$ with period $[\eta_1, \eta_2] = [1, 1]$ (the covariant coordinates), allowing us to reduce the range of variation for covariant coordinates $[\eta_1, \eta_2]$ within $[0, 1] \times [0, 1]$. This periodic box in $\eta$ is conventionally called the first Brillouin zone in crystallography.

As a consequence, Equation (31) allows us to set up a larger periodicity than that of the initial unit cell C thus enabling the possibility to extend the stability analysis to the case of period multiplying instabilities that have been shown to happen in periodic architectured materials.

In the context of this article using a finite element formulation, Equation (31.2) must be applied at the boundary nodes of our unit cell C:
$$ \begin{gathered} \tilde{q}^7 = z_1 \tilde{q}^4, \quad \tilde{q}^6 = z_2 \tilde{q}^3, \quad \tilde{q}^5 = z_1 z_2 \tilde{q}^2, \\ \text{with } z_i = \exp(2i\pi\eta_i \cdot t_i) = \exp(2i\pi\eta_i). \end{gathered} \tag{32} $$

In the same way as for the Equation (11), the displacement field of all the nodes in a unit cell can be written as a function of the independent node displacements [5]:

$$ \begin{gathered} \tilde{Q} = \tilde{Q}^* Z_{\eta}, \\ Z_{\eta} = \begin{bmatrix} \mathbb{I} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \mathbb{I} & 0 & 0 & z_1 z_2 \mathbb{I} & 0 & 0 \\ 0 & 0 & \mathbb{I} & 0 & 0 & z_2 \mathbb{I} & 0 \\ 0 & 0 & 0 & \mathbb{I} & 0 & 0 & z_1 \mathbb{I} \end{bmatrix}. \end{gathered} \tag{33} $$

The purpose of these Bloch periodicity conditions is to find, for each loading case, the critical load and the wave-vector of the patterns created by the instability.

Modelling the buckling requires a large displacements formulation [2]. As a consequence the previous stiffness matrix of Equation (4) is replaced by the tangent stiffness matrix which involves the geometrical non linear effect:

$$ K_t = \frac{EI}{L^3} \begin{bmatrix} \Gamma & 0 & 0 & -\Gamma & 0 & 0 \\ 0 & 12+36f & 6+3f & 0 & -12-36f & 6+3f \\ 0 & 6+3f & 4+4f & 0 & -6-3f & 2-f \\ -\Gamma & 0 & 0 & \Gamma & 0 & 0 \\ 0 & -12-36f & -6-3f & 0 & 12+36f & -6-3f \\ 0 & 6+3f & 2-f & 0 & -6-3f & 4+4f \end{bmatrix}, \tag{34} $$

where

$$ f = \frac{NL^2}{30EI}, $$

where N is the axial force in the beam. Note that, at the onset of instability, the internal forces stated in Equation (2) are still relevant. Thus, the axial force N depends linearly of the loading factor λ of Equation (27). As a consequence, the tangent stiffness matrix depends also upon λ : K_t(λ).
Instability occurs as soon as an eigenvalue of the assembled tangent stiffness matrix $K_t$ is null. Using the Bloch periodicity conditions, the periodic buckling criterion writes [5, 18, 37]:

$$ \det[Z_{\eta}K_t(\lambda)Z_{\eta}^H] = 0, \qquad (35) $$

where $Z_{\eta}^H$ is the conjugate transpose matrix of $Z_{\eta}$ and $K_t(0)$ is positive definite. A point of the researched periodic buckling surface is defined as the apparition of the first instability i.e. for the smallest $\lambda$ and also depends upon the wave-vector $\eta$:

$$ (\lambda_P, \eta_P) = \operatorname*{argmin}_{\eta} (\lambda(\eta) | \det[Z_{\eta}K_t(\lambda)Z_{\eta}^H] = 0), \qquad (36) $$

where index $P$ stands for Periodic buckling, $\eta_P$ provides the periodicity of the critical buckling mode and $\lambda_P$ is the critical loading factor which gives access to the periodic buckling surface. Even if it is not the objective of this article, the values of $(\lambda_P, \eta_P)$ give access to $Z_{\eta}K_t(\lambda)Z_{\eta}^H$ thus to the eigenvectors $\tilde{Q}^*$ associated with the null eigenvalues. From them one can derive the whole displacement field $U$ of the buckled mode using Equations (33) and (1) thus draw the whole pattern created by the instability.

Due to the non linear nature of the tangent matrix $K_t$ and of Equation (36), we were unable to get general analytical results. Thus the points of the periodic buckling surface were computed in the specific case of our example, by using a finite number of loading directions $e$. The discrete values of $(\alpha, \varphi)$ defined in Equation (28) were chosen so that the density of the points remains constant on the unit sphere [24]. The obtained periodic buckling surface is shown on Figure (5). It is interesting to note that periodic instabilities do not appear under pure tension ($E_{11} > 0$ and $E_{22} > 0$).

In order to compare our predictions with both numerical and experimental results obtained by [22], we have studied the buckling modes in the case of the equi-biaxial compression $E_{11} = E_{22}$ and $E_{12} = 0$ which corresponds to $\varphi = -3\pi/4$ and $\alpha = \pi/2$. The solution of Equation (36) is $\lambda_P = 0.00445$, $\eta_P \simeq (1/3, 1/3)$. In this case the matrix $Z_{\eta}K_t(\lambda)Z_{\eta}^H$ has two null eigenvalues, leading to a two-dimensional eigenspace. Using the method of [10], it has been possible to project this eigenvectors $\tilde{Q}$ on a peculiar basis defined by elementary symmetry groups of the deformed states. As expected, this basis, constituted of two buckling modes, presented in Figure 6, is in concordance with the results obtained by [22].

Interestingly, most of tested loading cases lead to a critical buckling mode over a very large number of cells (the covariant components $\eta_1, \eta_2$ of $\eta_P$ are not simple fractions). For
<figure><img src="image_5.png" /><figcaption>Figure 5: Limit surface associated with the periodic buckling.</figcaption></figure>

example, for uniaxial loading along **e**₁ axis, the critical wave-vector **η**<sub>P</sub> ≈ (0, 71039; 0, 71031) leads to a buckling mode involving 547 and 434 cells along **t**₁ and **t**₂, respectively.

This possibility was excluded in article [20] since they only considered buckling modes
including only one or two cells. For a comparison, our method predicts a critical buckling
deformation of 0.00445, where the method of [20] predicts 0.00489 for the equibiaxial loading.
Their method overestimates by almost 10% the critical buckling deformation and, as such,
is overly conservative. Our more complete Bloch-wave analysis shows that any method
relying on a strong hypothesis that predetermines in advance the number of cells to be
considered for periodic buckling presents a risk of overestimating the critical buckling limit
of 2D triangular architectured materials. As a consequence, even though the limit surface of
2D equilateral triangular architectured material seems to have been already presented in [20],
we are doubtful on the validity of this proposed limit surface. Our proposed limit surface
appears to be the one to be considered for those materials, since it appropriately captures
the critical buckling mode without predetermining the number of buckling cells involved.

4. Aperiodic buckling limit surface

Aperiodic buckling is associated with two kind of macroscopic instabilities [34]. The first
one corresponds to shear banding and the second one corresponds to a global buckling of the
whole structure (Figure 7). In both cases, an infinite number of cells is involved but neither
of them presents periodic pattern buckling. As a consequence, this type of instabilities is
<figure><img src="image_6.png" /><figcaption>Figure 6: Buckling modes of the regular triangular cell architecture under equi-biaxial compression. In dashed red line: contour of the 3x3 periodic pattern. In solid red line: contour of specific mode shape identified by [22].</figcaption></figure>

referred to as aperiodic buckling. A post-bifurcation analysis would be necessary to determine
which of the two modes appears but this is out of the scope of this article.

<figure><img src="image_7.png" /><figcaption>Figure 7: Aperiodic instabilities: (left) shear band along a band of normal **k**, (right) global buckling of the structure.</figcaption></figure>

The aperiodic instability criterion is based on the homogenized continuum tensor $\mathbb{H}$,
calculated the same way as equation (22) but considering the geometrical non-linearities of
the beam i.e. using $K_t$ instead of $K$ and $F$ instead of $E$. Here $F$ refers to the non symmetric deformation gradient tensor.

$$ \mathbb{H} = \frac{1}{2\mathcal{A}} \frac{\partial^2}{\partial \mathbf{F}^2} Q K_t Q^T. \qquad (37) $$

This criterion is the rank-one convexity criterion [13, 22], also called strong ellipticity criterion [38]:

$$ (\mathbf{k} \odot \mathbb{H} \cdot \mathbf{k}) \cdot \mathbf{g} = 0, \qquad (38) $$

where $\mathbf{k}$ and $\mathbf{g}$ are any unitary vectors and the operator $\odot$ is such that $(\mathbf{k} \odot \mathbb{H})_{ikl} = k_j H_{ijkl}$. Since $\mathbb{H}$, defined in Equation (37), have the major symmetry $H_{ijkl} = H_{klij}$, Equation (38) is equivalent to the Mandel-Rice criterion [25, 33]:

$$ \det(\mathbf{k} \odot \mathbb{H} \cdot \mathbf{k}) = 0. \qquad (39) $$

The tensor formed by $\mathbf{k} \odot \mathbb{H} \cdot \mathbf{k}$ is called the acoustic tensor $\mathbf{A}$ and its components are calculated by [24, 26]:

$$ A_{ik} = k_j H_{ijkl} k_l. \qquad (40) $$

In case of a shear band, $\mathbf{k}$ is the normal to the shear band. This criterion is also equivalent to the periodic buckling criterion (36) when $\|\eta\| \to 0$ [13].

Again, a point of the researched aperiodic buckling surface is defined by the apparition of the first instability i.e. for the smallest $\lambda$:

$$ (\lambda_A, \mathbf{k}_A) = \argmin (\lambda(\mathbf{k}) | \det(\mathbf{k} \odot \mathbb{H}(\lambda) \cdot \mathbf{k}) = 0), \qquad (41) $$

where index $A$ stands for Aperiodic buckling and $\lambda_A$ is the critical loading factor which gives access to the aperiodic buckling surface.

In our example of application, this analysis was carried out numerically. The obtained aperiodic buckling surface is presented in Figure 8. We can notice that, again, aperiodic buckling does not occurs when the architectured material is under pure tension.

## 5. Resulting domains of linear elasticity

In general, the domain of linear elasticity is bounded by the previous three surfaces. More precisely, for each loading direction $e$, it is defined by the lowest $\lambda = \min(\lambda_M, \lambda_P, \lambda_A)$.

In this example, the aperiodic buckling surface is always outside of the two others as pictured
<figure><img src="image_8.png" /><figcaption>Figure 8: Limit surface associated with the aperiodic buckling.</figcaption></figure>

on Figure 9. This figure presents the three limit surfaces and the stability domain in the section ($E_{11}$, $E_{22}$) of the deformation space.

<figure><img src="image_9.png" /><figcaption>Figure 9: Comparison of the three limit surfaces in the plane $E_{12} = 0$ for the equilateral triangle. The dots are the calculated points and the line is the interpolation between two dots. The blue line is the material non linearity limit surface, the red line is the periodic buckling limit surface and the green line is the aperiodic buckling limit surface. The area shaded in gray represents the domain of linear elasticity. The considered unit cell for this calculation is represented in the top right corner.</figcaption></figure>
The computed domain of linear elasticity is presented in the 3D whole space of the deformations in Figure 10. This domain can also be represented in a harmonic deviatoric and spheric basis [23] as in Figure 11. In this Figure, components are calculated with :

$$
\begin{align}
E_{d_1} &= \frac{\sqrt{2}}{2} (E_{11} - E_{22}), \nonumber \\
E_{d_2} &= \sqrt{2} E_{12}, \nonumber \\
E_s &= \frac{\sqrt{2}}{2} (E_{11} + E_{22}). \tag{42}
\end{align}
$$

It is interesting to note that this domain presents a symmetry of order 3 around the spheric axis visible on Figure 11. We have no explanation for this observation at that time.

<figure><img src="image_10.png" /><figcaption>Figure 10: Domain of linear elasticity for the equilateral triangle lattice limited by the material non linearity limit surface (■) and the periodic buckling limit surface (■) in the canonic basis.</figcaption></figure>

All the results presented previously can be easily calculated for any triangular unit cell.
For example, for the unit cell in Figure 12, the analytical elasticity tensor $C_{\text{tri2}}$ is:

$$
C_{\text{tri2}} = \frac{bhE}{L^3} \begin{bmatrix} \frac{\sqrt{2}b^2 + 2(4 + \sqrt{2})L^2}{8} & -\frac{b^2 - 2L^2}{4\sqrt{2}} & 0 \\ -\frac{b^2 - 2L^2}{4\sqrt{2}} & \frac{b^2 + 2L^2}{4\sqrt{2}} & 0 \\ 0 & 0 & \frac{\sqrt{2}b^2 + 2(4 + \sqrt{2})L^2}{4 + 8\sqrt{2}} \end{bmatrix}, \quad (43)
$$
<figure><img src="image_11.png" /><figcaption>Figure 11: Isovalues of the limit surface of the linear elasticity domain in the deviatoric deformation space (x1000). Isovalues are taken each 0.0015 along the spheric axis.</figcaption></figure>

which correspond to an orthotropic behavior. The three surfaces corresponding to the three criteria are represented in Figure 12. As for the equilateral triangle, the aperiodic buckling limit surface does not intervene in the construction of the linear elasticity domain. This domain is bounded by the material non linearity limit surface and the periodic buckling limit surface. For this example, thickness of the beams is calculated so that the relative density is the same between the equilateral triangle and the isosceles, right-angled triangle lattice materials.

Our method also allows to easily calculate the elasticity tensor and the domain of linear elasticity of more general 2D periodic architectured materials such as square (Figure 13) and hexagonal (Figure 14) lattice. Constitutive material is supposed to be the same as the equilateral triangle lattice and thickness of the beams are chosen so that the relative density remains the same in all the considered lattice materials. The elasticity tensor calculated for square lattice $C_{\text{Square}}$ shows tetragonal symmetry while $C_{\text{Hexa}}$ exhibits isotropic symmetry.
<figure><img src="image_12.png" /><figcaption>Figure 12: Comparison of the three limit surfaces in the plane $E_{12} = 0$ for the isosceles, right-angled triangle. The dots are the calculated points and the line are the interpolation between two dots. The blue line is the material non linearity limit surface, the red line is the periodic buckling limit surface and the green line is the aperiodic buckling limit surface. The area shaded in gray represents the domain of linear elasticity. The considered unit cell for this calculation is represented in the top right corner. The length L corresponds to the length of the horizontal beam going from the central node to the edge of the unit cell.</figcaption></figure>

$$ \mathcal{C}_{\text{Square}} = \frac{bhE}{2L} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \frac{b^2}{32L^2} \end{bmatrix}, \qquad (44) $$

$$ \mathcal{C}_{\text{Hexa}} = \frac{bhE}{\sqrt{3L(b^2 + 4L^2)}} \begin{bmatrix} 3b^2 + 4L^2 & -b^2 + 4L^2 & 0 \\ -b^2 + 4L^2 & 3b^2 + 4L^2 & 0 \\ 0 & 0 & \frac{b^2}{32L^2} \end{bmatrix}. \qquad (45) $$

In the case of square lattice (Figure 13), the periodic and aperiodic buckling limit surfaces are superimposed. This results highlights the fact that these two criteria are equivalent when $||\eta|| \to 0$ and that the periodic buckling criteria can detect long wavelength (macroscopic) instabilities. It is interesting to note that with the chosen material and geometry, the domain of linear elasticity for the hexagonal lattice (Figure 14) is only bounded by the material non-linearity limit surface. This is explained by the fact that, in order to have a relative density comparable with that of our equilateral triangular architectured material, the hexagonal one has to be composed with three times thicker beams that would not favor buckling. Considering Figures 9, 12, 13, 14, lattice with hexagonal cells has the greatest admissible
<figure><img src="image_13.png" /><figcaption>Figure 13: Comparison of the three limit surfaces in the plane $E_{12} = 0$ for the square lattice. The dots are the calculated points and the line is the interpolation between two dots. The blue line is the material non linearity limit surface, the red line is the periodic buckling limit surface. The green dots is the aperiodic buckling limit surface. The area shaded in gray represents the domain of linear elasticity. The considered unit cell for this calculation is represented in the top right corner. The length $L$ corresponds to the length of the horizontal beam going from the central node to the edge of the unit cell.</figcaption></figure>

deformation in compression because no instability occurs. With the three other considered lattices, instabilities occurs in compression. Among them, the equilateral triangle lattice shows the greatest linear elasticity domain in deformation.

The critical load on the domain of linearity allows to define the maximal and minimal strength of our lattice materials. These points are represented on an Ashby's diagram in order to compare the strength of our lattices with other known materials in Figure 15. Our materials are as resistant as wood along its fibers while being lighter.
<figure><img src="image_14.png" /><figcaption>Figure 14: Comparison of the three limit surfaces in the plane $E_{12} = 0$ for the hexagonal lattice. The dots are the calculated points and the line is the interpolation between two dots. The blue line is the material non linearity limit surface, the red line is the periodic buckling limit surface and the green line is the aperiodic buckling limit surface. The area shaded in gray represents the domain of linear elasticity. The considered unit cell for this calculation is represented in the top right corner. The length $L$ corresponds to the half length of the horizontal beam.</figcaption></figure>

## 6. Conclusion

In this article, the elastic behavior and the limit of linearity is given for 2D architectured materials. The components of the elasticity tensor are calculated analytically and the effective behavior at macroscopic scale is isotropic for equilateral triangular and hexagonal architectured materials, tetragonal for square and orthotropic for isosceles right angled triangular ones. Moreover, three criteria are established in order to build a limit surface of linearity. The first one detects the linearity limit of the constitutive material and is established analytically. The second one is the Floquet-Bloch criterion, that detects the buckling modes of the architectured medium on a finite number of cell. This criterion uses numerical resolution and is thus limited to an integer number of cells. However, it is possible that the very large number of cells actually corresponds to incommensurate wavelength buckling modes. The possibility of incommensurate wavelength buckling modes has been studied in a one-directional case by [39]. The last one, also obtained numerically, is the Rank-One convexity criterion, that deals with instabilities at macroscopic scale that propagates over an infinite number of cell. Once established, these three criteria can be represented as surfaces in the strain space. With the chosen architectured and constitutive material, it appears that the aperiodic buckling
<figure><img src="image_15.png" /><figcaption>Figure 15: Comparison of the strength vs density ratio of bulk aluminum (red dot) and the range of the linear elasticity domain of the equilateral triangle (black segment), the isosceles, right-angled triangle (red segment), square (blue segment) and hexagonal (yellow segment) lattice material with the same relative density in Ashby's diagram (from Ansys Granta Selector with permission). The upper and lower bounds represent the resistance of the material in traction and compression.</figcaption></figure>

criterion does not take part in the generation of the domain of linear behavior, meaning that whatever the macroscopic applied loading, the linearity limit will be reached by the limit of constitutive material linearity or by the periodic buckling. This is a constructive advantage for engineers. When the example 2D equilateral triangular architecture is subjected to pure macroscopic tension, its linearity behavior is limited by the material criterion while in pure macroscopic compression, the linear behavior is limited by mesoscopic instabilities, creating patterns in the mesoscopic structure. In between (mixed traction and compression), the limit of the material is due either to the constitutive material criterion or to the mesoscopic criterion. Thanks to the methodology explained in this article, it is possible to determine all the admissible macroscopic sollicitations ensuring the linear response of any 2D periodic architecture.

This work can be extended by considering different thickness in the beams and other simple geometries for the unit cell in order to change the degree of anisotropy in the Cauchy
elasticity. This will also change the limit surface of linearity providing new admissible
sollicitation space. This calculation can be put in optimization algorithm to get the most
adapted anisotropy while ensuring the linear response of the material and even reach non
accessible areas in Ashby’s diagram.

Acknowledgments

This work is found by the French National Research Agency and is realized as part of
the project ANR-19-CE10-0005, Measurement and Optimization of Periodic Architectured
Materials.
References

References

[1] S. Arabnejad and D. Pasini. Mechanical properties of lattice materials via asymptotic homogenization and comparison with alternative homogenization methods. *International Journal of Mechanical Sciences*, 77:249–262, 2013.

[2] K. Bathe and S. Bolourchi. Large displacement analysis of three-dimensional beam structures. *International Journal for Numerical Methods in Engineering*, 14(7):961–986, 1979.

[3] P. Bechterew. Analytical study of the generalized hooke's law. application of the method of coordinate transformation. *Zh. Russ. Fiz.-Khim. Obshch. Leningrad. Univ., Fizika*, 58(3):415–416, 1926.

[4] K. Bertoldi, M. Boyce, S. Deschanel, S. Prange, and T. Mullin. Mechanics of deformation-triggered pattern transformations and superelastic behavior in periodic elastomeric structures. *Journal of the Mechanics and Physics of Solids*, 56(8):2642–2668, 2008.

[5] G. Bordiga, L. Cabras, A. Piccolroaz, and D. Bigoni. Dynamics of prestressed elastic lattices: Homogenization, instabilities, and strain localization. *Journal of the Mechanics and Physics of Solids*, 146:104198, 2021.

[6] M. Born, K. Huang, and M. Lax. Dynamical theory of crystal lattices. *American Journal of Physics*, 23(7):474–474, 1955.

[7] L. Brillouin. *Wave propagation in periodic structures: electric filters and crystal lattices*, volume 2. Dover publications, 1953.

[8] Y. Bréchet and J. Embury. Architectured materials: Expanding materials space. *Scripta Materialia*, 68(1):1–3, 2013.

[9] T. Bückmann, N. Stenger, M. Kadic, J. Kaschke, A. Frölich, T. Kennerknecht, C. Eberl, M. Thiel, and M. Wegener. Tailored 3d mechanical metamaterials made by dip-in direct-laser-writing optical lithography. *Advanced Materials*, 24(20):2710–2714, 2012.

[10] C. Combescure, P. Henry, and R. Elliott. Post-bifurcation and stability of a finitely strained hexagonal honeycomb subjected to equi-biaxial in-plane loading. *International Journal of Solids and Structures*, 88-89:296–318, 2016.
[11] V. Deshpande, M. Ashby, and N. Fleck. Foam topology: bending versus stretching dominated architectures. *Acta Materialia*, 49(6):1035–1040, 2001.

[12] V. Deshpande, N. Fleck, and M. Ashby. Effective properties of the octet-truss lattice material. *Journal of the Mechanics and Physics of Solids*, 49(8):1747–1769, 2001.

[13] R. Elliott, N. Triantafyllidis, and J. Shaw. Stability of crystalline solids—i: Continuum and atomic lattice considerations. *Journal of the Mechanics and Physics of Solids*, 54(1):161–192, 2006.

[14] L. Euler. *Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes, sive solutio problematis isoperimetrici lattissimo sensu accepti*. Marcum Michaelum Bousquet, Lausanne & Geneva, 1744.

[15] M. François. *Tenseurs en Mécanique X3PM040 Anisotropie et composites M2 Mécanique et Fiabilité des Structures*. 2019.

[16] M. François, L. Chen, and M. Coret. Elasticity and symmetry of triangular lattice materials. *International Journal of Solids and Structures*, 129:18 – 27, 2017.

[17] G. Friesecke and F. Theil. Validity and failure of the cauchy-born hypothesis in a two-dimensional mass-spring lattice. *Journal of nonlinear Science*, 12(5), 2002.

[18] G. Geymonat, S. Müller, and N. Triantafyllidis. Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity. *Archive for Rational Mechanics and Analysis*, 122(3):231–290, 1993.

[19] L. Gibson and M. Ashby. *Cellular Solids: Structure and Properties*. Cambridge Solid State Science Series. Cambridge University Press, 2 edition, 1997.

[20] B. Haghpanah, J. Papadopoulos, D. Mousanezhad, H. Nayeb-Hashemi, and A. Vaziri. Buckling of regular, chiral and hierarchical honeycombs under a general macroscopic stress state. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 470(2167):20130856, 2014.

[21] R. Hill. Elastic properties of reinforced solids: Some theoretical principles. *Journal of the Mechanics and Physics of Solids*, 11(5):357 – 372, 1963.

[22] S. Kang, S. Shan, A. Košmrlj, W. Noorduin, S. Shian, J. Weaver, D. Clarke, and K. Bertoldi. Complex ordered patterns in mechanical instability induced geometrically frustrated triangular cellular structures. *Physical Review Letters*, 112(9):098701, 2014.
[23] N. Kesmia, N. Auffray, and B. Desmorat. Towards a versatile polynomial criterion for
archituctured materials. In Proc. of the 25th Congrès Français de Mécanique. AFM,
2022.

[24] M. A. Kotob, C. Combescure, M. Mazière, T. Rose, and S. Forest. A general and efficient
multistart algorithm for the detection of loss of ellipticity in elastoplastic structures.
International Journal for Numerical Methods in Engineering, 121(5):842-866, 2020.

[25] J. Mandel. Contribution théorique à l'étude de l'écrouissage et des lois de l'écoulement plastique. In H. Görtler, editor, Applied Mechanics, pages 502–509, Berlin, Heidelberg, 1966. Springer Berlin Heidelberg.

[26] N. Ohno, D. Okumura, and T. Niikawa. Long-wave buckling of elastic square honeycombs subject to in-plane biaxial compression. International journal of mechanical sciences, 46(11):1697-1713, 2004.

[27] N. Ohno, D. Okumura, and H. Noguchi. Microscopic symmetric bifurcation condition of cellular solids based on a homogenization theory of finite deformation. Journal of the Mechanics and Physics of Solids, 50(5):1125-1153, 2002.

[28] D. Okumura, N. Ohno, and H. Noguchi. Post-buckling analysis of elastic honeycombs subject to in-plane biaxial compression. International Journal of Solids and Structures, 39(13):3487-3503, 2002.

[29] J. Overvelde, S. Shan, and K. Bertoldi. Compaction through buckling in 2d periodic, soft and porous structures: effect of pore shape. Advanced Materials, 24(17):2337-2342, 2012.

[30] S. Papka and S. Kyriakides. In-plane compressive response and crushing of honeycomb. Journal of the Mechanics and Physics of Solids, 42(10):1499-1532, 1994.

[31] S. Papka and S. Kyriakides. Experiments and full-scale numerical simulations of in-plane crushing of a honeycomb. Acta Materialia, 46(8):2765-2776, 1998.

[32] S. Papka and S. Kyriakides. Biaxial crushing of honeycombs: —part 1: Experiments.
International Journal of Solids and Structures, 36(29):4367-4396, 1999.

[33] J. Rice. The localization of plastic deformation. In *in*: W.T. Koiter (Ed.), Theoretical and Applied Mechanics, pages 207-220. North-Holland Publishing Company, 1976.
[34] M. Santisi d'Avila, N. Triantafyllidis, and G. Wen. Localization of deformation and loss of macroscopic ellipticity in microstructured solids. *Journal of the Mechanics and Physics of Solids*, 97:275-298, 2016.

[35] J. Shim, S. Shan, A. Košmrlj, S. Kang, E. Chen, J. Weaver, and K. Bertoldi. Harnessing instabilities for design of soft reconfigurable auxetic/chiral materials. *Soft Matter*, 9(34):8198-8202, 2013.

[36] F. Trentadue, D. De Tommasi, and N. Marasciuolo. Stability domain and design of a plane metamaterial made up of a periodic mesh of rods with cross-bracing cables. *Applications in Engineering Science*, 5:100036, 2021.

[37] N. Triantafyllidis and M. Schraad. Onset of failure in aluminum honeycombs under general in-plane loading. *Journal of the Mechanics and Physics of Solids*, 46(6):1089-1124, 1998.

[38] C. Truesdell and W. Noll. *The Non-Linear Field Theories of Mechanics*, pages 1-579. Springer Berlin Heidelberg, Berlin, Heidelberg, 2004.

[39] L. Truskinovsky and A. Vainchtein. Quasicontinuum modelling of short-wave instabilities in crystal lattices. *Philosophical Magazine*, 85(33-35):4055-4065, 2005.

[40] A. Vyatskikh, S. Delalande, A. Kudo, X. Zhang, C. Portela, and J. Greer. Additive manufacturing of 3D nano-architected metals. *Nature Communications*, 9(1):1-8, 2018.

[41] M. Zhang, Z. Yang, Z. Lu, B. Liao, and X. He. Effective elastic properties and initial yield surfaces of two 3D lattice structures. *International Journal of Mechanical Sciences*, 138-139:146-158, 2018.